eeof.zip_EEOF

function [E,V,A,C]=eeof(X, M, convert) % Syntax: [E,V,A,C]=eeof(X, M); [E,V,A,C]=eeof(X, M, 1); % This function performs an extended empirical orthogonal % function (EEOF) analysis of matrix 'X', for embedding dimension 'M'. % Each of the L columns of X is a time series of length N. % % Returns: E - eigenfunction matrix. (LM by LM) % V - vector containing variances (unnormalized eigenvalues). % A - matrix of principal components. % C - lag-covariance matrix. % % V is ordered from large to small: E and A are sorted accordingly. % % Note that X is assumed to be centered. To center the data, use % the commands: % [r,c]=size(X); X=X-ones(r,1)*mean(X); before running EEOF. % If you also want to standardize the data, use: % X=X./(ones(r,1)*std(X));. % % If a third argument is supplied, the eigenfunctions/values will % be reordered into the same format as MSSA output - i. e. L blocks % of size M rather than M blocks of size L. % % This function provides the same output, within numerically determined % limits, as MSSA methods using Broomhead-King type covariance estimation: % it is intended as a check on those functions. % % Note that this function is *extremely* computationally intensive % for large matrices and lags. For example, if X is 1000 by 1000, % and M = 5, EEOF will take about 10 hours on a Cray YMP! Inputting % a subset of the PCs of X rather than the full data matrix can % substantially reduce the computational load. % % Written by Eric Breitenberger. Version date 1/11/96 % Please send comments and suggestions to eric@gi.alaska.edu % [N,L]=size(X); if M*L>=N-M+1, disp('Warning: Covariance matrix may be ill-conditioned.'), end % Create the extended matrix: T=zeros(N-M+1,M*L); for i=1:M T(:,L*(i-1)+1:L*i)=X(i:N-M+i,:); end % Compute the eigenvectors/values of the covariance matrix: C=(T'*T)/(N-M+1); clear X [E,V]=eig(C); V=diag(V); A=T*E; % compute principal components if nargin==3 % Prepare MSSA-style output: % sort E,V,C, and A from M blocks of L to L blocks of M. ind=1:L:(M-1)*L+1; for i=1:L, index=[index ind+i-1]; end E=E(index,index); V=V(index); % sort the covariance matrix and PCs: C=C(index,index); A=A(:,index); end % Sort eigenvalues/vectors/PCs in descending order: [V,ind]=sort(-V); V=-V'; E=E(:,ind); A=A(:,ind);