脉动风时程matlab程序

根据风的记录，脉动风可作为高斯平稳过程来考虑。

观察n 个具有零均值的平稳高斯过 程，其谱密度函数矩阵为：

_Sii ^）気临）...% （灼）］

〜、

S 21（国）S 22（⑷）...S 2n （⑷）

S （CO ）= ±1（00）乳儉）…Snn （G0）_

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

其中,

T

H C ）为H （「）的共轭转置。

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

j N

V j ⑴=送 Z ‘H jm ®)|

cosb l t 理 jm ® )P ml ]

m=! l ± j =1,2,3..., n （12）

其中，风谱在频率范围内划分成 N 个相同部分，△⑷=⑷/N 为频率增量，H jm （⑷丨）为 上述

下三角矩阵的模，jm （打）为两个不同作用点之间的相位角， r ml 为介于0和2二之间

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

文中模拟开孔处的来流风，因而只作单点模拟。即式（

4）可简化为：

N v （t ）=送 |H （创）72M cos b |t +d 】

（13）

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

2 4kx 2

S v （n ）二 V 10 2 473

（ 14） n （1 x ）

式中，S v （n ）――脉动风速功率谱; n ——脉动风频率（Hz ）； k ——地面粗糙度系数；

S( ) = H( J H C )T (10)

H （；：；■）= 旳11心）

|H

21（豹）

H

22 C 0 ... Bnlg) H n2( ) ... H nn®) 一 (11)

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');

对源程序的修改:

z=xcorr(v);

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

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

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

power=abs(Y(1:m)).A 2/(le ngth(Y).A

2);%% 计算功率谱 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)求相关系数，再进行傅立叶变换 ：

10

freq 10

10

10

t(s)

freq

S=fft(xcorr(X))

楼主的方法好像是这两个方法的混合。欢迎大家拍砖

freq

i