谐波叠加法 可运行程序，分享不谢

也要自己用。改变一天后，谐波叠加法可以操作程序。如果你使用它，复制过去。~~

采用%谐波叠加法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; Pm平均风速

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