% E04_5.m             28apr17   12apr22
% inverse filters and recursive deconvolution
clear all;
close all;

% ze = [.9*exp(j*.1*pi), .9*exp(-j*.1*pi), -.4]';  % minimum phase, invertible
% ze = [.9*exp(j*.1*pi), .9*exp(-j*.1*pi), -.97]'; % risky...
 ze = [.9*exp(j*.1*pi), .9*exp(-j*.1*pi), -1.41]'; % noninvertible

%ze = [1*exp(j*.1*pi), 1*exp(-j*.1*pi), 0]'; % [inv.syst is a sort of oscillator]
% ze = [1*exp(j*.55*pi), 1*exp(-j*.55*pi), 0]'; % [...even at a different freq.]
                          % [oscillates only if triggered by added input noise]
po = [.7+j*.6, .7-j*.6, 0]';
[b,a] = zp2tf(ze,po,1);

% alpha=.8;  b=[1-alpha]; a=[1, -alpha]; % simple IIR lowpass filter
   % e.g. a temperature transducer having finite thermal capacity
   % inv. filter is an FIR, h = {1, -alpha}/(1-alpha)
   
% b = ones(1,7)/7;  a = 1;  % moving average filter (noninvertible)

% fvtool(b,a);  fvtool(a,b); return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(5);
zplane(ze,po)

n = [0:1999];
x =  sin(.1*pi*n) + cos(.77*pi*n) + randn(size(n))*1;
w = filter(b,a,x);      % (.1*pi) component is strongly attenuated

% w(53) = w(53) + .1;
% w = w + randn(size(n))*.01;  % distortion includes noise: w = h*x + n
% b = round(b*1e3) /1e3;  a = round(a*1e3) /1e3;  % imprecise model

y = filter(a,b,w);      %  y=x if inverse system is stable
   %%%%%%%%%%%% note this script masks that y is a DELAYED version of x
figure(1);
subplot(411); plot(n,x); axis([0,length(n)-1,-2,2]);
subplot(412); plot(n,w); axis([0,length(n)-1,-2,2]);
subplot(413); plot(n,y); axis([0,length(n)-1,-2,2]);
subplot(414); plot(n,y-x); 
figure(2);
subplot(311); X = abs(fft(x,256)); plot(X(1:129));
subplot(312); W = abs(fft(w,256)); plot(W(1:129));
subplot(313); Y = abs(fft(y,256)); plot(Y(1:129));
return;


%% recursive deconvolution

[x,fs] = audioread('chimes.wav');
%fs = fs / 2;   
x = x(:,1)';   % left channel
L = fs*4;     % L <= length(chimes)=13921;
x = x(1:L);
x = x / max(x);   % amplify to get max=1
X = fft(x,16384);
figure(11);   plot(abs(X(1:8193)));
%soundsc(x,fs);    pause(3); 

 ze = [.9*exp(j*.1*pi), .9*exp(-j*.1*pi), -.99]';  % minimum phase
% ze = [.9*exp(j*.1*pi), .9*exp(-j*.1*pi), -1.01]'; % noninvertible
 po = [.7+j*.7, .7-j*.7, 0]';
 [b,a] = zp2tf(ze,po,1);
 
 b = [1 zeros(1,2500) -1]/2; a = 1;
% fvtool(a,b); return; % do not try to find the poles!

h = impz(b,a,L)';  
% figure(22);  plot(h); return;
w = conv(h,x);       %  note: length(w)=L+L-1
% w = w + randn(size(w))*.03;  % noise added: w = h*x + n

W = fft(w,16384);
figure(22); plot(abs(W(1:8193)));
 soundsc(w,fs);  pause(2.5); 

% h = round(h*1e2)/1e2;  % rounding error in coeff. estimate (try 1e2 - 1e5)

y = zeros(size(x));
y(1) = w(1) / h(1);
for k=2:L
    y(k) = ( w(k) - fliplr(h(2:k)) * y(1:k-1)' ) / h(1);
%    if (k==5000); y(k) = round(y(k)*10)/10; end;     % rounding in the recursion
%    if (k>2000); y(k) = round(y(k)*10)/10; end;    
end;

figure(33);
 subplot(311); plot(x); ylabel('x');
 subplot(312); plot(w(1:L)); ylabel('w(1:L)');
 subplot(313); plot(y-x); ylabel('error');

 soundsc(y,fs);

   %   note: the Matlab command
% [y,r] = deconv(w,h); 
   %   calculates
   %   y = filter(w,h,[1 zeros(1,L-1)]);
   %         (which is formally identical to rec. deconv.) and
   %   r = w-conv(h,y); % the remainder