mssaeig.m

多通道奇异谱分析程序

  function [E,V,C]=mssaeig(X, M,method,calc)
%  Syntax: [E,V]=mssaeig(X, M, 'large');  [E,V,C]=mssaeig(X, M, 'BK');  
% MSSAEIG starts an MSSA of matrix 'X', for embedding dimension 'M'.
% It builds the covariance matrix and calculates the eigenvectors/values.
%
% Input:  X - Data matrix, organized such that each of the L columns of
%             X is a time series of length N.
%
%         M - Embedding dimension M.
%
%    method - (Optional) The method of calculating the covariance matrix.
%             Can be 'unbiased' - N-k weighting (default)
%                      'biased' - N weighting
%                          'BK' - (Broomhead/King) based on trajectory matrix
%
%      calc - (Optional) Can be either 'small' or 'large' (default).
%             If calc='large', the FFT-based COVAR routine will be used, which
%             is faster for problems with long series, many channels, or large
%             lags. If calc='small', the direct summations in AC and XC will
%             be used, which are often faster for problems with short series,
%             few channels, or small lags. If BK is chosen for the covariance
%             method, this argument does not apply.
%
% Output: E - eigenfunction matrix in standard form, where the
%             i-th column contains, in order, the L i-th order 
%             eigenvectors (T-EOFs). 
%         V - vector containing variances (unnormalized eigenvalues).
%         C - lag-covariance matrix.
%
% See section 2 of Vautard, Yiou, and Ghil, 1992, Physica D 58, 95-126;
% and Plaut and Vautard, 1994, J. Atm. Sciences 51, 210-236.
%
% Written by Eric Breitenberger.    Version date 1/11/96
% Please send comments and suggestions to [e-mail]   
%

if nargin==2
  method='unbiased';
  calc='large';
elseif nargin==3
  arg=findstr(method, 'BK unbiased biased large small');
%                      ^  ^        ^      ^     ^
%                      1  4        13     20    26
  if arg==1
    method='BK';
    calc='none';
  elseif arg==4
    method='unbiased';
    calc='large';
  elseif arg==[6 13]
    method='biased';
    calc='large';
  elseif arg==20
    method='unbiased';
    calc='large';
  elseif arg==26
    method='unbiased';
    calc='small';
  elseif calc~=1
    error('Improper input arguments.')
  end
elseif nargin==4
  arg=findstr(calc, 'large small');
  %                  ^     ^
  %                  1     7
  if arg==1
    calc='large';
  elseif arg==7
    calc='small';
  else 
    error('Improper input arguments.')
  end
elseif nargin~=1
  error('Incorrect number of input arguments.')
end

[N,L]=size(X);
if M>=N-M+1, disp('Lag M is too large - covariance matrix may be ill-conditioned.'), end
if strcmp(method,'BK') & ~strcmp(calc,'none'), error('Cannot specify ''large'' or ''small'' with ''BK''.'), end

if ~strcmp(calc,'none')
  disp(['Using ' method ' covariance estimation with ' calc ' matrix computational method.'])
else
  disp(['Using ' method ' covariance estimation.'])
end

% Calculate the covariance matrix according to the options:
if strcmp(method,'BK')
  C=bk(X,M);

elseif strcmp(calc,'large')  % Use COVAR
  c=covar(X,M-1,method);
  % Build the covariance matrix T, block by block:
  C=zeros(L*M,L*M);

  % First, build the diagonal blocks:
  for i=1:L
    q=c(M:2*M-1, (i-1)*(L+1)+1);
    Tij=toeplitz(q);    % create Toeplitz block
    C((i-1)*M+1:i*M, (i-1)*M+1:i*M)=Tij;
  end

  % Now build all the other blocks:
  for i=1:L
    for j=i+1:L
      q=c(:,(i-1)*L+j);
      Tij=toeplitz(q(M:-1:1),q(M:2*M-1)); 
      C((i-1)*M+1:i*M, (j-1)*M+1:j*M)=Tij;  % Above diagonal
      C((j-1)*M+1:j*M, (i-1)*M+1:i*M)=Tij'; % Below diagonal
    end
  end

elseif strcmp(calc,'small')   % Use XC and AC

  % Build the covariance matrix T, block by block:
  C=zeros(L*M,L*M);

  % First, build the diagonal blocks:
  for i=1:L
    c=ac(X(:,i), M-1, method);  % calculate autocovariance estimates
    Tij=toeplitz(c);    % create Toeplitz block
    C((i-1)*M+1:i*M, (i-1)*M+1:i*M)=Tij;
  end

  % Now build all the other blocks:
  for i=1:L
    for j=i+1:L   % Above diagonal
      q=xc(X(:,i),X(:,j),M-1,method); % calculate cross-covariance estimates
      Tij=toeplitz(q(M:-1:1),q(M:2*M-1)); 
      C((i-1)*M+1:i*M, (j-1)*M+1:j*M)=Tij;  % Above diagonal
      C((j-1)*M+1:j*M, (i-1)*M+1:i*M)=Tij'; % Below diagonal
    end
  end
end

% Finally, calculate and sort the eigenvectors/values of C
[E,V]=eig(C);
% Sort eigenvalues/vectors
[V,i]=sort(diag(-V));
V=-V';  
E=E(:,i);