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

观察
个具有零均值的平稳高斯过程，其谱密度函数矩阵为：
（9）
将
进行Cholesky分解，得有效方法。

（10）
其中，
（11）
为
的共轭转置。

根据文献[8]，对于功率谱密度函数矩阵为
的多维随机过程向量，模拟风速具有如下形式：
（12）
其中，风谱在频率范围内划分成
个相同部分，
为频率增量，
为上述下三角矩阵的模，
为两个不同作用点之间的相位角，
为介于
和
之间均匀分布的随机数，
是频域的递增变量。

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

即式（4）可简化为：
（13）
本文采用Davenport水平脉动风速谱：
（14）
式中，
脉动风速功率谱；
脉动风频率(Hz)；
地面粗糙度系数；
标准高度为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*v10^2*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/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).^(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)).^2/(length(Y).^2);%%计算功率谱
freq=10*(1:m)/length(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');
对源程序的修改：
z=xcorr(v);
Y=fft(z);%%对数值解作傅立叶变换
Y(1)=[];%%去掉零频量
m=length(Y)/2;%%计算频率个数；
power=abs(Y(1:m)).^2/(length(Y).^2);%%计算功率谱
freq=10*(1:m)/length(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))^2
一是先对X(t)求相关系数，再进行傅立叶变换:
S=fft(xcorr(X))
楼主的方法好像是这两个方法的混合。

欢迎大家拍砖^_^
N=1024;
d=0.01;
n=d:d:N/100;
v10=16;
k=0.03;
x=1200*n/v10;
s1=4*k*v10^2*x.^2./n./(1+x.^2).^(4/3);%%Davenport谱for i=1:1:N
H(i)=chol(s1(i));%%Cholesky分解
end
rand('state',0);
thta=2*pi*rand(1,N);%%介于0和2pi之间均匀分布的随机数i=sqrt(-1);
B=H.*exp(i*thta);
G=fft(B,2*N);
for p=1:1:2*N
v(p)=2*sqrt(d)*real(G(p)*exp(i*2*pi*d*p*0.1/3));
end
[power,freq]=psd(v,1024*2,10,boxcar(1024),0,'mean'); power = power * 2 *0.1;%规一化修正
loglog(freq,power,'r',n,s1,'b')。