% This test file compares fft-based estimation
% with spectral envelopes and subspace methods.


% Setting up signal
N=250; mag0=2; mag1=1.5; o1=1.3;
               mag2=2; o2=1.35;
               t=0:N-1; t=t(:); y=mag0*randn(1,N)'.*exp(j*randn(1,N)');
               f1=2*pi*rand; y=y+mag1*exp(j*(o1*t+f1));
               f2=2*pi*rand; y=y+mag2*exp(j*(o2*t+f2));


% plotting true and fft spectra
figure(1), clf, hold on, th=0:N-1; th=th*2*pi/N;
               plot(th,abs(fft(y.*hamming(N)))/N,'g')
               plot(th,abs(fft(y))/N,'g-.')
               arrow([o1 o2],[mag1 mag2])
               plot( [0 2*pi],[mag0 mag0],'r')


% setting up filter parameters and plotting response
thetamid=1.325; [Ah,bh]=cjordan([5],[0.9*exp(thetamid*j)]);
              [sv,w]=rsigma(Ah,bh,linspace(0,2*pi,256));
              Mag=max([mag0 mag1 mag2 max(sv)])*1.2; %axis([0 2*pi 0 Mag])
              plot(w(:),sv(:),'r--')

% obtaining and plotting spectral envelope
[sqrtrho,theta]=env(Ah,bh,y);
              plot(theta,sqrtrho.^2,'b'),
              xlabel('blue:envelope, green:fft, green--:filter, red:true')

% obtaining and plotting subspace analysis results
y=y(:); y=y-mean(y);
              [x ,t]=rlsim(Ah,bh,y); Pf=cov(conj(x));
              [xr,t]=rlsim(Ah,bh,conj(flipud(y))); Pr=cov(conj(xr));
[Wf1,wpolyf1,Pf1,errorf]=commutant(Pf,Ah,eye(size(Ah)));
[Wr1,wpolyr1,Pr1,errorb]=commutant(Pr,Ah,eye(size(Ah)));
if errorf<errorb, P=Pf1; errorf, else, P=Pr1; errorb, end
[omega_ss,residues_ss]=sm(P,Ah,bh);
arrowb(2*pi-omega_ss,residues_ss)

%print -depsc mag9whole
%axis([1.1 1.7 0 14])
%print -depsc mag9
