%********************* 谐波叠加法模拟风速时程 修改1.0(采用Kaimal谱)**************

clc

clear

%************************* 风速时程参数设定 *************************

m=32;         %模拟风速点数

N=2^11;       %频率采样点数

dt=0.25;      %时间间隔

omegaup=3*pi; %上限频率

%************************  设定风速谱参数 *************************

L=514.5;      %斜拉桥跨度

z=45;         %风速测点离地面高度

z0=0.03;      %地面粗糙度

Uz=20;        %平均风速

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

K=0.4;        %Kaimal常数=0.4

M=2*N;        %nfft傅里叶变换长度，取采样个数的2倍

%形成风速时程矩阵

v=zeros(m,M*m);    %创建m行，M*m列的时程矩阵

u=zeros(m,M*m);

v1=zeros(M*m,m);   %创建m行，M*m列的时程矩阵

u2=zeros(M*m,m);

t=dt*(0:1:(M*m-1));%创建时程横坐标时间点

domega=(omegaup-0.001)/N;%频率间距

D=zeros(m,m,N);

U1=K*Uz/log(z/z0); %测点位置的摩阻速度

%形成目标谱

omega1=0.01:domega:omegaup;%形成频率列表，初始化频率列表，初始频率为omegaup/N

Sw1=200*U1^2.*z/Uz./(1+50.*omega1.*z./(2*pi*Uz)).^(5/3);%Kaimal谱密度表达式 水平向

Su1=3.36*U1^2.*z/Uz./(1+10.*omega1.*z./(2*pi*Uz)).^(5/3);%L-P谱密度表达式 竖直向

delta=12.9;

%模拟风速测点间的距离(第一个间隔不取为零，否则会出现S不正定的情况)

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*U1^2.*z/Uz./(1+50.*omega.*z./(2*pi*Uz)).^(5/3);

Su=3.36*U1^2.*z/Uz./(1+10.*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));%相关系数计算，采用Davenport形式。

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

U(j1,l,k)=Su(k)*Coh(j1,l,k);

end

end

end

%进行Cholesky分解

for i=1:1:N

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

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

Hu(:,:,i)=chol(U(:,:,i));

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

end

%填充谱数据矩阵D

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

E(:,j,:)=Hu(:,j,:);

i=sqrt(-1);

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

Bu1=sqrt(2*domega).*E(j,:,:);

for ii=1:j

for jj=1:N

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

Bu2(ii,jj)=Bu1(1,ii,jj);

end

end

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

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

for jj=1:j

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

Gu(jj,1:M)=fft(Bu2(jj,:),M);

for jjj=2:m

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

Gu(jj,((jjj-1)*M+1):(jjj*M))=Gu(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));

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

end

end

end %第一个for对应

v1=v';u1=u';

%风速时程绘图

figure(1)

subplot(3,1,1);

plot(t(1:1:4096),v(1,1:1:4096));

xlabel('t(s)');

ylabel('v(t)');

title('Point 1');

subplot(3,1,2);

plot(t(1:1:4096),v(2,1:1:4096));

xlabel('t(s)');

ylabel('v(t)');

title('Point 2');

subplot(3,1,3);

plot(t(1:1:4096),v(m,1:1:4096));

xlabel('t(s)');

ylabel('v(t)');

title('Point 32');

figure(2)

subplot(3,1,1);

plot(t(1:1:4096),u(1,1:1:4096));

xlabel('t(s)');

ylabel('u(t)');

title('Point 1');

subplot(3,1,2);

plot(t(1:1:4096),u(2,1:1:4096));

xlabel('t(s)');

ylabel('u(t)');

title('Point 2');

subplot(3,1,3);

plot(t(1:1:4096),u(m,1:1:4096));

xlabel('t(s)');

ylabel('u(t)');

title('Point 32');

[power1,freq1]=pwelch(v(1,:),hamming(3024),10,M*m,20);

[power2,freq2]=pwelch(v(2,:),hamming(3024),10,M*m,20);

[power5,freq5]=pwelch(v(5,:),hamming(3024),10,M*m,20);

%功率谱检验

figure(3)                                                             %第1点模拟自功率谱和计算自功率谱比较(对数坐标)

semilogx(freq1,freq1.*power1,'r',omega1,omega1.*Sw1,'b')

xlabel('Frequency w(rad)')

ylabel('Power wS(w)(rad.m2/s)')

title('Point 1')

figure(7)

loglog(freq5,freq5.*power5,'r',omega1,omega1.*Sw1,'b')

xlabel('Frequency w(rad)')

ylabel('Power wS(w)(rad.m2/s)')

title('Point 5')

figure(8)

plot(freq5,freq5.*power5,'r',omega1,omega1.*Sw1,'b')

xlabel('Frequency w(rad)')

ylabel('Power wS(w)(rad.m2/s)')

title('Point 5')