Global Index (all files) (short | long) | Local Index (files in subdir) (short | long)
[ev,ew,amps,rec]=getmssa(data,llag,ldim);
function [ev,ew,amps,rec]=getmssa(data,llag,[ldim]);
calculate (multichannel) singular spectrum analysis MSSA
for SSA simply give a vector as data-array
the method is similar to extended EOFs (EEOF)
NaN's are possible and internally replaced by random values
they should not exceed 10% of the series
input : data - the data-array (max 2-dim)
llag - length of lag (should not exceed one third
of the length of the series)
[ldim] [1] - dimension to be lagged
output : ev - 'vertical' eigen-vectors of lagged 'data'
normalized to 1
MATLAB 5: dim1=lag dim2=series dim3=EOF#
MATLAB 4: dim1=lag+series dim3=EOF#
ew - corresponding eigen-values
ew/sum(ew) gives explained variance
amps - amplitudes of ev
sorted in columns (like 'ev')
rec - reconstructed fields (ev*amps)
only available for MATLAB 5 or higher
Warning: this is very slow !
version 0.1.1 last change 15.09.1997
| This function calls | |
|---|---|
function [ev,ew,amps,rec]=getmssa(data,llag,ldim);
% Gerd Krahmann, LODYC Paris, Sep 1997
%
% Ref.: Y.Sezginer Unal & M.Ghil, Clim. Dyn., 1995, 11:255-278
% added 'rec' G.Krahmann, LODYC/Paris Sep 1997, 0.1.0-->0.1.1
% check if NaN's exist in data
% if yes, replace them by white-noise random values
if isnan(sum(data(:)))
dummy = sum(data);
dummy = find(isnan(dummy));
disp('WARNING: filling gaps with white noise')
for n=1:length(dummy)
bad = find(isnan(data(:,dummy(n))));
va = nstd(data(:,dummy(n)));
me = nmean(data(:,dummy(n)));
if length(bad)>length(data)/10
disp('WARNING: more than 10% gaps in series')
end
data(bad,dummy(n)) = randn(length(bad),1)*va + me;
end
end
% check for length of lag in relation to length of time series
if llag>size(data,1)/3
disp('WARNING: lag is longer than 1/3 of time series')
end
% remove mean of time series
data = data - repmat(mean(data),[size(data,1),1]);
% reorder data to lagged series
ind = zeros(llag,size(data,1)-llag+1,size(data,2));
dummy = [1:llag]';
ind = ind +...
repmat(dummy,[1,size(data,1)-llag+1,size(data,2)]);
dummy = [0:size(data,1)-llag];
ind = ind +...
repmat(dummy,[llag,1,size(data,2)]);
clear dummy
dummy(1,1,[1:size(data,2)]) = [0:size(data,2)-1]*size(data,1);
ind = ind +...
repmat(dummy,[llag,size(data,1)-llag+1,1]);
ind = permute(ind,[2,1,3]);
ind = ind(:,:)';
ndata = data(ind);
sn = size(ndata);
% transpose the data matrix if sn(1)>sn(2)
% to get increased speed
if sn(1)>sn(2)
ndata=ndata';
transp=1;
else
transp=0;
end
% calculate covariance matrix
% and divide by number of realizations
covmat=ndata*ndata'/size(ndata,2);
% calculate Singular Value Decomposition
[u,s,v]=svd(covmat);
ev=u;
ew=s;
ew=diag(ew);
% normalize Eigenvectors to rms=1
rev=rms(ev);
ev=ev./(ones(size(ev,1),1)*rev);
% project to get amplitudes
amps=ndata'/ev';
% if the data matrix was transposed, reorder output arguments
% and correct for rms(ev)=1
% at the same time check for time series full of NaN's
if transp==1
ew=ew/sn(2)*sn(1);
evn=amps(:,1:length(ew));
ampsn=ev(:,1:length(ew));
ev=evn;
amps=ampsn;
rev=rms(ev);
ev=ev./(ones(size(ev,1),1)*rev);
amps=amps.*(ones(size(amps,1),1)*rev);
end
% reshape results
if intvers>4
ev = reshape(ev,llag,size(ev,1)/llag,size(ev,2));
end
% reconstructed fields
if intvers>4 & nargout>3
rec=zeros(size(data,1),size(data,2),size(ev,3));
count=zeros(size(data,1),size(data,2),size(ev,3));
for n=1:size(ev,3)
for m=1:size(ev,1)
for k=1:size(ev,3)
rec(n+[0:size(ev,1)-1],:,k)=rec(n+[0:size(ev,1)-1],:,k)+...
ev(:,:,n)*amps(n,k);
count(n+[0:size(ev,1)-1],:,k)=count(n+[0:size(ev,1)-1],:,k)+1;
end
end
end
rec=rec./count;
end