%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