matlab程序.docx

1、matlab程序clearclcfid=fopen(CO2.txt,r);tline=fgetl(fid);co=fscanf(fid,%f,13,57);Fa=2;n=1; for i=1:1:53 %把矩阵改写为时间序列 if i=1 for j=4:1:13 tco(n)=co(j,i); t(n)=co(1,i)+j./12; n=n+1; end else for j=2:1:13 tco(n)=co(j,i); t(n)=co(1,i)+j./12; n=n+1; end endend %趋势项处理X=ones(n-1,2);X(:,2)=t;b=regress(tco,X);p2

2、=polyfit(t,tco,2);for i=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); %二次拟合endfigure(1)plot(t,tco)title(CO_2浓度的时间变化曲线,fontsize,15)xlabel(Time,fontsize,15)ylabel(CO_2浓度/ppm,fontsize,15)hold on%plot(t,q1,linestyle,-,color,1 0 0)plot(t,q2,linestyle,-,color,1 0 0)legend(实测曲线,趋

3、势线,location,northwest)legend(boxoff)%周期项处理%这里用谐波分析，假设T=634，这样处理有问题但结果没有受到太大影响for i=1:1:n-1 c(i)=tco(i)-q2(i); %剔除趋势项endcc=smooth(c,12); % 12点平滑a0=mean(c);w=2*pi/(n-1); for i=1:1:floor(n-1)/2) for j=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,:)

4、; s(i)=(a(i).2+b(i).2)/2; %计算方差贡献endi=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;while kkFa kk=kk+1; else break endendfor i=1:1:n-1 p(i)=a0; for j=1:1:kk-1 p(i)=p(i)+a(I(j)*cos(I(j)*w*i)+b(I(j)*

5、sin(I(j)*w*i); %周期项拟合 ppppp=(n-1)/I(j) %把波数变为周期 endendfigure(3)plot(t,c)xlabel(Time,fontsize,15)ylabel(CO_2浓度/ppm,fontsize,15)hold onplot(t,cc,linestyle,-,color,0 0 0)plot(t,p,linestyle,-,color,1 0 0)legend(剔除趋势项后的变化曲线,12点平滑后的曲线,拟合曲线,location,northeast)legend(boxoff)%随机项处理for i=1:1:n-1 sj(i)=c(i)-p(

6、i); %剔除周期项endu=mean(sj);for i=1:1:n-1 sj(i)=sj(i)-u; %得到平稳距平序列endfigure(4)i=1:1:n-1;plot(i,sj)title(平稳随机序列曲线,fontsize,15)xlabel(Time,fontsize,15)ylabel(CO_2浓度/ppm,fontsize,15)for to=1:1:floor(n/4) K(to)=0; for i=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(

7、to)=K(to)/(var(sj); %计算相关函数，to是时间差endfor i=1:1:floor(n/4) for j=1:1:floor(n/4) if i=j xro(i,j)=1; elseif idata2(1,j)&(height(i)data2(1,43)|(j=1)&(height(i)v(max)&(k(i)=1) max=i; end if(min=1)&(k(i)=0)&(k(1)=1)|(v(i)v(min)&(k(i)=0) min=i; endendif(l3)&(L+1F) disp( 引入第, num2str(max), 个变量); k(max)=0; L

8、=L+1; l=l+1; r=matdel(max,m+1,r); q=1; endelse F2=v(min)/(r(m+1,m+1)/(n-l-1); if(F2F) disp( 引入第, num2str(max), 个变量); k(max)=0; L=L+1; l=l+1; r=matdel(max,m+1,r); q=1; end endendenddisp(没有可剔除或引入的变量，逐步回归结束);a=zeros(L);j=1;for i=1:m if (k(i)=0) a(j)=i; j=j+1; end;end;xx=x(:,a(1);for i=2:Lxx=xx x(:,a(i)

9、;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);for i=2:L ss=ss,num2str(a(i);endss=ss,个变量;ss1=y=,num2s

10、tr(b(1),+（,num2str(b(2),x,num2str(a(1),);for i=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,%sn,ss);fprintf(fid,%sn,ss1);fprintf(fid,%sn,ss2);fprintf(fid,%sn,ss3);fprintf(fid,%sn,ss4);fprintf(fid

11、,%s,ss5);fclose(fid);end clear;fid=fopen(D:gongjumatlabmfiledicengshuju.txt,r);data=fscanf(fid,%f,27,148);t=data(1,:);for k=1:26plot(t,data(k+1,:);switch k case 1 title(Pinus); case 2 title(Tsuga); case 3 title(Picea); case 4 title(Betula); case 5 title(Alnus); case 6 title(Carpinus); case 7 title(Corylus); case 8 title(Quercus（落）endsaveas(gcf,D:gongjumatlabmfiletutu,num2str(k),.jpg);end;