資源描述：
《脈動(dòng)風(fēng)時(shí)程matlab程序.doc》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在教育資源-天天文庫(kù)。
1、根據(jù)風(fēng)的記錄，脈動(dòng)風(fēng)可作為高斯平穩(wěn)過程來考慮。觀察個(gè)具有零均值的平穩(wěn)高斯過程，其譜密度函數(shù)矩陣為：（9）將進(jìn)行Cholesky分解，得有效方法。（10）其中，（11）為的共軛轉(zhuǎn)置。根據(jù)文獻(xiàn)[8]，對(duì)于功率譜密度函數(shù)矩陣為的多維隨機(jī)過程向量，模擬風(fēng)速具有如下形式：（12）其中，風(fēng)譜在頻率范圍內(nèi)劃分成個(gè)相同部分，為頻率增量，為上述下三角矩陣的模，為兩個(gè)不同作用點(diǎn)之間的相位角，為介于和之間均勻分布的隨機(jī)數(shù)，是頻域的遞增變量。文中模擬開孔處的來流風(fēng)，因而只作單點(diǎn)模擬。即式（4）可簡(jiǎn)化為：（13）本文采用Davenp
2、ort水平脈動(dòng)風(fēng)速譜：（14）式中，脈動(dòng)風(fēng)速功率譜；脈動(dòng)風(fēng)頻率(Hz)；地面粗糙度系數(shù)；標(biāo)準(zhǔn)高度為10m處的風(fēng)速(m/s)。Matlab程序:N=10;d=0.001;n=d:d:N;%%頻率區(qū)間（0.01～10）v10=16;k=0.005;x=1200*n/v10;s1=4*k*v10^2*x.^2./n./(1+x.^2).^(4/3);%%Davenport譜subplot(2,2,1)loglog(n,s1)%%畫譜圖axis([-10015-1001000])xlabel('freq');yla
3、bel('S');fori=1:1:N/dH(i)=chol(s1(i));%%Cholesky分解endthta=2*pi*rand(N/d,1000);%%介于0和2pi之間均勻分布的隨機(jī)數(shù)t=1:1:1000;%%時(shí)間區(qū)間（0.1～100s）forj=1:1:1000a=abs(H);b=cos((n*j/10)+thta(:,j)');c=sum(a.*b);v(j)=(2*d).^(1/2)*c;%%風(fēng)荷載模擬endsubplot(2,2,2)plot(t/10,v)%%顯示風(fēng)荷載xlabel('
4、t(s)');ylabel('v(t)');Y=fft(v);%%對(duì)數(shù)值解作傅立葉變換Y(1)=[];%%去掉零頻量m=length(Y)/2;%%計(jì)算頻率個(gè)數(shù)；power=abs(Y(1:m)).^2/(length(Y).^2);%%計(jì)算功率譜freq=10*(1:m)/length(Y);%%計(jì)算頻率，因?yàn)椴介L(zhǎng)為0.1，而不是1，故乘以10subplot(2,2,3)loglog(freq,power,'r',n,s1,'b')%%比較axis([-10015-1001000])xlabel('fre
5、q');ylabel('S');對(duì)源程序的修改：z=xcorr(v);Y=fft(z);%%對(duì)數(shù)值解作傅立葉變換Y(1)=[];%%去掉零頻量m=length(Y)/2;%%計(jì)算頻率個(gè)數(shù)；power=abs(Y(1:m)).^2/(length(Y).^2);%%計(jì)算功率譜freq=10*(1:m)/length(Y);%%計(jì)算頻率，因?yàn)椴介L(zhǎng)為0.1，而不是1，故乘以10subplot(2,2,3)loglog(freq,power,'r',n,s1,'b')%%比較axis([-10015-1001000
6、])xlabel('freq');ylabel('S');樓主的修改使模擬得到的功率譜與源譜的數(shù)量級(jí)對(duì)上了，但是吻合不是太好。但是好像這樣做是不對(duì)的。????求信號(hào)x(t)的功率譜有兩種方法，一是對(duì)X(t)做傅立葉變換，再平方??????????S=abs(fft(x))^2????一是先對(duì)X(t)求相關(guān)系數(shù)，再進(jìn)行傅立葉變換:??????????S=fft(xcorr(X))樓主的方法好像是這兩個(gè)方法的混合。歡迎大家拍磚^_^
