Skip to content

Instantly share code, notes, and snippets.

@werediver werediver/ssa-octave.m
Last active May 5, 2016

Embed
What would you like to do?
ssa.sci ( https://gist.github.com/werediver/9785544 ) adapted for GNU Octave
% Force GNU Octave to interpret this file as a script
1;
function [y_out, m, sig] = ssa_normalize(y)
% http://en.wikipedia.org/wiki/Student%27s_t-statistic
m = mean(y);
sig = sqrt(variance(y)) % stdev(y);
y_out = (y - m) / sig;
endfunction
function [y_out] = ssa(y, L, I)
[LAMBDA, U, V] = ssa_decompose(y, L);
y_out = ssa_reconstruct(LAMBDA, U, V, I);
endfunction
function [LAMBDA, U, V, X] = ssa_decompose(y, L)
% y expected to be a column vector
% Stage 1: Decomposition
% Step 1: Embedding
N = length(y);
if (L > N / 2)
L = N - L;
end
K = N - L + 1;
% Time-delayed embedding of y, the trajectory matrix
X = zeros(L, K);
for i = 1 : K
X(:, i) = y(i : i + L - 1);
end
% Step 2: Singular value decomposition
% Direct approach
% _V will never be used
%[U, LAMBDA, _V] = svd(X)
%LAMBDA = diag(LAMBDA)
% Covariance matrix
C = X * X' / K;
% For stationary series
%C = toeplitz(corr(X, L));
% _V will never be used
[U, LAMBDA, _V] = svd(C);
% The eigenvalues of C are the squared eigenvalues of X
LAMBDA = sqrt(diag(LAMBDA));
% Principal components
V = X' * U;
for i = 1 : L
V(:, i) = V(:, i) / LAMBDA(i);
end
endfunction
function [y] = ssa_reconstruct(LAMBDA, U, V, I)
% I expected to be a row vector
[K, L] = size(V);
N = K + L - 1;
% Stage 2: Reconstruction
% Step 3: Grouping
for i = I
LAMBDA_U(:, i) = LAMBDA(i) * U(:, i);
end
% Reconstructed components
%X = LAMBDA_U(:, I) * V_transposed(I, :);
X = LAMBDA_U(:, I) * V(:, I)';
% Step 4: Diagonal averaging
y = zeros(N, 1);
for i = 1 : K + L - 1
v = adiag(X, i);
y(i) = sum(v) / length(v);
end
y = real(y);
endfunction
function [v] = adiag(x, z)
% 2D matrix antidiagonal
[N, M] = size(x);
% The total diagonals count
Z = N + M - 1;
% The current diagonal length
z_len = min([z, Z - z + 1, N, M]);
i = max([z, N * (z - N + 1)]);
step = N - 1;
v = x(i : step : i + step * z_len - 1);
endfunction
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.