%AR Spectral Estimation Example done during Session 29.
%Yule-Walker Method with autocorrelation matrix
%formed from biased estimate of autocorrelation sequence versus
%Unconstrained LS Method with forward-backward averaging.
clear all
clf
set(0,'defaultaxesfontsize',22);
%P=input('Number of poles of AR process');
%poles=input('Input P pole locations');
%M=input('Number of zeroes');
%zeroes=input('Input M zero locations');
%P=input('length of AR sequence');
P=3;  Phat=3
Nl=256;
poles=[.9*exp(-j*.4*pi) .95*exp(j*.6*pi) .95*exp(j*.7*pi)];
arpoly=poly(poles).';
%
%create AR(P) process:
%
x=zeros(Nl+P,1);
nu=randn(Nl+P,1);
for n=1:Nl+P;
x(n+P,1)=nu(n);
for k=1:P;
x(n+P,1)=x(n+P,1)-arpoly(k+1)*x(n+P-k);
end
end
%compute biased estimate of autocorrelation sequence
%
xr=x(2*P+1:Nl+P);
br=xcorr(xr,'biased');
[Y,I]=max(br);
%
%Solve Yule-Walker Equations for estimate of AR coefficients
%
Rc=toeplitz(conj(br(I(1):I(1)+Phat-1,1)));
arYW=-inv(Rc)*br(I(1)+1:I(1)+Phat);
%
%compute unconstrained LS estimate of 
%correlation matrix -- NOT Toeplitz
%
Rf=zeros(Phat,Phat);
rf=zeros(Phat,1);
for l=1:Phat;
for k=1:l;
for n=2:Nl-2*Phat+1;
Rf(l,k) = Rf(l,k) + xr(Phat+n-k)*conj(xr(Phat+n-l));
end
Rf(k,l)=conj(Rf(l,k));
end
for n=1:Nl-2*Phat;
rf(l,1) = rf(l,1) + xr(Phat+n)*conj(xr(Phat+n-l));
end
end
arLS=-inv(Rf)*rf;
%
%compute estimate of spectrum, plot and compare
%with true spectrum
%
specestYW=1./abs(fftshift(fft([1 ; arYW],1024))).^2;
specestLS=1./abs(fftshift(fft([1 ; arLS],1024))).^2;
truespec=1./abs(fftshift(fft(arpoly,1024))).^2;
freqomega=linspace(-pi,pi,1024);
fmagYW=10*log10(specestYW);
fmagLS=10*log10(specestLS);
fmag=10*log10(truespec);
Fmax=max(fmag);
fmag=fmag-Fmax;
fmagYW=fmagYW-Fmax;
fmagLS=fmagLS-Fmax;
plot(freqomega,fmagYW,'k','LineWidth',4);
axis([-pi pi -40 10]);
title('AR Spectral Estimates'),...
xlabel('Omega'), ylabel('Magnitude'), grid
hold
plot(freqomega,fmagLS,'b','LineWidth',4);
plot(freqomega,fmag,'r','LineWidth',4);
legend('Yule-Walker','Unconstrained LS','True Spectrum')
hold off
diffYW=real(specestYW-truespec)/1024;
ErrorYW=mean(diffYW.*diffYW)
specestLS=1./abs(fftshift(fft([1 ; arLS],1024))).^2;
diffLS=real(specestLS-truespec)/1024;
ErrorLS=mean(diffLS.*diffLS)
pause
clf
z=roots(arpoly);
zYW=roots([1; arYW]);
zLS=roots([1; arLS]);
polar(0,1,'.')
hold on
plot(real(z),imag(z),'rx','MarkerSize',18,'Linewidth',4);
plot(real(zYW),imag(zYW),'kx','MarkerSize',18,'Linewidth',4);
plot(real(zLS),imag(zLS),'bx','MarkerSize',18,'Linewidth',4);
hold off

