风振分析命令.doc

！画对比曲线n=0.001:0.001:0.5v10=34x=1200*n/v10s1=0.667*x.2./n./(1+x.2).(4/3) %规格化功率谱plot(n,s1)a=(x.2-0.77292).2+(2*1.477*0.7729*x).2b=800*2.182*x.2/v10s2=b./a %近似谱plot(n,s2)plot(n,s1,n,s2)plot(n,s1,k-,n,s2,b:)clearPij= 8.6719 8.9066 11.2417 8.3669 7.1518h=16.7 49.6 85 112.8 146for i=1:5 for j=1:5 Spij(i,j)=exp(-abs(h(i)-h(j)/60).*Pij(i).*Pij(j) endendgoal=chol(Spij)goal=goalN=4001afa=2.182Wf=0.1376kexi=1.477dw=0.01256V=34x=sqrt(-1)fai=2*pi*rand(1,N)for time=1:4001 t=(time-1)*0.05 for j=1:5 sum1=0 for m=1:j sum2=0 for k=1:N w=(k-1)*dw F1=sqrt(afa*pi*V/900)*w*x F2=(Wf2-w2)+2*kexi*Wf*w*x GS=F1/F2 GD(j,m)=sqrt(2*dw)*goal(j,m)*abs(GS) if k1 theta=atan2(imag(GS),real(GS) end Result=GD(j,m)*cos(w*t+theta+fai(k) sum2=sum2+Result end xie(j,m)=sum2 sum1=sum1+xie(j,m) end ft(time,j)=sum1 end ok(time)=tendfor i=1:4001 tu1(i)=ft(i,1) tu2(i)=ft(i,2) tu3(i)=ft(i,3) tu4(i)=ft(i,4) tu5(i)=ft(i,5)endn=0.001:0.001:10v10=34x=1200*n/v10a=(x.2-0.77292).2+(2*1.477*0.7729*x).2b=800*2.182*x.2/v10s2=b./a %近似谱c=2.182*v10*pi*pi*pi*n.2/225d=(4*pi*pi*n.2-0.13762).2+ 6.5226*n.2s3=c./dpxx,f=psd(tu1,1002,20,hanning(1002),501)f=f*0.2loglog(f,pxx,n,s2)clearPij= 8.6719 8.9066h=16.7 49.6 for i=1:2 for j=1:2 Spij(i,j)=exp(-abs(h(i)-h(j)/60).*Pij(i).*Pij(j) endendgoal=chol(Spij)goal=goalN=501afa=2.182Wf=0.1376kexi=1.477dw=0.01256V=34x=sqrt(-1)fai=2*pi*rand(1,N)for time=1:501 t=(time-1)*0.05 for j=1:2 sum1=0 for m=1:j sum2=0 for k=1:N w=(k-1)*dw F1=sqrt(afa*pi*V/900)*w*x F2=(Wf2-w2)+2*kexi*Wf*w*x GS=F1/F2 GD(j,m)=sqrt(2*dw)*goal(j,m)*abs(GS) if k1 theta=atan2(imag(GS),real(GS) end Result=GD(j,m)*cos(w*t+theta+fai(k) sum2=sum2+Result end xie(j,m)=sum2 sum1=sum1+xie(j,m) end ft(time,j)=sum1 end ok(time)=tendfor i=1:501 tu1(i)=ft(i,1) tu2(i)=ft(i,2) tu3(i)=ft(i,3) tu4(i)=ft(i,4) tu5(i)=ft(i,5)endclearPij= 8.6719 8.9066h=16.7 49.6 for i=1:2 for j=1:2 Spij(i,j)=exp(-abs(h(i)-h(j)/60).*Pij(i).*Pij(j) endendgoal=chol(Spij)goal=goalN=1001afa=2.182Wf=0.1376kexi=1.477dw=0.01256V=34x=sqrt(-1)fai=2*pi*rand(1,N)for time=1:1001 t=(time-1)*0.05 for j=1:2 sum1=0 for m=1:j sum2=0 for k=1:N w=(k-1)*dw F1=sqrt(afa*pi*V/900)*w*x F2=(Wf2-w2)+2*kexi*Wf*w*x GS=F1/F2 GD(j,m)=sqrt(2*dw)*goal(j,m)*abs(GS) if k1 theta=atan2(imag(GS),real(GS) end Result=GD(j,m)*cos(w*t+theta+fai(k) sum2=sum2+Result end xie(j,m)=sum2 sum1=sum1+xie(j,m) end ft