ssa

1
function [C,LBD,RC] = SSA(N,M,X,nET)
2
%% Calculate covariance matrix C (Toeplitz approach)
3
%    covX  = xcorr(X,M-1,'unbiased');
4
%    Ctoep = toeplitz(covX(M:end));
5
%% Calculate covariance matrix (trajectory approach)
6
%  it ensures a positive semi-definite covariance matrix
7
   Y = zeros(N-M+1,M);
8
   for m=1:M
9
     Y(:,m) = X((1:N-M+1)+m-1);
10
   end
11
   Cemb = Y'*Y/(N-M+1);
12
%% Choose covariance estimation
13
%    C = Ctoep;
14
   C = Cemb;
15
%% Calculate eigenvalues LAMBDA and eigenvectors RHO
16
% Function eig returns two matrices,
17
% the matrix RHO with eigenvectors arranged in columns,
18
% the matrix LAMBDA with eigenvalues along the diagonal
19
   [RHO,LBD] = eig(C);
20
   LBD       = diag(LBD);           % extract the diagonal elements
21
   [LBD,ind] = sort(LBD,'descend'); % sort eigenvalues
22
   RHO       = RHO(:,ind);          % and eigenvectors
23
%% Calculate principal components PC
24
% The principal components are given as the scalar product
25
% between Y, the time-delayed embedding of X, and the eigenvectors RHO
26
   PC = Y*RHO;
27
%% Calculate reconstructed components RC
28
% In order to determine the reconstructed components RC,
29
% we have to invert the projecting PC = Y*RHO;i.e. RC = Y*RHO*RHO'=PC*RHO'
30
% Averaging along anti-diagonals gives the RCs for the original input X
31
   RC = zeros(N,nET);
32
   for m=1:nET
33
     buf = PC(:,m)*RHO(:,m)'; % invert projection
34
     buf = buf(end:-1:1,:);
35
% Anti-diagonal averaging
36
     for n=1:N
37
       RC(n,m) = mean(diag(buf,-(N-M+1)+n));
38
     end
39
end
40