根据风的记录，脉动风可作为高斯平稳过程来考虑。观察n个具有零均值的平稳高斯过

程，其谱密度函数矩阵为：

_Sii^)気临)...% (灼)］

?、S21(国)S22(⑷)...S2n(⑷)

S(CO)=

±1(00)乳儉)…Snn(G0)_

将SC )进行Cholesky分解，得有效方法。

S( ) = H( J

S( ) = H( J

H C )T

(10)

其中,

H (；：；■)=旳11心)|

H (；：；■)=

旳11心)

|H 21(豹)

H 22 C 0 ...

Bnlg)

Hn2( )...

Hnn?) 一

(11)

T

H C )为H (「)的共轭转置。

根据文献［8］，对于功率谱密度函数矩阵为SC )的多维随机过程向量， 模拟风速具有如

F形式：

j N

Vj ⑴=送 Z ‘H jm?)|cosblt 理 jm? )Pml ]

m=! l ±

j =1,2,3..., n(12)

其中，风谱在频率范围内划分成 N个相同部分，△⑷=⑷/N为频率增量，H jm (⑷丨)为 上述下三角矩阵的模，jm (打)为两个不同作用点之间的相位角，rml为介于0和2二之间

均匀分布的随机数， j =|是频域的递增变量。

文中模拟开孔处的来流风，因而只作单点模拟。即式(4)可简化为：

N

TOC \o "1-5" \h \z v(t)=送 |H (创)72M cosb|t+d】(13)

im

本文采用Davenport水平脉动风速谱：

24kx2

Sv(n)二 V102 473( 14)

n (1 x )

式中，Sv (n)――脉动风速功率谱; n——脉动风频率(Hz)； k——地面粗糙度系数；

x =1200_

V10

v10标准高度为10m处的风速(m/s)。

Matlab 程序：

N=10;

d=0.001;

n=d:d:N;%% 频率区间(0.01 ?10)

v10=16;

k=0.005;

x=1200* n/v10;

s1=4*k*v10A2*x.A2./n./(1+x42).A(4/3);%%Dave nport谱

subplot(2,2,1)

loglog(n,s1)%% 画谱图

axis([-100 15 -100 1000])

xlabel('freq');

ylabel('S');

for i=1:1:N/d

H(i)=chol(s1(i));%%Cholesky 分解

end

thta=2*pi*rand(N/d,1000);%% 介于0和2pi之间均匀分布的随机数 t=1:1:1000;%% 时间区间(0.1 ?100s)

for j=1:1:1000

a=abs(H);

b=cos(( n*j/10)+thta(:,j)');

c=sum(a.*b);

v(j)=(2*d).A(1/2)*c;%% 风荷载模拟

end subplot(2,2,2)

plot(t/10,v)%%显示风荷载

xlabel('t(s)');

ylabel('v(t)');

Y=fft(v);%%对数值解作傅立叶变换

Y(1)=[];%%去掉零频量

m=length(Y)/2;%%计算频率个数；

power=abs(Y(1:m)).A2/(le ngth(Y).A2);%% 计算功率谱

freq=10*(1:m)/le ngth(Y);%% 计算频率，因为步长为0.1，而不是1，故乘以10 subplot(2,2,3)Ioglog(freq,power,'r',n,s1,'b')%% 比较 axis([-100 15 -100 1000]) xlabel('freq');

ylabel('S');freqt(s)freq

10

10

10

10 freq

t(s)

freq

对源程序的修改:

z=xcorr(v);

Y=fft(z);%%对数值解作傅立叶变换

Y(1)=[];%%去掉零频量

m=length(Y)/2;%%计算频率个数；

power=abs(Y(1:m)).A2/(le ngth(Y).A2);%% 计算功率谱

freq=10*(1:m)/le ngth(Y);%% 计算频率，因为步长为 0.1，而不是1，故乘以10 subplot(2,2,3) loglog(freq,power,'r',n,s1,'b')%% 比较 axis([-100 15 -100 1000]) xlabel('freq');

ylabel('S');

楼主的修改使模拟得到的功率谱与源谱的数量级对上了，但是吻合不是太好。但是好像这样

做是不对的。

求信号x(t)的功率谱有两种方法，一是对X(t)做傅立叶变换，再平方

S=abs(fft(x))A2

一是先对X(t)求相关系数，再进