%  SIGNAL = sinusoids + noise
%  Setting up the signal parameters and time history
N=100; n=15;
a0=1; a1=1.5; theta1=1.3; a2=2; theta2=1.35;
k=1:N; k=k(:);
y=a1*exp(1i*(theta1*k+2*pi*rand))+a2*exp(1i*(theta2*k+2*pi*rand));
% MA coefficient vector
b=[1.0000    0.4318   -0.4496   -0.5440];
u_real=randn(1,N+length(b)-1);
u_imag=randn(1,N+length(b)-1);
u=a0*(u_real+i*u_imag)/sqrt(2);
% generate noise
noise=b*toeplitz(u(length(b):-1:1),u(length(b):end));
y=y'+noise;
Ac=compan(eye(1,n+1));
Bc=eye(n,1);
That=dlsim_complex(Ac,Bc,y);
NN=2048; th=linspace(0,2*pi,NN);
T=ToepProj(That);
[Ts, Tma]=MA_decomp(T,3);
figure(1);clf
h=freqresp(f,th);
        subplot(3,1,1);
        plot(theta1,a1,'r^','MarkerSize',12);hold on;
        plot(th,vec(abs(h)).^2,'b','LineWidth',1.2);
        legend('true spectral line','true MA spectrum','Location','North');
        arrow([theta1 theta2],[a1 a2],12);
        set(gca,'xlim',[0 2*pi]);
        set(gca,'YTick',[]);


r_ma=[fliplr(Tma(1,:)) Tma(1,2:end)];
p_ma=abs(fft(r_ma,length(th)));
[w_s,r_s]=sm(Ts,Ac,Bc,rank(Ts,0.001));
r_s=r_s(w_s<pi); w_s=w_s(w_s<pi);
subplot(3,1,2);
        plot(w_s(1),r_s(1),'r^','MarkerSize',12);hold on;
        plot(th,p_ma,'b','LineWidth',1.2);
        legend('estimated spectral line','estimated MA spectrum','Location','North');
        arrow(w_s,r_s,12);
        set(gca,'xlim',[0 2*pi]);
        set(gca,'YTick',[]);


[Ts0, Tn]=MA_decomp(T,0);
[w_s0,r_s0]=sm(Ts0,Ac,Bc,rank(Ts0,0.001));
r_s0=r_s0(w_s0<pi); w_s0=w_s0(w_s0<pi);
subplot(3,1,3);
        plot(w_s0(1),r_s0(1),'r^','MarkerSize',12);hold on;
        plot(th,Tn(1,1)*ones(size(th)),'b','LineWidth',1.2);
        legend('estimated spectral line','estimated noise variance','Location','North');
        arrow(w_s0,r_s0,12);
        set(gca,'xlim',[0 2*pi]);
        set(gca,'YTick',[]);