脈動風(fēng)時程matlab程序

1、精選優(yōu)質(zhì)文檔-傾情為你奉上根據(jù)風(fēng)的記錄，脈動風(fēng)可作為高斯平穩(wěn)過程來考慮。觀察個具有零均值的平穩(wěn)高斯過程，其譜密度函數(shù)矩陣為： （9）將進行Cholesky分解，得有效方法。 （10）其中， （11）為的共軛轉(zhuǎn)置。根據(jù)文獻8，對于功率譜密度函數(shù)矩陣為的多維隨機過程向量，模擬風(fēng)速具有如下形式： （12）其中，風(fēng)譜在頻率范圍內(nèi)劃分成個相同部分，為頻率增量，為上述下三角矩陣的模，為兩個不同作用點之間的相位角，為介于和之間均勻分布的隨機數(shù)，是頻域的遞增變量。文中模擬開孔處的來流風(fēng)，因而只作單點模擬。即式（4）可簡化為： （13）本文采用Davenport水平脈動風(fēng)速譜： （14） 式中，脈動風(fēng)速功率譜；

2、脈動風(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.0110）v10=16;k=0.005;x=1200*n/v10;s1=4*k*v102*x.2./n./(1+x.2).(4/3);%Davenport譜subplot(2,2,1)loglog(n,s1)%畫譜圖axis(-100 15 -100 1000)xlabel('freq');ylabel('S');for i=1:1:N/dH(i)=chol(s1(i);%Cholesky分解endthta=2*pi

3、*rand(N/d,1000);%介于0和2pi之間均勻分布的隨機數(shù)t=1:1:1000;%時間區(qū)間（0.1100s）for j=1:1:1000 a=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('t(s)');ylabel('v(t)');Y=fft(v);%對數(shù)值解作傅立葉變換Y(1)=;%去掉零頻量m=length(Y)/2;%計算頻率個數(shù)；power=abs(Y(

4、1:m).2/(length(Y).2);%計算功率譜freq=10*(1:m)/length(Y);%計算頻率，因為步長為0.1，而不是1，故乘以10subplot(2,2,3)loglog(freq,power,'r',n,s1,'b')%比較axis(-100 15 -100 1000)xlabel('freq');ylabel('S');對源程序的修改：z=xcorr(v);Y=fft(z);%對數(shù)值解作傅立葉變換Y(1)=;%去掉零頻量m=length(Y)/2;%計算頻率個數(shù)；power=abs(Y(1:m).2/(l

5、ength(Y).2);%計算功率譜freq=10*(1:m)/length(Y);%計算頻率，因為步長為0.1，而不是1，故乘以10subplot(2,2,3)loglog(freq,power,'r',n,s1,'b')%比較axis(-100 15 -100 1000)xlabel('freq');ylabel('S');樓主的修改使模擬得到的功率譜與源譜的數(shù)量級對上了，但是吻合不是太好。但是好像這樣做是不對的。     求信號x(t)的功率譜有兩種方法，一是對X(t)做傅立葉變換，再平方           S=abs(fft(x)2     一是先對X(t)求相關(guān)系數(shù)，再進行傅立葉變換: