%Minimum Variance Spectral Estimation for Sum of
%Sinusoids vs MUSIC Example done during Session 10.
%Forward-Backward Averaged Correlation Matrix Employed.
clear all
clf
set(0,'defaultaxesfontsize',20);
P=3;
Nl=64;
M=12;
%
%P=input('no. of complex sinewaves = ');
%omegas=input('frequencies of the sinewaves = ');
%amps=input('amplitudes of the sinewaves = ');
%omegas=[.55*pi .6*pi .65*pi];
omegas=[-.2*pi .1*pi .2*pi];
poles=exp(j*omegas);
amps=[.5 1 1];
noisepwr=1;
x=zeros(Nl,1);
%
%create Nl values of SOS data
%
for i=1:P;
x(1:Nl,1)=x(1:Nl,1)+amps(1,i)*exp(j*omegas(1,i)*(0:Nl-1)');
end;
x=x+sqrt(noisepwr)*randn(size(x));
%
%compute unconstrained LS estimate of autocorrelation matrix
%
Rb=zeros(M,M);
for n=1:Nl-M+1;
Rb=Rb+x(n:n+M-1,1)*x(n:n+M-1,1)';
end
Rff=Rb(M:-1:1,:);  Rf=conj(Rff(:,M:-1:1));
Rfb=0.5*(Rf+Rb)/(Nl-M+1)+.001*eye(M);
Rinv=inv(Rfb);
%Eigenanalysis
[E,D,V]=svd(Rfb);
Es=E(:,1:P);
Pn=E(:,P+1:M)*E(:,P+1:M)';
%
%plot poles and zeroes
%
%ESPRIT algorithm:
Es1=Es(1:M-1,:); Es2=Es(2:M,:);
Psi=Es1\Es2;
[E,Phi]=eig(Psi);
Phivec=diag(Phi);
%plot eigenvalues from ESPRIT and compare with true frequencies
polar(0,1,'.')
hold on
plot(real(poles),imag(poles),'x','MarkerSize',16,'Linewidth',3);
plot(real(Phivec),imag(Phivec),'o','MarkerSize',16,'Linewidth',3);
hold off
freqests=sort(angle(Phivec))
omegas
pause
% 
%compute estimate of spectrum, plot and compare
%with true spectrum
%
freqomega=linspace(-pi,pi,1024);
for n=1:1024;
s=exp(j*freqomega(n)*[0:M-1]).';
MinVar(n)=1/(s'*Rinv*s);
MUSIC(n)=1/(s'*Pn*s);
end
fmag=10*log10(abs(MinVar));
fmag2=10*log10(abs(MUSIC));
fmag=fmag-max(fmag);
fmag2=fmag2-max(fmag2);
plot(freqomega,fmag,'k','LineWidth',4);
axis([-pi pi -60 0]);
title('Minimum Variance Vs. MUSIC'),...
xlabel('Omega'), ylabel('Magnitude'), grid
hold on
plot(freqomega,fmag2,'b','LineWidth',4);
%specfft=20*log10(abs(fftshift(fft(x,1024))));
%fmag1=specfft-max(specfft);
%freqomega=linspace(-pi,pi,1024);
%plot(freqomega,fmag1,'r','LineWidth',2);
%pause
freqest=angle(Phivec);
for p=1:P;
plot([omegas(1,p) omegas(1,p)],[-60 0],'r','LineWidth',2);
plot([freqest(p) freqest(p)],[-60 0],'m--','LineWidth',3);
end
legend('Min Var Estimate','MUSIC Estimate');
hold off