%AR Spectral Estimation Example done during Session 27.
%Solve Yule-Walker Equations with autocorrelation matrix
%formed from biased estimate of autocorrelation sequence
clear all
clf
set(0,'defaultaxesfontsize',20);
PQ=3;
Nl=512;
%
%P=input('no. of complex sinewaves = ');
%omegas=input('frequencies of the sinewaves = ');
%amps=input('amplitudes of the sinewaves = ');
P=3;
omegas=[-.4*pi .6*pi .7*pi];
poles=exp(j*omegas);
amps=[.25 1 1];
x=zeros(Nl,1);
%
%create Nl values of SOS data
%
for i=1:P;
x(1:Nl,1)=x(1:Nl,1)+amps(1,i)*exp(j*omegas(1,i)*(0:Nl-1)');
end;
%
%compute biased estimate of autocorrelation sequence
%
xrm=x(P:Nl);
r=xcorr(xrm,'biased');
[Y,I]=max(r);
%
%Solve Yule-Walker Equations for estimate of AR coefficients
%
Rc=toeplitz(conj(r(I(1):I(1)+PQ-1,1)));
arest=-inv(Rc)*r(I(1)+1:I(1)+PQ);
%
%plot poles and zeroes
%
polar(0,1,'.')
hold on
plot(real(poles),imag(poles),'x','MarkerSize',16,'Linewidth',3);
hold off
pause
% 
%compute estimate of spectrum, plot and compare
%with true spectrum
%
specest=abs(1./fftshift(fft([1 ; arest],1024)));
freqomega=linspace(-pi,pi,1024);
fmag=20*log10(abs(specest));
fmag=fmag-max(fmag);
plot(freqomega,fmag,'k','LineWidth',4);
axis([-pi pi -60 0]);
title('AR Spectral Estimate'),...
xlabel('Omega'), ylabel('Magnitude'), grid
hold
for p=1:P;
plot([omegas(1,p) omegas(1,p)],[-60 0],'b','LineWidth',3);
end
hold off
pause
clf
zest=roots([1; arest]);
polar(0,1,'.')
hold on
plot(real(poles),imag(poles),'kx','MarkerSize',18,'Linewidth',4);
plot(real(zest),imag(zest),'rx','MarkerSize',18,'Linewidth',4);
hold off

