echo on
format short
format compact
clc
pkg load signal
graphics_toolkit("gnuplot");
%
% Inverse filtering -- continuation of nHomework3a.m ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pause
clc
% Suppose we have a causal LTI system S1,
%     with impulse response h,
%    whose z transform, H(z), is a ratio of polynomials
%
%                       B_poly/A_poly
%
% Now suppose we have another causal LTI system S2,
%      with impulse response g,
%     whose z transform, G(z), happens to be exactly the inverse, namely
%
%                       A_poly/B_poly
%
% If we convolve some input x with h, and then convolve the results with g
%
%                       ( x # h ) # g
%
%    the associative property of convolution
%        means that this is equivalent to
%
%                       x # ( h # g )
%
pause
clc
% More interestingly, if we already have the output S1(x) = x # h
%      we can get the original input back via
%
%      x = S2(S1(x))
%
%  This is called "inverse filtering".
%
pause
% Let's try this with a simple one-pole recursive filter,
%       with a center frequency half the sampling rate,
%       excited by an impulse:
%
pause
r = .95*exp(i*.5*pi);
A= poly([r r']);
impulse=zeros(256,1); impulse(15)=1;
B = [1 0 0];
out1 = filter(B,A,impulse);
%
% (The referenced lecture notes should help you understand
%      what all that is doing --
%  for now, just consider it to be some useful magical incantations...)
%
out1(1:30)'
%
pause
clc
% Now call filter() again, with A and B reversed...
out2 = filter(A,B,out1);
out2(1:30)'
%
pause
% Success! we got our impulse back...
%
pause
clc
%
% Now let's try the recusive filter estimated by the kind of LPC modeling
%    we did in nHomework3a.m
%
% We could apply this filter to a known input, say an impulse,
% and then inverse filter, to get the input back just as we did above.
%
% But that's kind of boring --
%  it's more interesting to inverse filter the original speech.
%
% Here we don't actually know what the excitation was.
% Of course, the signal wasn't generated by this type of filter at all,
%    but rather by a physical process that is something like --
%    or at least can be usefully modeled as -- such a filter.
% So the result of this inverse filtering will be a hypothetical signal,
%    what the input would have been
%         if the output really has been created by our modeled filter.
%
pause
clc
% Do the FFT and LPC analysess as before...
% Load SX133.wav
[X,FS]=audioread("SX133.wav");
Nyquist = FS/2;
% Extract samples 40234-40444
X1=X(40234:40444);
% Apply Hamming window X1
X1a=X1.*hamming(length(X1));
% Put windowed signal in the start of a 512-element vector
nfft=512; nfft2=nfft/2;
X2=zeros(nfft,1); X2(1:length(X1a))=X1a;
% Now calculate the log power spectrum:
F1=fft(X2);
F2=20*log10(abs(F1(1:nfft2)));
F2=F2-max(F2);
Freqs=(FS/2)*(0:(nfft2-1))/nfft2;
%
pause
function A = makeA(S,N)
 % S is a vector of length M whose ith sample will be predicted
 % as a linear combination of the N previous samples.
 % Thus we want to set up a matrix A,
 % with M-N rows and N columns, 
 % whose jth row consists of the samples from S(j) to S(j+N-1).
 % Obviously M must be greater than N.
 AA = hankel(S,S(1:N));
 A = AA(1:length(S)-N,:);
endfunction
%
Order=14;
X3=X1-mean(X1);
%
A=makeA(X3,Order);
%
pause
% b (in Ax=b) is X3 starting from the sample number Order+1
b = X3((Order+1):length(X3));
% ...so now to  solve for x, left-multiplying both sides by 1/A:
x = pinv(A)*b;
%
bcoef = zeros(Order+1); bcoef(1)=1;
acoef = [1 -x(Order:-1:1)'];
%
pause
clc
% OK, now we have the a and b coefficients for filter() --
%    and to inverse-filter, we just reverse their order in the call.
%
X4 = filter(acoef, bcoef, X3);
%
pause
clc
% Now calculate the log power spectrum of the inverse-filtered signal:
X4a = X4.*hamming(length(X4));
%
X4b=zeros(nfft,1); X4b(1:length(X4a))=X4a;
FX4b=fft(X4b);
FX4c=20*log10(abs(FX4b(1:nfft2)));
FX4c=FX4c-max(FX4c);
%
pause
% And compare it to the original:
figure(1)
plot(Freqs,F2,'r-',Freqs,FX4c,'b-');
axis("auto");
xlabel("Frequency"); ylabel("Power (dB)");
title("Original specrum in red, Inverse-filtered spectrum in blue")
%
pause
% What's the difference, and why?
%
