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);