%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

