Global Index (all files) (short | long) | Local Index (files in subdir) (short | long)
[evector,ew,amps,fac,covmat]=geteof(data,col,filtpar,area);
function [evector,ew,amps,fac,covmat]=geteof(data,[col],[filtpar],[area]);
calculates the eigen-vectors and -values of the covariance matrix of 'data'
the EOFs 'ev' are the solution of (data*data')/(length(dim2)-1) * ev = ev * ew
the 'amps' solve data = ( amps' * ev' )'
NaN's are possible but the amplitudes are not any longer orthogonal!
input : data - the data set (no data : NaN)
assumed to have 1 temporal and 1 or more
spatial dimensions
col [1] - column of the temporal dimension
filtpar [0] - -1 for no tranposing
-2 for no normalization
>0 for filtering covariance matrix
filtpar is the 3rd argument for filt2d
DON'T DO THIS !!!
the next is experimental !!!
area [1] - one weight for EACH grid point (in time
and space !)
output : evector - 'vertical' (temporal) eigen-vectors of 'data'
normalized to 1
ew - corresponding eigen-values
ew/sum(ew) gives explained variance
amps - (spatial) amplitudes of evector
sorted in columns (as 'evector')
fac - number of realisations of pairs
covmat - covariance matrix
uses : nans.m filt2d.m
version 1.3.1 last change 30.07.2002
| This function calls | |
|---|---|
function [evector,ew,amps,fac,covmat]=geteof(data,col,filtpar,area);
% Gerd Krahmann, IfM Kiel, Feb 1993
%
% Ref.: U.Send's EOF notes
% and
% http://www.dkrz.de/forschung/reports/report11-1.1/pingo-45.html
% for the area stuff
% added filtering mode G.Krahmann, Sep 1994
% changed filtering mode G.Krahmann, Jan 1995
% changed output order G.Krahmann, Sep 1995
% transponing for speed G.Krahmann, Feb 1996
% removed problem with amps
% and introduced 'gap_par' G.Krahmann, Mar 1996
% added compatibility to MATLAB 5 G.Krahmann, LODYC Paris, Aug 1997
% changed multi-dim handling G.Krahmann, LODYC Paris, Aug 1997
% added area handling
% and removed comp. to MATLAB 4 G.K., LDEO Jul 2000
% bugfix G.K., LDEO Jul 2002 1.3.0-->1.3.1
% add weights, if demanded
if nargin>3
disp('ask Gerd, if this fails')
data = data.*sqrt(area);
end
% check for MATLAB version and reshape/reorder if necessary
nd=ndims(data);
if nargin>1
if isempty(col)
col=1;
end
else
col=1;
end
dimorder=[1:nd];
dimorder(col)=nan;
dimorder=dimorder(find(~isnan(dimorder)));
dimorder=[col,dimorder];
if nd>2
data=permute(data,dimorder);
oldsize=size(data);
data=data(:,:);
end
% transpose the data matrix if s(1)>s(2)
% to get increased speed
sd=size(data);
if sd(1)>sd(2)
data=data';
transp=1;
if nargin>2
if filtpar==-1
data=data';
transp=0;
else
disp('Matrix transposed for increased speed!')
end
else
disp('Matrix transposed for increased speed!')
end
else
transp=0;
end
% check if NaN's exist in data
% if yes, calculate the number of data-pairs for calculating the
% covariance matrix separately for each pair, else use the
% length of the series as number of pairs
h1=~isnan(data);
if (length(data(:))~=sum(h1(:)))
nan_exist=1;
s=size(data);
fac=zeros(s(1));
for i=1:s(1)
for k=1:s(1)
fac(i,k)=length( find( h1(i,:).*h1(k,:) ) ) ;
end
end
allnan=find(sum(h1)==0);
fac = fac-1;
if ~isempty(find(fac<=0))
disp('ask Gerd')
end
else
nan_exist=0;
s=size(data);
fac=s(2);
allnan=[];
end
% replace NaNs with 0, that gives no contribution when calculating the
% covariance matrix
[data,nan_ind]=nans(data,0,nan);
% calculate covariance matrix
% and divide by number of realizations
covmat=data*data';
covmat=covmat./fac;
h=find(isnan(covmat));
if length(h) > 0
covmat(h)=zeros(1,length(h));
end
% filter covariance matrix, if wanted
if (nargin>2)
if filtpar<1 & filtpar~=-1 & filtpar~=-2 & filtpar>0
covmat=filt2d(covmat,2,filtpar);
end
end
% calculate Singular Value Decomposition
[u,s,v]=svd(covmat);
evector=u;
ew=s;
ew=diag(ew);
% normalize Eigenvectors to rms=1
if transp==0
norms=0;
if exist('filtpar')
if filtpar==-2
norms=1;
end
end
if norms==0
rev=rms(evector);
evector=evector./(ones(size(evector,1),1)*rev);
end
end
% project to get amplitudes
amps=data'/evector';
% if the data matrix was transposed, reorder output arguments
% and correct for rms(evector)=1
% at the same time check for time series full of NaN's
if transp==1
ew=ew/sd(2)*sd(1);
evn=amps(:,1:length(ew));
ampsn=evector(:,1:length(ew));
evector=evn;
amps=ampsn;
if ~isempty(allnan)
evector(allnan,:)=evector(allnan,:)*nan;
end
rev=rms(evector);
evector=evector./(ones(size(evector,1),1)*rev);
amps=amps.*(ones(size(amps,1),1)*rev);
else
if ~isempty(allnan)
amps(allnan,:)=amps(allnan,:)*nan;
end
ew=ew(1:sd(1));
end
% check for MATLAB version and reshape if necessary
if exist('oldsize')
if length(oldsize)>2
if transp==0
amps=reshape(amps',[size(amps,2),oldsize(2:length(oldsize))]);
else
amps=reshape(amps',[size(amps,1),oldsize(2:length(oldsize))]);
end
[dummy,ind]=sort(dimorder);
amps=permute(amps,ind);
end
end
% add weights, if demanded
if nargin>3
disp('ask Gerd, if this fails')
if ndims(amps)==2
amps = amps./sqrt(area');
else
amps = amps./sqrt(area);
end
end