ssa.sci ( https://gist.github.com/werediver/9785544 ) adapted for GNU Octave

ssa-octave.m

% 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