%自己编的EOF程序
clear; clc
fid=fopen('f:\ga\pt\wave\jijie\pt.dat','r');
hgtn=fread(fid,'float');
aay=reshape(hgtn,246,12,58);
%a=mean(aay,3);
bb=aay(:,6,:)+aay(:,7,:)+aay(:,8,:);
ab=squeeze(bb);
%资料的标准化
an=ab';
avn=mean(an);
ssn=std(an,1);
n=58;
jp=an-repmat(avn,n,1);
bn=jp./repmat(ssn,n,1);
a=bn;
%提取主分量!
%先进行时空转换
s=a*a';
[v,d]=eig(s);
j=n; %将特征值按从大到小排序
for i=1:n
dd(i,i)=d(j,j);
j=j-1;
end
j=n; %将特征向量从大到小排序
for i=1:n
for k=1:n
vv(k,i)=v(k,j);
end
j=j-1;
end
sum=0.; %求累计方差
for i=1:n
sum=sum+dd(i,i);
end
for i=1:n %求方差贡献
c(i)=dd(i,i)/sum;
end
%求得eof的距阵
x1=a'*vv(:,1); %x1为第一特征向量空间场
tri=vv(:,1); %tri为时间系数
fid=fopen('f:\ga\pt\wave\eof\1000v.dat','w');
out=fwrite(fid,x1,'float');
fclose(fid);
fid=fopen('f:\ga\pt\wave\eof\ri.dat','w');
out=fwrite(fid,tri,'float');
fclose(fid);