EOF.zip_EOF；经验正交分解_matlab eof_正交分解_气象_海洋

clear; close all;clc; load('75-13.mat'); [m,it]=size(Untitled); X=Untitled M=5; onoff=0; % method=0 原始场 % 1 距平 % 2 中心化 % 3 对数中心化 % 4 标准化 % 5 均值化 method = 4; % n=length(x(1,:)); % m=length(x(:,1)); [m,n] = size(X); % for iii=1:m % meana=mean(x(iii,:)); % mean value of each column % x(iii,:)=x(iii,:)-meana; %算出每行的距平值 % end switch method case 0 case 1 %距平 mF=mean(X); X=X-mF(ones(m,1),:); case 2 %中心化 for iii=1:m X(iii,:)=X(iii,:)-mean(X(iii,:)); end case 3 %对数中心化 F=log(F); for iii=1:m X(iii,:)=X(iii,:)-mean(X(iii,:)); end case 4 %标准化 for iii=1:m X(iii,:)=(X(iii,:)-mean(X(iii,:)))/std(X(iii,:)); end case 5 %均值化 for iii=1:m X(iii,:)=X(iii,:)/mean(X(iii,:)); end otherwise end sprintf('m=%d n=%d',m,n) if m>n S=X'*X/n; sub_x=X/sqrt(n); onoff=1; else S=X*X'/n; end [L,D]=eig(S,'nobalance'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R=X*X'/39; % 协方差矩阵R=X*X'是11280*11280的方阵~现% 要旋转定义矩阵R=X'*X是405*405的矩阵 这里没有旋转 [v,d]=eig(R); % 进行EOF分解~因为X'*X与X*X'的秩相同，所以特征值相同。d为x的特征值组成的对角阵。v为X*X'的特征向量~ PCA=v'*X % 时间本征函数 D=rot90(d,2); % 矩阵上下翻转后再左右翻转(查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr(flipud(d))=rot90(d,2) diagonal=diag(D); % 只保留特征值 EOF=fliplr(v) % 矩阵左右翻转 PCA=flipud(PCA); % 矩阵作上下翻转 % 取第一二三模态时间序列%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PCA1=PCA(1,:);PCA2=PCA(2,:);PCA3=PCA(3,:); % 取第一二三模态对应的空间向量~并写成格点格式ix240*iy47=11280%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f1=EOF(:,1); eof1=reshape(f1,28,1); f2=EOF(:,2); eof2=reshape(f2,28,1); f3=EOF(:,3); eof3=reshape(f3,28,1); %%%%%%%%%%xx=eof*pc; %检查eof的正交性 ae=EOF*EOF'; %检查pc的正交性 ap=PCA*PCA'; %回复到x的距平值 XX=EOF*PCA; %以上说明计算数据正确 % 回复到winter_juping的距平值 说明计算正确 % 求各模态的方差贡献率%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% sum_d=sum(diagonal); % 总的特征值 e1=(diagonal(1,1)/sum_d)*100%; %38.2920% e2=(diagonal(2,1)/sum_d)*100%; %9.2266% e3=(diagonal(3,1)/sum_d)*100%; %9.1010% %总的方差贡献率%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% count=0; for i=1:m; count=count+diagonal(i); G1(i)=count/sum_d; % G1(i)是累积方差贡献率 end for i=1:m; G2(i)=diagonal(i)/sum_d; % G2(i)是方差贡献率 end