matlab程序.docx
- 文档编号:9421427
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:14
- 大小:17.08KB
matlab程序.docx
《matlab程序.docx》由会员分享,可在线阅读,更多相关《matlab程序.docx(14页珍藏版)》请在冰豆网上搜索。
matlab程序
clear
clc
fid=fopen('CO2.txt','r');
tline=fgetl(fid);
co=fscanf(fid,'%f',[13,57]);
Fa=2;
n=1;
fori=1:
1:
53%把矩阵改写为时间序列
ifi==1
forj=4:
1:
13
tco(n)=co(j,i);
t(n)=co(1,i)+j./12;
n=n+1;
end
else
forj=2:
1:
13
tco(n)=co(j,i);
t(n)=co(1,i)+j./12;
n=n+1;
end
end
end
%%%
%%趋势项处理
X=ones(n-1,2);
X(:
2)=t';
b=regress(tco',X);
p2=polyfit(t,tco,2);
fori=1:
1:
n-1
q1(i)=b
(1)+b
(2)*t(i);%线性拟合
q2(i)=p2
(1)*t(i).^2+p2
(2)*t(i)+p2(3);%二次拟合
end
figure
(1)
plot(t,tco)
title('CO_{2}浓度的时间变化曲线','fontsize',15)
xlabel('Time','fontsize',15)
ylabel('CO_{2}浓度/ppm','fontsize',15)
holdon
%plot(t,q1,'linestyle','--','color',[100])
plot(t,q2,'linestyle','-','color',[100])
legend('实测曲线','趋势线','location','northwest')
legend('boxoff')
%%%
%%周期项处理
%%这里用谐波分析,假设T=634,这样处理有问题但结果没有受到太大影响
fori=1:
1:
n-1
c(i)=tco(i)-q2(i);%剔除趋势项
end
cc=smooth(c,12);%12点平滑
a0=mean(c);
w=2*pi/(n-1);
fori=1:
1:
floor((n-1)/2)
forj=1:
1:
n-1
ccc(i,j)=cos(i*w*j);
sss(i,j)=sin(i*w*j);
end
a(i)=2*mean(c.*ccc(i,:
));
b(i)=2*mean(c.*sss(i,:
));
s(i)=(a(i).^2+b(i).^2)/2;%计算方差贡献
end
i=1:
1:
floor((n-1)/2);
figure
(2)
plot(i,s)
title('谐波谱图','fontsize',15)
xlabel('波数','fontsize',15)
ylabel('方差贡献','fontsize',15)
[news,I]=sort(s,'descend');%方差贡献从小到大排列
F=0;
kk=1;
whilekk<=length(s)%找显著波数
F=s(kk)*(n-4)/2/(var(c)-s(kk));
ifF>Fa
kk=kk+1;
else
break
end
end
fori=1:
1:
n-1
p(i)=a0;
forj=1:
1:
kk-1
p(i)=p(i)+a(I(j))*cos(I(j)*w*i)+b(I(j))*sin(I(j)*w*i);%周期项拟合
ppppp=(n-1)/I(j)%把波数变为周期
end
end
figure(3)
plot(t,c)
xlabel('Time','fontsize',15)
ylabel('CO_{2}浓度/ppm','fontsize',15)
holdon
plot(t,cc,'linestyle','-','color',[000])
plot(t,p,'linestyle','--','color',[100])
legend('剔除趋势项后的变化曲线','12点平滑后的曲线','拟合曲线','location','northeast')
legend('boxoff')
%%%
%%随机项处理
fori=1:
1:
n-1
sj(i)=c(i)-p(i);%剔除周期项
end
u=mean(sj);
fori=1:
1:
n-1
sj(i)=sj(i)-u;%得到平稳距平序列
end
figure(4)
i=1:
1:
n-1;
plot(i,sj)
title('平稳随机序列曲线','fontsize',15)
xlabel('Time','fontsize',15)
ylabel('CO_{2}浓度/ppm','fontsize',15)
forto=1:
1:
floor(n/4)
K(to)=0;
fori=1:
1:
n-1-to
K(to)=K(to)+(sj(i)-mean(sj))*(sj(i+to)-mean(sj));
end
K(to)=K(to)/(n-1-to);
ro(to)=K(to)/(var(sj));%计算相关函数,to是时间差
end
fori=1:
1:
floor(n/4)
forj=1:
1:
floor(n/4)
ifi==j
xro(i,j)=1;
elseifi xro(i,j)=ro(j-i); else xro(i,j)=ro(i-j); end end yro(i)=ro(i); end fork=1: 1: floor(n/4) fi=xro(1: k,1: k)\yro(1: k)'; FPE(k)=0; fori=1: 1: k FPE(k)=FPE(k)+fi(i)*ro(i); end FPE(k)=(1-FPE(k))*(n+k)/(n-k-2); end [newFPE,I2]=sort(FPE,'descend'); zyp=I2(length(FPE));%定阶,可能有问题 fi=xro(1: zyp,1: zyp)\yro(1: zyp)';%计算自回归模型的系数 %%% %%预报 nn=n; fori=54: 1: 57%把2011-2014的数据改写为序列形式 ifi==57 forj=2: 1: 12 tco(nn)=co(j,i); t(nn)=co(1,i)+j./12; nn=nn+1; end else forj=2: 1: 13 tco(nn)=co(j,i); t(nn)=co(1,i)+j./12; nn=nn+1; end end end fori=n: 1: nn-1%计算随机项的拟合值 sj(i)=0; forj=i-zyp: 1: i-1 sj(i)=sj(i)+fi(i-j)*sj(j); end sj(i)=sj(i)+u; end fori=n: 1: nn-1 p(i)=a0; forj=1: 1: kk-1 p(i)=p(i)+a(I(j))*cos(I(j)*w*i)+b(I(j))*sin(I(j)*w*i); end fco(i)=p2 (1)*t(i).^2+p2 (2)*t(i)+p2(3)+p(i)+sj(i);%预报结果 end figure(5) scatter(tco(n: nn-1),fco(n: nn-1)) title('实测值与拟合值散点图','fontsize',15) xlabel('实测值/ppm','fontsize',15) ylabel('拟合值/ppm','fontsize',15) xx=tco(n: nn-1); yy=fco(n: nn-1); r=corrcoef(xx',yy')%计算相关系数 figure(6) plot(t(n: nn-1),tco(n: nn-1)) holdon plot(t(n: nn-1),fco(n: nn-1),'linestyle','--','color',[100]) title('实测曲线与拟合曲线','fontsize',15) xlabel('Time','fontsize',15) ylabel('CO_{2}浓度/ppm','fontsize',15) legend('实测曲线','拟合曲线','location','northwest') legend('boxoff') 花粉 fid=fopen('D: \gongju\matlab\mfile\dibiaoshuju.txt','r'); data=fscanf(fid,'%f',[10,74]); data=data'; fclose(fid); fid=fopen('D: \gongju\matlab\mfile\ttt.txt','r'); data3=fscanf(fid,'%f',[5,148]); fclose(fid); time=data3(1,: )'; x=data3(2: 5,: )'; x=[ones(148,1)x]; b=[10.1093-0.011847-0.079571-0.035224-0.015373]'; y=x*b; b1=[1456.25790.370572.44331.15970.45866]'; y1=x*b1; plot(time,y); title('温度变化'); saveas(gcf,['D: \gongju\matlab\mfile\tu\temperature.jpg']); plot(time,y1); title('降水量变化'); saveas(gcf,['D: \gongju\matlab\mfile\tu\rainfall.jpg']); fid=fopen('D: \gongju\matlab\mfile\height.txt','r'); height=fscanf(fid,'%f'); fclose(fid); fid=fopen('D: \gongju\matlab\mfile\tp.txt','r'); data2=fscanf(fid,'%f',[3,43]); fclose(fid); fori=1: 74 forj=1: 42 if((height(i)>data2(1,j))&&(height(i)<=data2(1,j+1)))||((j==42)&&(height(i)>data2(1,43)))||((j==1)&&(height(i)<=data2(1,1))) t(i)=data2(2,j)+(height(i)-data2(1,j))/(data2(1,j+1)-data2(1,j))*(data2(2,j+1)-data2(2,j)); p(i)=data2(3,j)+(height(i)-data2(1,j))/(data2(1,j+1)-data2(1,j))*(data2(3,j+1)-data2(3,j)); end end end t=t'; p=p'; 逐步回归 functionstepregress(x,y,F) %x=zscore(x,1); %y=zscore(y,1); r=corrcoef([x,y]); l=0; L=0; [n,m]=size(x); k=ones(m); q=1; while(q==1) q=0; fori=1: m v(i)=r(i,m+1)^2/r(i,i); end max=1; min=1; fori=1: m if((max==1)&&(k(i)==1)&&(k (1)==0))||((v(i)>v(max))&&(k(i)==1)) max=i; end if((min==1)&&(k(i)==0)&&(k (1)==1))||((v(i) min=i; end end if(l<3)&&(L+1<=m) F1=v(max)/((r(m+1,m+1)-v(max))/(n-l-2)); if(F1>F) disp(['引入第',num2str(max),'个变量']); k(max)=0; L=L+1; l=l+1; r=matdel(max,m+1,r); q=1; end else F2=v(min)/(r(m+1,m+1)/(n-l-1)); if((F2 disp(['剔除第',num2str(min),'个变量']); k(min)=1; L=L-1; l=l+1; r=matdel(min,m+1,r); q=1; else F1=v(max)/((r(m+1,m+1)-v(max))/(n-l-2)); if(F1>F) disp(['引入第',num2str(max),'个变量']); k(max)=0; L=L+1; l=l+1; r=matdel(max,m+1,r); q=1; end end end end disp('没有可剔除或引入的变量,逐步回归结束'); a=zeros(L); j=1; fori=1: m if(k(i)==0) a(j)=i; j=j+1; end; end; xx=x(: a (1)); fori=2: L xx=[xxx(: a(i))]; end; xx=[ones(n,1)xx]; b=regress(y,xx);%回归系数 R=sqrt(1-r(m+1,m+1));%复相关系数 yyy=xx*b;%y的估计值 ymean=mean(y);%y平均值 Q=(y-yyy)'*(y-yyy);%剩余平方和 U=(yyy-ymean)'*(yyy-ymean);%回归平方和 rs=Q/(n-L-1);%剩余方差 f=U/L/(Q/(n-L-1));%F统计量 fid=fopen('result','w'); ss=['引入第',num2str(a (1))]; fori=2: L ss=[ss,',',num2str(a(i))]; end ss=[ss,'个变量']; ss1=['y=',num2str(b (1)),'+(',num2str(b (2)),'x',num2str(a (1)),')']; fori=2: L ss1=[ss1,'+(',num2str(b(i+1)),'x',num2str(a(i)),')']; end; ss2=['复相关系数=',num2str(R)]; ss3=['剩余方差=',num2str(rs)]; ss4=['F统计量=',num2str(f)]; ss5=['剩余平方和=',num2str(Q)]; fprintf(fid,'%s\n',ss); fprintf(fid,'%s\n',ss1); fprintf(fid,'%s\n',ss2); fprintf(fid,'%s\n',ss3); fprintf(fid,'%s\n',ss4); fprintf(fid,'%s',ss5); fclose(fid); end clear; fid=fopen('D: \gongju\matlab\mfile\dicengshuju.txt','r'); data=fscanf(fid,'%f',[27,148]); t=data(1,: ); fork=1: 26 plot(t,data(k+1,: )); switchk case1 title('Pinus'); case2 title('Tsuga'); case3 title('Picea'); case4 title('Betula'); case5 title('Alnus'); case6 title('Carpinus'); case7 title('Corylus'); case8 title('Quercus(落)' end saveas(gcf,['D: \gongju\matlab\mfile\tu\tu',num2str(k),'.jpg']); end;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 程序