eof.m

多通道奇异谱分析程序

function [F,L,B]=eof(X,n,s);

%  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,'norm');
%
%  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.
%                s - Normalization option. If s='norm', then each
%                    column of X will be normalized (assigned
%                    unit variance). If s is not specified, the 
%                    data are not normalized. 
%     
%  Output:       F - eigenfunction matrix (columns are eigenvectors).
%                L - vector of eigenvalues.(all eigenvalues are returned)
%                B - principal components matrix.
%
%  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; s='none';
elseif nargin==2
  if isstr(n)
    s=n; n=c;
  else
    s='none';
  end
end

X=X-ones(r,1)*mean(X); % center the data
if s=='norm'
  X=X./(ones(r,1)*std(X)); % normalize
elseif s~='none'
  error('Improper normalization option. Please check inputs.') 
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)