mrc.m

多通道奇异谱分析程序

  function [R]=mrc(A,E,L,index)
% Syntax: [R]=mrc(A,E,L); [R]=mrc(A,E,L,index);
% This function calculates the 'reconstructed components' using the 
% eigenvectors (E, from MSSAEIG) and principal components (A, from MPC).
% Both of these matrices may have had columns trimmed but they must
% have the same number of columns. 
%
% Input:  A - principal components matrix (N-M+1 by k)
%         E - eigenvector matrix (L*M by k)
%         L - number of channels in the original data matrix.
%     index - (optional) index vector into A and E. For example, if 
%             index = [1 4], the RCs for EOF/PC pairs 1 and 4 will
%             be computed. The default is to calculate an RC for 
%             each input EOF/PC pair.
%
% Output: R - Matrix of RCs, N by L*k. Each PC/EOF pair has L RCs - one
%             for each channel. These are arranged in R as follows: the
%             i-th section (of L columns) contains the L RCs for the
%             i-th PC/EOF pair. 
%
%  See section 2.4 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]   

[ml,k]=size(E);
[ra, ka]=size(A);
if k~=ka, error('E and A must have the same number of columns.'), end
if nargin==3, index=1:k; end

M=ml/L;   %     These lines assume that E and A
N=ra+M-1; %     have the "right" row dimensions.

R=zeros(N,L*length(index));
Z=zeros(M-1,k);
A=[A' Z']; % zero pad A to make it N by k
A=A';
Ej=zeros(M,L);

% Calculate RCs
for j=1:length(index)
  Ej(:)=E(:,index(j)); % Convert the j-th EOF into a matrix of filters
  for i=1:L     % Compute the RCs for the j-th EOF/PC pair
    R(:,(j-1)*L+i)=filter(Ej(:,i),M,A(:,index(j)));
  end
end

% Adjust first M-1 rows and last M-1 rows
for i=1:M-1
  R(i,:)=R(i,:)*(M/i);
  R(N-i+1,:)=R(N-i+1,:)*(M/i);
end