clear;close all;clc

tic

% -------------------------------------------------------------------------

% -------------------------------------------------------------------------

% 求距平

path_0 = '*****';

path = [path_0,'sstsst_30s_30n_120e_80w.nc'];

file = netcdf(path,'nowrite');

sst = file{'sst'}(:);

[a,b,c] = size(sst);

sst1 = zeros(b*c,a);

for i = 1:a

   sst1(:,i) = reshape(sst(i,:,:),b*c,1);

end

% sst1列向量表示时间序列，如sst1(:,1)表示第一个月

% 1 -----------------------------------------------------------------------

sst2 = zeros(b*c,12);

for i = 1:a/12

   for j = 1:12

       k = (i-1)*12+j;

       sst2(:,j) = sst2(:,j)+sst1(:,k);

   end

end

sst_month_mean = sst2/(a/12);  % month mean

% 距平值

X = sst1;

for i = 1:a/12

   for j = 1:12

       k = (i-1)*12+j;

       X(:,k) = X(:,k)-sst_month_mean(:,j);

   end

end

X1 = X;

clear X

% 2 -----------------------------------------------------------------------

X = zeros(size(sst1));                            

for i = 1:12

   X(:,i:12:end) = sst1(:,i:12:end)-repmat(mean(sst1(:,i:12:end),2),1,size(sst1(:,i:12:end),2));

end

X2 = X;

% NOTE:

% (1) 两个方法都一样，都能够求出距平。很明显，方法2比较简单。

% -------------------------------------------------------------------------

% -------------------------------------------------------------------------

% eof计算

load('******X');

[m,n] = size(X);

% 1 -----------------------------------------------------------------------

% EOF

[EOF,E] = eig(X*X'/n);

PC = EOF'*X;

% reverse the order

E = fliplr(flipud(E));

% lambda = diag(E); % retain eigenvalues only

EOF_0 = fliplr(EOF);

PC_0 = flipud(PC);

% 2 -----------------------------------------------------------------------

% EOF时空转换

[EOF1,E1] = eig(X'*X);

EOF1 = fliplr(EOF1);

E1 = fliplr(flipud(E1));

EOFa = X*EOF1;

for i = 1:780

   EOF_k(:,i) = EOFa(:,i)/sqrt(E1(i,i));

end

EOF_1 = EOF_k;

PC_1 = EOF_k'*X;

tic

% 3 -----------------------------------------------------------------------

% EOF

[EOF2,E2] = eig(X*X');

PC2 = EOF2'*X;

% reverse the order

E2 = fliplr(flipud(E2));

% lambda = diag(E); % retain eigenvalues only

EOF_2 = fliplr(EOF2);

PC_2 = flipud(PC2);

% NOTE:

% (1) 方法1和3：程序差别在有没有除以总的时间序列长度（n），但是结果都一样，整个空间patern还是不变的

% (2) 方法2：当空间样本m远大于时间序列长度n时，计算m*m矩阵的特征根很困难，可以考虑对其进行时空转换

% (3) fliplr(flipud(E)) = rot90(E,2);

toc