covar.m

多通道奇异谱分析程序

function c=covar(X,Y,maxlag,method)

% COVAR: Covariance function estimates.
% Syntax: c=covar(X,Y,lag,method), covar(X)
%
% COVAR calculates auto- or cross-covariance estimates for time series. It
% can handle matrices or vectors. A maximum lag and various estimation
% methods can also be specified.
% 
% Input: The first two arguments X and Y are handled as follows:
%
% 1) If a single vector X is input, COVAR calculates the auto-covariance
%    function for series X. The result is output in c, which has 2k+1 elements.
%    The k+1th element is the covariance at lag zero; succeeding elements
%    k+1:2k+1 are the covariances at lags 1 to k. The vector is symmetric
%    about the midpoint k+1. Seems silly, but it's this way to be consistent
%    with the matrix input/output.
%
% 2) If two vectors X and Y (of the same length) are input, COVAR calculates
%    the cross-covariance between the series. The result is output in c, 
%    which has 2k+1 elements. c(k+1) corresponds to zero lag. The first k
%    elements contain the positive lags, the last k contain the negative lags.
%    By convention, X LEADS Y FOR POSITIVE LAGS.
%
% 3) If a single matrix X is input, where the columns are time series, the
%    lagged covariances for all column pairs will be calculated. If X has 
%    L columns, the output c will be a matrix of size 2k+1 by L^2. The first
%    L columns of c will contain the covariances between column 1 of X and
%    the other columns, the next L columns will contain the covariances for
%    column 2 of X, and so on. Each of these columns is of length 2k+1,
%    arranged as described above for cross-covariances.
%
% Optionally, the maximum lag may be specified in the next argument, after
% X and/or Y. By default, all possible lags are computed.
%
% In the final input argument, the method of covariance estimation 
% may be specified as one of: 
%                 'biased'   - (N-weighted or Yule-Walker)  
%                 'unbiased' - (N-k weighted - this is the default)
%                 'VG'       - (Vautard/Ghil, same as 'unbiased')
%
%  Written by Eric Breitenberger.   Version date 1/11/96
%  Please send comments and suggestions to [e-mail]   
% 

% Parse the input arguments:
% Flags in this section:
Xv=0; % X is a vector 
Xm=0; % X is a matrix 
Yv=0; % Y is a vector 

[N,L]=size(X);
if min(N,L)==1   % X is a vector
  X=X(:); 
  [N,L]=size(X); 
  Xv=1; 
else
  Xm=1;          % X is a matrix
end


if nargin==1
  method='unbiased';
  maxlag=N-1;
elseif nargin==2
[Ny,Ly]=size(Y);
  if isstr(Y)          % covar(X,method) has been called
    method=Y;
    maxlag=N-1;
  else
    if max(Ny,Ly)==1   % covar(X,maxlag) has been called
      maxlag=Y;
      method='unbiased';
    else               % covar(X,Y) has been called
      if min(Ny,Ly)~=1, error('X and Y cannot both be matrices.'), end
      Yv=1;
      Y=Y(:);
      if length(Y)~=N, error('X and Y must be the same length.'), end
      maxlag=N-1;
      method='unbiased';
    end
  end
elseif nargin==3
  [Ny,Ly]=size(Y);
  if min(Ny,Ly)~=1, error('X and Y cannot both be matrices.'), end
  if max(Ny, Ly)==1  % covar(X,lag,method) has been called
    method=maxlag;
    maxlag=Y;
  else
    Yv=1;
    Y=Y(:);
    if length(Y)~=N, error('X and Y must be the same length.'), end
    if isstr(maxlag)   % covar(X,Y,method) has been called
      method=maxlag;  
      maxlag=N-1;
    else               % covar(X,Y,maxlag) has been called
      method='unbiased';
    end  
  end
elseif nargin==4
  [Ny,Ly]=size(Y);
  if min(Ny,Ly)~=1, error('X and Y cannot both be matrices.'), end
  Yv=1;
  Y=Y(:);
  if length(Y)~=N, error('X and Y must be the same length.'), end
else
  error('Too many input arguments!')
end

if maxlag>=N, error('Too big a lag!'), end
if strcmp(method, 'VG'), method='unbiased'; end 

nfft=2.^(fix(log(N-1)/log(2))+2); % twice the next power of two


% Use FFT to compute the covariances:
if (Xv & Yv)       % Cross-covariance of X and Y.
  x=fft(X,nfft);
  y=fft(Y,nfft);
  c=ifft(y.*conj(x));
  c=real(fftshift(c));  % Zero lag is at 1+nfft/2
  c=c(1+nfft/2-maxlag:1+nfft/2+maxlag);
elseif Xv          % Autocovariance of X
  x=fft(X,nfft);
  c=ifft(x.*conj(x));
  c=real(fftshift(c));
  c=c(1+nfft/2-maxlag:1+nfft/2+maxlag);
elseif Xm          % Form covariance matrix of X
  c=zeros(2*maxlag+1,L^2);
  xtmp=zeros(nfft,L);
  for i=1:L
    xtmp(:,i)=fft(X(:,i),nfft);
  end
  ind=1+nfft/2-maxlag:1+nfft/2+maxlag;
  for i=1:L
    for j=i:L
      y=ifft(xtmp(:,j).*conj(xtmp(:,i))); 
      y=real(fftshift(y));
      col=(i-1)*L+j;
      c(:,col)=y(ind);
    end
  end
  ind=2*maxlag+1:-1:1;
  for i=2:L
    for j=1:i-1
      col=(i-1)*L+j;
      colnew=(j-1)*L+i;
      c(:,col)=c(ind,colnew);
    end
  end
end

if strcmp(method, 'VG'), method='unbiased'; end

if strcmp(method, 'unbiased')
  nrm=[N-maxlag:N  N-1:-1:N-maxlag]';
elseif strcmp(method,'biased')
  nrm=N*ones(2*maxlag+1,1);
end

if Xm
  c=c./(nrm*ones(1,L^2));
elseif Xv | Yv
  c=c./nrm;
end