clear;close all;clc
tic
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
% 求距平
path_0 = '*****\';
path = [path_0,'sst\sst_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
评论0