%AR Spectral Estimation Example done during Session 27. %Solve Yule-Walker Equations with autocorrelation matrix %formed from biased estimate of autocorrelation sequence clear 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; Nl=512; poles=[.9*exp(-j*.4*pi) .9*exp(j*.6*pi) .9*exp(j*.7*pi)]; arpoly=poly(poles).'; % %create AR(P) process: % x=zeros(Nl+P,1); nu=randn(Nl,1); for n=1:Nl; 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 unbiased or biased estimate of %autocorrelation sequence xrm=x(P:Nl); r=xcorr(xrm,'unbiased'); [Y,I]=max(r); % %Solve Yule-Walker Equations for estimate of AR coefficients % R=toeplitz(conj(r(I(1):I(1)+P-1,1))); arest=-inv(R)*r(I(1)+1:I(1)+P); % %plot poles % polar(0,1,'.') hold on plot(real(poles),imag(poles),'x','MarkerSize',18,'Linewidth',4); 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 truespec=abs(1./fftshift(fft(arpoly,1024))); fmag=20*log10(abs(truespec)); fmag=fmag-max(fmag); plot(freqomega,fmag,'b','LineWidth',4); freqomega=linspace(-pi,pi,2048); fmag=20*log10(abs(fftshift(fft(xrm,2048)))); fmag=fmag-max(fmag); pause plot(freqomega,fmag,'r','LineWidth',3); legend('AR spectral estimate','true spectrum','sample spectrum') hold off pause clf z=roots(arpoly); zest=roots([1; arest]); polar(0,1,'.') hold on plot(real(z),imag(z),'kx','MarkerSize',18,'Linewidth',4); plot(real(zest),imag(zest),'rx','MarkerSize',18,'Linewidth',4); hold off