资源描述：
《EOF分解程序说课材料.doc》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、fid=fopen('HadISST1_SST_1961-1990.txt','r');Num=360;data=zeros(360,180,Num);fori=1:Numaaa=fscanf(fid,'%s',7);data(:,:,i)=fscanf(fid,'%f',[360,180]);endsst1=data(1:90,11:70,1:Num);%选取所需要区域的数据sst2=data(311:360,11:70,1:Num);sst3=zeros(140,60,Num);sst3(90:-1:1,1:60,1:Num)=sst1;sst3(140:-1:91,1:60,1:Num
2、)=sst2;sst=sst3;fori=1:140forj=1:60fork=1:Numif(sst(i,j,k)==-1000)
3、
4、((sst(i,j,k)==-32768))sst(i,j,k)=NaN;endendendendsst_area1=zeros(Num,8400);%zeros全零数组fori=1:Num;squ=squeeze(sst(:,:,i));%执行该指令后sst数据转换为二维数组sst_area1(i,:)=reshape(squ,1,8400);%将数据转变为二维endsst_nan=isnan(sst_area1);i=0;forj=1:8400ifsum
5、(sst_nan(:,j))==0;i=i+1;sst_region(:,i)=sst_area1(:,j);endend%求距平~注意季节的变换X=zeros(size(sst_region));%学者给的程序fori=1:12X(i:12:Num-12+i,:)=sst_region(i:12:end,:)-repmat(mean(sst_region(i:12:end,:),1),size(sst_region(i:12:end,:),1),1);endR=X*X';%协方差矩阵R=X*X'是8400*8400的方阵~现定义矩阵R=X'*X是156*156的矩阵[v,d]=eig(R)
6、;%进行EOF分解~因为X'*X与X*X'的秩相同所以特征值相同~d为x的特征值组成的对角阵~v为X*X'的特征向量~[diagonal,I]=sort(diag(d),'descend');v=v(:,I);G=diagonal/sum(diagonal);Gs=0;fori=1:length(G)Gs(i)=sum(G(1:i));ifGs(i)>0.8break;endendN=length(Gs)%v=fliplr(v);%矩阵作左右翻转%d=rot90(d,2);%矩阵上下翻转后再左右翻转(查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr(flipud(d)
7、)=rot90(d,2))%diagonal=diag(d);spacef=X'*v;fori=1:Num;spacef(:,i)=spacef(:,i)/sqrt(diagonal(i));%空间本征函数endtimef=X*spacef;%时间本征函数sum_d=sum(diagonal);count=0;fori=1:Num;count=count+diagonal(i);G1(i)=count/sum_d;%G1(i)是累积方差贡献率endfori=1:Num;G2(i)=diagonal(i)/sum_d;%G2(i)是方差贡献率end%**********************
8、****************************************************%将删去的陆地与冰点的填充值补回sst_area2=zeros(Num,8400);sst_area2(:,:)=NaN;i=1;spacef2=spacef';forj=1:8400ifsum(sst_nan(:,j))==0;sst_area2(:,j)=spacef2(:,i);i=i+1;endendsst_area3=sst_area2';%**********************************************************************
9、****%画图%subplot(2,1,2)%将绘图窗口划分为2*1个子窗口，并在第2个子窗口中绘图figure(1)x=1:Num;plot(x,timef(:,1),'g');%ylim([-8080]);%xlabel('1980-1992年156个月','fontsize',15,'fontname','隶书')ylabel('INDEX','fontsize',12,'fontname','黑体')
