echo on
format short
format compact
clc
pkg load signal
graphics_toolkit("gnuplot");
% Load SX133.wav
[X,FS]=audioread("SX133.wav");
Nyquist = FS/2;
% Extract samples 40234-40444
X1=X(40234:40444);
% Make a time-domain plot
Time=[40234:40444]/FS;
figure(1); plot(Time,X1,'b');
pause
clc
% Apply Hamming window to X1
X1a=X1.*hamming(length(X1));
% Plot the windowed audio samples
figure(2); plot(Time,X1a);
pause
clc
%
% 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;
pause
% Now calculate the log power spectrum, and plot it:
F1=fft(X2);
F2=20*log10(abs(F1(1:nfft2)));
F2=F2-max(F2);
Freqs=(FS/2)*(0:(nfft2-1))/nfft2;
figure(3); plot(Freqs,F2,'r-');
axis("auto");
xlabel("Frequency"); ylabel("Power (dB)");
pause
clc
% Now calculate the (Nth-order) auto-regression coefficients
%     for the original sample.
pause
%
% We'll do this by solving for x in the equation
%
%  A*x = b
%
%  where b is the sequence we're trying to predict
%        x is the set of auto-regression coefficients doing the prediction
%   and  A is a matrix where each row has N samples
%               used to predict the next sample
%
% We know b, we can use b to create A,
%    and then we can solve for x.
pause
clc
%  So first, write a simple function for setting up an A matrix
%
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
%
pause
% To see what this is doing, consider
makeA(1:10,3)
% Look up "hankel" if you want the math -- otherwise just accept the magic...
%   (...and you can ignore the warnings)
pause
%
clc
% So in this case, (remember that) our input is X1.
%    We'll try an order-14 prediction (up to 7 complex poles...)
Order=14;
%
% First we remove the mean
%      (because our regression equation has a zero intercept,
%        and we're not interested in any additive term...)
%
X3=X1-mean(X1);
%
pause
% Now, to solve for x in A*x=b
%     we'll multiply both sides of the equation by the (pseudo-)inverse of A.
%
%          and we create the A matrix from our mean-corrected samples X3:
%
A=makeA(X3,Order);
%
pause
% b (in Ax=b) is X3 starting from the sample number Order+1
b = X3((Order+1):length(X3));
pause
% ...so now to  solve for x, left-multiplying both sides by 1/A:
x = pinv(A)*b;
%
%  We should probably check that A is full rank
% 
rank(A)
%
pause
clc
% So what are the resulting predictor coefficients?
x
%
pause
%  Not easy to interpret in that form...
%
pause
%
% So to start with,
%  Let's ask how good a job our auto-regressive prediction is doing.
%     Let's plot the time-domain version and see:
figure(4); plot(A*x,'bx',b,'ro');
%
% (Original values are red o's, predicted values are blue x's...)
pause
% Not bad!  What's the r^2 (percent of variance accounted for)?
r = corr(A*x,b); r^2
%
% Again, pretty good!
pause
clc
% Now let's look at the impulse response of the recursive filter
%     whose coefficients are x ---
% Well, OK, whose coefficients are 1 followed by -1 times x backwards...
%   we have to obey the conventions derived from the LCDE formulation
%  ...see https://www.ling.upenn.edu/courses/ling525/NewFilters1.html
%
bcoef = zeros(Order+1); bcoef(1)=1;
Imp = zeros(512,1); Imp(Order)=1;
acoef = [1 -x(Order:-1:1)'];
ImpResp = filter(bcoef, acoef, Imp);
figure(5)
plot(ImpResp);
pause
% And its spectrum...
figure(6)
FImpResp = 20*log10(abs(fft(ImpResp)(1:nfft2)));
FImpResp = FImpResp-max(FImpResp);
plot(Freqs,FImpResp)
axis("auto");
xlabel("Frequency"); ylabel("Power (dB)");
%
pause
% And now let's add the spectrum of the original signal --
%
hold on; plot(Freqs,F2); hold off
pause
% Again, not bad, given that we're approximating 256 numbers by 14 numbers...
%
pause
clc
% A 14th-order analysis yields 7 complex poles
%  (though we might have gotten a couple of real poles in place of one of them)
% The roots:
roots(acoef)
pause
% Now the pole frequencies and amplitudes:
PFreqs = (angle(roots(acoef))*Nyquist/(pi))'
PAmps = abs(roots(acoef))'
% Let's make a table of of the seven pole frequencies and amplitudes,
%       sorted by frequency:
Poles = sortrows([F=PFreqs(1:2:Order); PAmps(1:2:Order)]')
% ...or in easier-to-read form:
printf("%.0f %.3f\n",Poles')
%
% That's the first step in modern "formant trackers"
%   which identify some of the LPC poles with "formants".
%
%   Though you should keep in mind that the LPC anaysis
%                           %     is NOT "formant analysis",
%     but rather a technique for approximating a signal in the time domains,
%        using a radically reduced number of values.
%
%     And in the freqency domain,
%         the analysis is modeling not only "formants"
%     but also the source spectrum,
%     the room acoustics,
%     the microhone frequency response,
%     any background noise,
%     etc. etc.
