function [V,EOFs,EC,error]=EOF(D,p)

% function [V,EOFs,EC,error]=EOF(D,p)
%
% This function decomposes a data set D into its EOFs分解一个数据集 D为EOFS 
% and eigenvalues V.特征值V The most effecient algorithm proposed by
% von Storch and Hannostock (1984) is used.
%
% D = each row is assumed to be a sample每一行假定为一个样本. Each column a
% variable.每一列为一变量
% Thus a column represents a time series of one variable因此一列代表一个时间序列的变量
% p = is an optional  indicating the number of 可选项 指示EOFS的数量
% EOFs the user wants to take into account (smaller than
% the number of time steps and the number of samples).
%
% EOFs= matrix with EOFs in the columns
% V = vector with eigenvalues associated with the EOFs与EOFS有关的特征值向量
% EC = EOF Coefficients系数
% error = the reconstruction error (L2-norm)重建的误差
% 
% This function uses svds from Matlab version 5
%
% Written by Martijn Hooimeijer (1998)


% first compute zero-averaged data sets
%DS=zeroavg(D);

% Determine size of the data matrix (m=number of samples, n= number of variables)
[m,n]=size(D);
q=min(m,n);

% Determine singular value decomposition of the problem
if nargin < 2
   [U,S,F]=svds(D,min(m,n));
else
   [U,S,F]=svds(D,min(q,p));
end

% rewrite the large eigenvalue matrix to a vector and 
% apply the appropriate normalisation把大的矩阵特征值写成一个变量并标准化
V=diag(S).^2/(m-1);

% Define the EOFs (with EOFs in the columns)
EOFs=F;

% Determine the EOF coefficients (form the zero-averaged data)
EC=U*S;

% Determine the difference between the original data and the 
% reconstructed data确定原始数据和重建数据的不同
diff=(D-EC*EOFs');

% determine the L2 error norm for each variable
error=sqrt(sum(abs(diff.^2)));