搜索
%谐波叠加法模拟风速时程（采用Kaimal谱）clcclear%风速时程参数设定m=10;%模拟风速个数N=2^8;%频率采样点数dt=0.5;%时间间隔omegaup=2*pi;%上限频率%风速谱参数设定L=1000;%跨度z=50;%模拟点离地面高度z0=0.03;%地面粗糙度Uz=40;%50m处的平均风速delta=100;%模拟点间距(m)lambda=10;%空间相关函数中的系数K=0.4;%Kaman常数0.4M=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);forj=1:m%对m条风速循环rand('state',0);%模拟随机数序列thet=2*pi*rand(j,N);%生成随机相位forl=1:N%生成频率序列omega(l)=(l-1)*domega+j/m*domega;endSw=200*U^2.*z/Uz./(1+50.*omega.*z./(2*pi*Uz)).^(5/3);%计算数据库矩阵，功率谱计算，Kaimal谱forj1=1:mforl=1:mfork=1:NCoh(j1,l,k)=(exp(-lambda*omega(k)*delta/(2*pi*Uz)))^(abs(j1-l));S(j1,l,k)=Sw(k)*Coh(j1,l,k);endendend%Cholsky分解fori=1:1:NH(:,:,i)=chol(S(:,:,i));H(:,:,i)=H(:,:,i)';end%填充谱数据矩阵DD(:,j,:)=H(:,j,:);i=sqrt(-1);B1=sqrt(2*domega).*D(j,:,:);forii=1:jforjj=1:NB2(ii,jj)=B1(1,ii,jj);endendB2=B2.*exp(i.*thet);forjj=1:jG(jj,1:M)=fft(B2(jj,:),M);forjjj=2:mG(jj,((jjj-1)*M+1):(jjj*M))=G(jj,1:M);endend%谐波叠加生成模拟点的风速时程forp=1:M*mfork=1:jv(j,p)=v(j,p)+real(G(k,p)*exp(i.*k./m.*domega.*(p-1).*dt));endendend%显示风速时程figure(1)%第1点风速时程plot(t,v(1,:))figure(2)%第5点风速时程plot(t,v(5,:))[power1,freq1]=psd(v(1,:),M*m,2,boxcar(512),0,'mean');%计算第1点的模拟自功率谱[power5,freq5]=psd(v(5,:),M*m,2,boxcar(512),0,'mean');%计算第5点的模拟自功率谱[power15,freq15]=cpsd(v(1,:),v(5,:),[],512,512,2);%计算第1,5点的目标互功率谱Sw15=[];%第1,5点的目标互功率谱fori=1:Ns15=S(1,5,i);Sw15=[Sw15,s15];end%功率谱检验figure(3)%第1点模拟自功率谱和计算自功率谱比较（对数坐标）loglog(freq1,power1,'r',omega1,Sw1,'b')figure(4)%第5点模拟自功率谱和计算自功率谱比较（对数坐标）loglog(freq5,power5,'r',omega1,Sw1,'b')figure(5)%第1,5点模拟自功率谱和计算互功率谱比较（一般坐标）plot(freq15,abs(power15),'r',omega1,Sw15,'b')