%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


