EOF分解程序

fid=fopen('','r');

Num=360;

data=zeros(360,180,Num);

for i=1:Num

aaa=fscanf(fid,'%s',7);

data(:,:,i)=fscanf(fid,'%f',[360,180]);

end

sst1=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)=sst2;

sst=sst3;

for i=1:140

for j=1:60

for k=1:Num

if(sst(i,j,k)==-1000)||((sst(i,j,k)==-32768))

sst(i,j,k)=NaN;

end

end

end

end

sst_area1=zeros(Num,8400); % zeros全零数组

for i=1:Num;

squ=squeeze(sst(:,:,i)); % 执行该指令后sst数据转换为二维数组

sst_area1(i,:)=reshape(squ,1,8400); % 将数据转变为二维

end

sst_nan=isnan(sst_area1);

i=0;

for j=1:8400

if sum(sst_nan(:,j))==0;

i=i+1;

sst_region(:,i)=sst_area1(:,j);

end

end

% 求距平~注意季节的变换

X=zeros(size(sst_region)); % 学者给的程序

for i=1:12

X(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);

end

R=X*X'; % 协方差矩阵R=X*X'是8400*8400的方阵~现定义矩阵R=X'*X是156*156的矩阵

[v,d]=eig(R); % 进行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;

for i=1:length(G)

Gs(i)=sum(G(1:i));

if Gs(i)> break;

end

end

N=length(Gs)

%v=fliplr(v); % 矩阵作左右翻转

%d=rot90(d,2); % 矩阵上下翻转后再左右翻转(查看生成的对角阵是由小到大排列的~此指令可使其由大到小排列~fliplr(flipud(d))=rot90(d,2))

%diagonal=diag(d);

spacef=X'*v;

for i=1:Num;

spacef(:,i)=spacef(:,i)/sqrt(diagonal(i)); % 空间本征函数

end

timef=X*spacef; % 时间本征函数

sum_d=sum(diagonal);

count=0;

for i=1:Num;

count=count+diagonal(i);

G1(i)=count/sum_d; % G1(i)是累积方差贡献率

end

for i=1:Num;

G2(i)=diagonal(i)/sum_d; % G2(i)是方差贡献率

end

%********************************************************** ****************

% 将删去的陆地与冰点的填充值补回

sst_area2=zeros(Num,8400);

sst_area2(:,:)=NaN;

i=1;

spacef2=spacef';

for j=1:8400

if sum(sst_nan(:,j))==0;

sst_area2(:,j)=spacef2(:,i);

i=i+1;

end

end

sst_area3=sst_area2';

%**********************************************************

**************** % 画图

% subplot(2,1,2) % 将绘图窗口划分为2*1个子窗口，并

在第2个子窗口中绘图

figure(1)

x=1:Num;

plot(x,timef(:,1),'g');

%ylim([ -80 80 ]);

% xlabel('1980-1992年156个月','fontsize',15,'fontname','隶书')

ylabel('INDEX','fontsize',12,'fontname','黑体')

set(gca,'xtick',(1:6:Num))

set(gca,'xticklabel',{'1980','','1981','','1982','','1983', '','1984','','1985','','1986','','1987','','1988','','1989' ,'','1990','','1991','','1992','','1993'})

title('北太平洋第1模态1980至1992年SST时间序列', 'color', 'k','fontsize',15,'fontname','幼圆')

grid on

hold off

% %subplot(2,1,1)

% lon=[:];

% lat=[:];

% m_proj('Equidistant Cylindrical','lat',[ ],'lon',[ ]);

%