eofcent.m

多通道奇异谱分析程序

function [F,L,B]=eofcent(X,n);
%  EOF calculates the empirical orthogonal functions
%  and amplitudes (principal components) of the data matrix 'X'.
% Syntax: [F,L,B]=eof(X); [F,L,B]=eof(X,.9);
%
%  Input:        X - data matrix. For a standard (S-mode) EOF analysis,
%                    the columns of X are time series, while the rows
%                    are spatial maps. The eigenfunctions in this case
%                    will be spatial patterns, and the principal
%                    components are time series.
%                n - number of eigenfunctions to return (optional).
%                    If n is less than 1, it is interpreted as
%                    a fractional variance (e. g. n=.9), and enough
%                    eigenvectors are returned to account for n*100%
%                    of the variance. The default is to return all EOFs.
%     
%  Output:       F - eigenfunction matrix (columns are eigenvectors).
%                L - vector of eigenvalues.(all eigenvalues are returned)
%                B - principal components matrix.
%
% EOFCENT does the same thing as EOF, but does not allow the data matrix to
% be modified within the function, thus avoiding the memory penalty of passing
% the large data matrix into the function. If you want to center or 
% standardize the data, you must do it in the main workspace before calling
% EOFCENT The commands "[r,c]=size(X); X=X-ones(r,1)*mean(X);" will center the
% data. If you then want to standardize the data, use "X=X./(ones(r,1)*std(X));".
%
%  Written by Eric Breitenberger.    Version date 1/11/96
%  Please send comments and suggestions to [e-mail]   
%
[r,c]=size(X);
if c>r, disp('Warning: Covariance matrix may be ill-conditioned.'), end 
if nargin==1
  n=c;
end

S=X'*X;   % compute the covariance matrix       
[F,L]=eig(S);
clear S

%  sort eigenvectors, eigenvalues
[L,i]=sort(diag(-L));
L=-L';  
F=F(:,i);

% figure out how many eigenvectors to keep:
if n<1 % if n is in the form of fractional variance, convert to an index
  var=n*sum(L);
  i=find(cumsum(L)>=var);
  n=i(1);
end  

if c>n, F=F(:,1:n); end  % keep only first n eigenvectors
B=X*F;                   % calculate principal components (first n)