%Minimum Variance Spectral Estimation for ARMA(p,q) %Process Example done during Session 30. %Forward-Backward Averaged Correlation Matrix Employed. clear all clf set(0,'defaultaxesfontsize',20); %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; PQ=16; Nl=512; Itilde=zeros(PQ+1,PQ+1); for n=1:PQ+1; Itilde(n,PQ-n+2)=1; end poles=[.9*exp(-j*.7*pi) .95*exp(j*.3*pi) .95*exp(j*.75*pi)]; zeeroes=[.75*exp(j*.15*pi) .75*exp(-j*.9*pi)]; Q=2; arpoly=poly(poles).'; arnum=poly(zeeroes).'; % %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 unconstrained LS estimate of forward- %backward averaged correlation matrix -- NOT Toeplitz % Rf=zeros(PQ+1,PQ+1); xr=x(P:Nl,1); for n=1:Nl-2*PQ+1; Rf=Rf+xr(n:n+PQ,1)*xr(n:n+PQ,1)'; end Rf=Rf/(Nl-PQ); Rbb=Rf(PQ+1:-1:1,:); Rb=conj(Rbb(:,PQ+1:-1:1)); Rfb=0.5*(Rf+Rb); Rinv=inv(Rfb); % %plot poles and zeroes % 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 MV estimate of spectrum, plot and compare %with true spectrum % freqomega=linspace(-pi,pi,1024); for n=1:1024; a=exp(-j*freqomega(n)*[0:PQ]); specest(n)=1/(a*Rinv*a'); end fmag1=10*log10(abs(specest)).'; truespec=abs(fftshift(fft(arnum,1024))./fftshift(fft(arpoly,1024))); fmag2=20*log10(abs(truespec)); offset=mean(fmag2-fmag1); fmag1=fmag1+offset; fmax=max([fmag1 ; fmag2]); fmag1=fmag1-fmax; fmag2=fmag2-fmax; plot(freqomega,fmag1,'k','LineWidth',4); axis([-pi pi -60 0]); title('Min Var Spectral Estimate'),... xlabel('Omega'), ylabel('Magnitude'), grid hold plot(freqomega,fmag2,'b','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('MinVar spectral estimate','true spectrum') hold off