自己也要用到，改了一天，谐波叠加法可运行程序，大家用到的话就复制过去吧，不谢~~

%谐波叠加法,采用Kaimal谱

clc

clear

close all

%风速时程参数设置

m=10;

N=2^8;

dt=0.5;

omegaup=2*pi;

% 风速谱参数设定

L=1000;

z=50;

z0=0.03;      %地面粗糙度

Uz=40;        %50m处的平均风速

delta=100;    %模拟点间距

lambda=10;    %空间相关函数中的系数

K=0.4;

M=2*N;

v=zeros(m,M*m);

t=0.5*(0:1:(M*m-1));

domega=omegaup/N;

D=zeros(m,m,N);

U=K*Uz/log(z/z0);    %U为摩擦速度

%形成目标谱

omega1=omegaup/N:omegaup/N:omegaup;

Sw1=200*U^2.*z/Uz./(1+50.*omega1.*z./(2*pi*Uz)).^(5./3);

for j=1:m

rand('state',0);

thet=2*pi*rand(j,N);

for l=1:N  %%???

omega(l)=(l-1)*domega+j/m*domega;

end

Sw=200*U^2.*z/Uz./(1+50.*omega.*z./(2*pi*Uz)).^(5/3);

%计算谱数据库矩阵，功率谱计算，kaimal谱

for j1=1:m

for l=1:m

for k=1:N

Coh(j1,l,k)=(exp(-lambda*omega(k)*delta/(2*pi*Uz)))^(abs(j1-l));

S(j1,l,k)=Sw(k)*Coh(j1,l,k);

end

end

end

%Cholesky分解

for i=1:1:N

H(:,:,i)=chol(S(:,:,i));

H(:,:,i)=H(:,:,i)';

end

%填充谱数据矩阵D

D(:,j,:)=H(:,j,:);

i=sqrt(-1);

B1=sqrt(2*domega).*D(j,:,:);

for ii=1:j

for jj=1:N

B2(ii,jj)=B1(1,ii,jj);

end

end

B2=B2.*exp(i.*thet);

for jj=1:j

G(jj,1:M)=fft(B2(jj,:),M);

for jjj=2:m

G(jj,((jjj-1)*M+1):(jjj*M))=G(jj,1:M);

end

end

%谐波叠加生成模拟点的风速时程

for p=1:M*m

for k=1:j

v(j,p)=v(j,p)+real(G(k,p)*exp(i.*k./m.*domega.*(p-1).*dt));

end

end

end

% 显示风速时程

figure(1)

plot(t,v(1,:))

figure(2)

plot(t,v(5,:))

[power1,freq1]=psd(v(1,:),M*m,2,boxcar(512),0,'mean');

[power5,freq5]=psd(v(5,:),M*m,2,boxcar(512),0,'mean');

[power15,freq15]=cpsd(v(1,:),v(5,:),[],512,512,2);

%计算第一、第五点的模拟互功率谱

Sw15=[];

for i=1:N

s15=S(1,5,i);

Sw15=[Sw15,s15]

end

%功率谱检验

figure(3)

loglog(freq1,power1,'r',omega,Sw1,'b')

figure(4)

loglog(freq5,power5,'r',omega,Sw1,'b')

figure(5)

loglog(freq15,abs(power15),'r',omega,Sw15,'b')