%ARMA Spectral Estimation Example done during Session 29. %Uses two-step estimation procedure developed in %Section 12.3.8 of the Proakis and Manolakis text %that averts solution of nonlinear system of equations. 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=8; Nl=512; % %create ARMA(P,Q) process: % x=zeros(Nl+3*P+1,1); y=zeros(Nl+3*P+1,1); nu=randn(Nl+3*P,1); for n=1:Nl+2*P; 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(2*P+1:Nl+2*P); %br=conj(xcorr(xrm,'biased')); br=xcorr(xrm,'biased'); [Y,I]=max(br); % %Solve Yule-Walker equations using lag values greater %than Q on the right hand side of the equation to estimate %the AR coefficients. I created an overdetermined set of %PQ=8 equations in P=3 unknowns. % imax=I(1)+Q; Rc=toeplitz(br(imax:imax+PQ-1,1),br(imax:-1:imax-P+1,1)); arest=-Rc\br(imax+1:imax+PQ); % %plot poles and zeroes % MA=conv(xrm,[1; arest]); mr=xcorr(MA(Q+1:Nl),'biased'); [Rm,I]=max(mr); rmm=mr(I(1)-Q:I(1)+Q); % polar(0,1,'.') hold on plot(real(poles),imag(poles),'x','MarkerSize',16,'Linewidth',3); plot(real(zeeroes),imag(zeeroes),'o','MarkerSize',16,'Linewidth',3); hold off pause % %compute estimate of spectrum, plot and compare %with true spectrum % specest=abs(fftshift(fft(rmm,1024))./abs(fftshift(fft([1 ; arest],1024))).^2); freqomega=linspace(-pi,pi,1024); fmag=10*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('ARMA 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