%Minimum Variance Spectral Estimation for Sum of
%Sinusoids Example done during Session 30.
%Forward-Backward Averaged Correlation Matrix Employed.
clear all
clf
set(0,'defaultaxesfontsize',20);
P=3;
Nl=128;
M=10;
%
%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 .3*pi];
poles=exp(j*omegas);
amps=[.5 1 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;
%
%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);
%
%plot poles and zeroes
%
polar(0,1,'.')
hold on
plot(real(poles),imag(poles),'x','MarkerSize',16,'Linewidth',3);
hold off
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]).';
specest(n)=1/(s'*Rinv*s);
end
fmag=10*log10(abs(specest));
fmag=fmag-max(fmag);
plot(freqomega,fmag,'k','LineWidth',4);
axis([-pi pi -60 0]);
title('Min Var Spectral Estimate'),...
xlabel('Omega'), ylabel('Magnitude'), grid
pause
hold
specfft=20*log10(abs(fftshift(fft(x,1024))));
fmag1=specfft-max(specfft);
freqomega=linspace(-pi,pi,1024);
plot(freqomega,fmag1,'r','LineWidth',2);
pause
for p=1:P;
plot([omegas(1,p) omegas(1,p)],[-60 0],'b','LineWidth',3);
end
legend('Min Var Estimate','Mag. Squared of FFT');
hold off
for p=1:P+1;
omega0=input('Input frequency of interest:')
a0=exp(j*omega0*[0:M-1]).';
h0=(Rinv*a0)/(a0'*Rinv*a0);
H0=20*log10(abs(fftshift(fft(h0,1024))));
H0mag=H0-max(H0);
plot(freqomega,H0mag,'k','LineWidth',4);
axis([-pi pi -50 10]);
title('DTFT of h[n;omega0]'),...
xlabel('Omega'), ylabel('Magnitude'), grid
hold on
for p=1:P;
plot([omegas(1,p) omegas(1,p)],[-50 10],'b','LineWidth',3);
end
hold off
pause
y=real(conv(x,h0));
plot(y,'LineWidth',4);
axis([0 (Nl+M-1) -1 1]);
title('output when x[n] is input to h(n;omega0)'),...
xlabel('Time'), ylabel('Amplitude'), grid
end

