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
9
Y(:,m) = X((1:N-M+1)+m-1);
12
%% Choose covariance estimation
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
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
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
33
buf = PC(:,m)*RHO(:,m)'; % invert projection
34
buf = buf(end:-1:1,:);
35
% Anti-diagonal averaging
37
RC(n,m) = mean(diag(buf,-(N-M+1)+n));