%Estimation of Power Spectrum of ARMA process via %Yule-Walker based AR Spectral Estimation %done during Session 27. 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; poles=[.9*exp(-j*.4*pi) .9*exp(j*.6*pi) .9*exp(j*.8*pi)]; zeeroes=[.95*exp(j*.25*pi) .95*exp(-j*.7*pi)]; Q=2; arpoly=poly(poles).'; arnum=poly(zeeroes).'; PQ=12; Nl=512; % %create ARMA(P,Q) process: % x=zeros(Nl+P,1); nu=randn(Nl+Q+P,1); for n=1:Nl; for q=1:Q+1; x(n+P,1)=x(n+P,1)+arnum(q)*nu(n+Q+1-q); end 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 % xrm=x(P:Nl); r=conj(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',18,'Linewidth',3); plot(real(zeeroes),imag(zeeroes),'o','MarkerSize',18,'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,'b','LineWidth',4); axis([-pi pi -60 0]); title('AR Spectral Estimate'),... xlabel('Omega'), ylabel('Magnitude'), grid hold truespec=abs(fftshift(fft(arnum,1024))./fftshift(fft(arpoly,1024))); fmag=20*log10(abs(truespec)); fmag=fmag-max(fmag); plot(freqomega,fmag,'k','LineWidth',4); %freqomega=linspace(-pi,pi,2048); %fmag=20*log10(abs(fftshift(fft(xrm,2048)))); %fmag=fmag-max(fmag); %plot(freqomega,fmag,'r','LineWidth',3); legend('AR spectral estimate','true 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); plot(real(zeeroes),imag(zeeroes),'ko','MarkerSize',18,'Linewidth',3); hold off