%singular spectral analysis
% program is written 17.02.2009 by L.V. Zotov
% do not forget about hank.m file


clear;
addpath('/home/leonid/lectures/labs/');% path to hank.m file
N_signal=1024;
% generating two-sin signal
for (k=1:1:N_signal)
      garm(k)=k/50+5*sin(2*pi/20*(k-1))+k/100*sin(2*pi/200*(k-1));
end;

plot(garm);


% ARMA process generating
ar(1)=0.5*randn(1);
ar(2)=-0.2*ar(1)+0.5*randn(1);

for (i=3:1:N_signal)
    ar(i)=0.9*ar(i-1)-0.7*ar(i-2)+0.5*randn(1);
end;

plot(ar);

%making a signal
signal=garm+5*ar;

plot(signal);

%--------------------------------------------------------------------------
cd '/home/leonid/lectures/labs/lab6';

fout=fopen('signal.dat','wt');
for(i=1:1:N_signal)
 fprintf(fout,'%12.6f%12.6f\n',i,signal(i)); 


end;
fclose(fout);

% parameter - caterpillar length
L=220;

%matrix
K=N_signal-L;
A=zeros(L,K);
  for ii=1:1:L    
   for j=1:1:K
    A(ii,j)=signal(ii+j-1);
   end;
  end;
  
 %decomposition
[U,S,V] = svd(A);
K_max=max(L,K)
L_min=min(L,K)


%10 components evaluation
for(i=1:1:10)
 clear SX;
 clear G;
 SX=zeros(L,N_signal-L);
 SX(i,i)=(S(i,i));
 G=U*SX*V';
 RX(i,:)=hank(G,L,N_signal,L_min,K_max);
end;

%weights calculation
for(i=1:1:N_signal)
    if(i<=L_min)
        w(i)=i;
    elseif(i<=K_max)
        w(i)=L_min;
    else
        w(i)=N_signal-i;
    end;
end;

WK=zeros(10,10);
for(i=1:1:10)
   for(j=1:1:10)
     for(k=1:1:N_signal)
      WK(i,j)=WK(i,j)+RX(i,k)*w(k)*RX(j,k);  
     end; 
   end;
end;

%WK(1,1)=0;
[I,J] = meshgrid([1:1:10]);
pcolor(I,J,WK)

%plot components
subplot(10,1,1); plot(RX(1,:)); title('1')
subplot(10,1,2); plot(RX(2,:)); title('2')
subplot(10,1,3); plot(RX(3,:)); title('3')
subplot(10,1,4); plot(RX(4,:)); title('4')
subplot(10,1,5); plot(RX(5,:)); title('5')
subplot(10,1,6); plot(RX(6,:)); title('6')
subplot(10,1,7); plot(RX(7,:)); title('7')
subplot(10,1,8); plot(RX(8,:)); title('8')
subplot(10,1,9); plot(RX(9,:)); title('9')
subplot(10,1,10); plot(RX(10,:)); title('10')

%grouping
 plot(signal,'black');
 hold on;
 plot(RX(1,:),'red')
 plot(RX(2,:)+RX(3,:))
 plot(RX(4,:)+RX(5,:),'green')
 hold off;

 all_ssa=RX(1,:);
 for(i=2:1:10)
     all_ssa=all_ssa+RX(i,:);
 end;
 
 ost=signal-all_ssa;
 

  hold on
  plot(ost,'black')
  plot(RX(6,:)+RX(7,:)+RX(8,:)+RX(9,:)+RX(10,:),'cyan')
  plot(ar,'magenta')