%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