function [eig_val,eig_vec,PC,RC,RCp]=SSA(X,varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function SSA : performs Singular Spectrum Analysis on time series X
%       with the method of Vautard and Ghil, Phys. D. 1989
%
% syntax : [eig_val,eig_vec,PC,RC,RCp]=SSA(X,[M, K])
%   Inputs 
%   X :  line vector
%   M : window length. Default value is M=length(X)/10.
%   K : number of EOFs used for reconstruction (0 by default)
%
% output : 
%   eig_val : eigenvalue spectrum
%   eig_vec : eigenvector matrix ("temporal EOFs")
%   PC      : matrix of principal components
%   RC      : matrix of RCs (N*M, K) (only if K>0)
%   RCp     : Reconstructed timeseries
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% first make sure the vector size is 1*N
if size(X,1)~= 1
    X=X';
end
N=length(X); 
Xr=(X-mean(X)); % center the series
% set default value for M
if nargin('SSA')==0
    M=round(N/10);
    K=0;
elseif nargin('SSA')==-1
        M=varargin{1};
        K=0;
elseif nargin('SSA')==-2
        M=varargin{1};   
        K=varargin{2};
end
Np=N-M+1;
% compute autocorelation
[gam,lags]=xcorr(Xr,M-1,'unbiased');
% Fill in Covariance matrix
C=toeplitz(gam(M:2*M-1)); %take postive half of autocorrelation diagram, hence M to 2M-1
%C=toeplitz(gam(1:M)); 

% apply SVD

% solve eigenvalue problem
[V,D]=eig(C);
eig_vec=fliplr(flipud(V)); % flip them to have them ranked in decreasing order
eig_val=fliplr(flipud(D));

% compute PCs
decal=zeros(Np,M);
for t=1:N-M+1
    decal(t,:)=Xr(t:M+t-1);
end
PC=decal*eig_vec; % the columns of this matrix are Ak(t), k=1 to M


% compute reconstructed timeseries if K > 0
if (K>=1)
    RC=zeros(N,K);
     % first M terms   
     for t=1:M-1
        Av=flipud(PC(1:t,1:K));
        eig_vec_red=eig_vec(1:t,1:K);
        RC(t,:)=1/t*sum(Av.*eig_vec_red,1);    
     end
     %middle of timeseries
     for t=M:Np
        Av=flipud(PC(t-M+1:t,1:K));
        eig_vec_red=eig_vec(1:M,1:K);
        RC(t,:)=1/M*sum(Av.*eig_vec_red,1);    
     end  
     % Last M terms
     for t=Np+1:N
        Av=flipud(PC(t-M+1:Np,1:K));
        eig_vec_red=eig_vec(t-N+M:M,1:K);
        RC(t,:)=1/(N-t+1)*sum(Av.*eig_vec_red,1);  
     end  
     %sum and restore the mean
     RCp=sum(RC,2)+mean(Xr);
 else 
     RC=[];
     RCp=[];
end