%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