%Importdata
fid=fopen('e:/uh/ci.dat','rb','ieee-le');
na=fread(fid,inf,'float');
na=na(1:444); %取特定列数组中的一段
nb=load('e:/uh/n.txt'); %读取txt文件
pra=load('e:/uh/ph_pr_are.txt'); %grads对某点降水异常做的面积平均
prep=ncread('e:/uh/precip.mon.mean.nc','precip');% 读取nc数据
prep=prep(:,:,1:444); %提取特定时间段
prep=reshape(prep,10368,12,37);
a=mean(prep,3); % 对全球144x72个格点分别做1-12月的37年平均
%pp是气候态,37years
pp=cat(3,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a,a);
prepa(:,:,:)=prep(:,:,:)-pp(:,:,:);% 全球降水异常
prepa1=reshape(prepa,144,72,444);
%根据corr选择相关最好的phase
nc=nb(:).*reshape(w(7:450),444,1);
%开始做回归%
X=[ones(size(nc)) nc na];
for i=1:144
for j=1:72
Y1=prepa1(i,j,1:444);
Y=reshape(Y1,444,1);
[b,bint,r,rint,stats] =regress(Y,X);
const(i,j)=b(1);
b_nc(i,j)=b(2);
b_na(i,j)=b(3);
end;
end;
%用回归系数算出重建场
for i=1:144
for j=1:72
for k=1:444
repra(i,j,k)=const(i,j)+b_nc(i,j)*nc(k)+b_na(i,j)*na(k);
end;
end;
end;
% 重建后的降水场
pp1=reshape(pp,144,72,444);
for i=1:144
for j=1:72
for k=1:444
repr(i,j,k)=repra(i,j,k)+pp1(i,j,k);
end;
end;
end;
%将数据输出dat文件格式
fid=fopen('e:/uh/repr1.dat','w');
for i=1:444
for j=1:72
for i=1:144
fwrite (fid,repra(i,j,k),'float32');
end;
end;
end;
fclose(fid);