echo on
format shortG
format compact
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Same as Homework3a, but using arburg() to calculate regression coefficients
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pause
clc
% The next two lines are for Octave:
% pkg load signal
% graphics_toolkit("gnuplot");
%
% Load SX133.wav
[X,FS]=audioread("SX133.wav");
Nyquist = FS/2;
% Extract samples 40234-40444
%   = 2.514 to  2.527 seconds
%   middle of the syllable "quick"
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, using the arburg() function.
pause
%    We'll try an order-14 prediction (up to 7 complex poles...)
Order=14;
arcoeffs = arburg(X1,Order)
%
pause
%  ...and let's compare the spectra again:
%
pause
bcoeffs = zeros(1,Order+1); bcoeffs(1)=1;
Imp = zeros(512,1); Imp(Order)=1;
ImpResp = filter(bcoeffs, arcoeffs, Imp);
figure(1)
FImpResp = abs(fft(ImpResp));
FImpResp = 20*log10(FImpResp(1:nfft2));
FImpResp = FImpResp-max(FImpResp);
plot(Freqs,FImpResp)
axis("auto");
xlabel("Frequency"); ylabel("Power (dB)");
title(sprintf("%dth order filter",Order));
hold on; plot(Freqs,F2); hold off
pause
% Roots, by factoring the polynomial coefficients
Roots = roots(arcoeffs)
pause
% The pole frequencies and bandwidths:
PFreqs = (angle(Roots)*Nyquist/(pi))'
PAmps = abs(Roots)'
Poles14 = sortrows([PFreqs(1:2:Order); PAmps(1:2:Order)]')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Same deal with a higher-order filter:
Order=18;
arcoeffs = arburg(X1,Order)
bcoeffs = zeros(1,Order+1); bcoeffs(1)=1;
Imp = zeros(512,1); Imp(Order)=1;
ImpResp = filter(bcoeffs, arcoeffs, Imp);
figure(2)
FImpResp = abs(fft(ImpResp));
FImpResp = 20*log10(FImpResp(1:nfft2));
FImpResp = FImpResp-max(FImpResp);
plot(Freqs,FImpResp)
axis("auto");
xlabel("Frequency"); ylabel("Power (dB)");
title(sprintf("%dth order filter",Order));
hold on; plot(Freqs,F2); hold off
pause
% Roots, by factoring the polynomial coefficients
Roots = roots(arcoeffs)
pause
% The pole frequencies and bandwidths:
PFreqs = (angle(Roots)*Nyquist/(pi))'
PAmps = abs(Roots)'
Poles18 = sortrows([PFreqs(1:2:Order); PAmps(1:2:Order)]')
pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% And with a lower-order filter:
Order=8;
arcoeffs = arburg(X1,Order)
bcoeffs = zeros(1,Order+1); bcoeffs(1)=1;
Imp = zeros(512,1); Imp(Order)=1;
ImpResp = filter(bcoeffs, arcoeffs, Imp);
figure(3)
FImpResp = abs(fft(ImpResp));
FImpResp = 20*log10(FImpResp(1:nfft2));
FImpResp = FImpResp-max(FImpResp);
plot(Freqs,FImpResp)
axis("auto");
xlabel("Frequency"); ylabel("Power (dB)");
title(sprintf("%dth order filter"),Order);
hold on; plot(Freqs,F2); hold off
pause
% Roots, by factoring the polynomial coefficients
Roots = roots(arcoeffs)
pause
% The pole frequencies and bandwidths:
PFreqs = (angle(Roots)*Nyquist/(pi))'
PAmps = abs(Roots)'
Poles8 = sortrows([PFreqs(1:2:Order); PAmps(1:2:Order)]')
