matlab计算风速时程命令流

%clearall;N=500;%采样点数wu=4*2*pi;%截断频率dm=wu/N;%频率步长dt=3.2*pi/(2*wu);%时间步长0.2k=0.00464;%地面粗糙度系数地面粗糙度等级ABCD：K=0.001290.002150.004640.01291d=0.001;f=d:d:10;%时间从0.001到10s，步进值为0.001v10=28.2%设计风速28.2m/s——50年一遇十分钟平均风速最大值x=1200*f/v10;s=4*k*v10*v10.*x.^2./f./(1+x.^2).^(4/3);z1=10;%取第一点为10米高度z2=52.8;%取第二点为52.08米高度r=0.2;%考虑地表粗糙度影响的无量纲幂指数，按中国规范取0.22-c类v5=33.33%计算n米高处的平均风速——52.8m处平均风速期望值C=10;%指数衰减系数(取经验值)v1=zeros(2*N,1);%产生一个全零矩阵v2=zeros(2*N,1);thta1=rand(N,1);thta2=rand(N,1);%随机矩阵node=1;forK=1:nodeforj=1:2*Nsum1=0;sum2=0;forl=1:Nm1=l*dm-0.5*dm;m2=l*dm;x1=1200*m1/(2*pi*v10);s11=2*pi*4*k*v10*v10*x1*x1./m1./(1+x1*x1).^(4/3);x2=1200*m2/(2*pi*v10);s22=2*pi*4*k*v5*v5*x2*x2./m2./(1+x2*x2).^(4/3);s12=sqrt(s11*s22).*exp(-2*m2*C*abs(z1-z2)./(2*pi*(v10+v5)));s21=sqrt(s11*s22).*exp(-2*m1*C*abs(z1-z2)./(2*pi*(v10+v5)));S=[s11s12;s21s22];H=chol(S);%丘拉斯基分解-因式分解a1=abs(H(1,1));H1=H';a21=abs(H1(2,1));a22=abs(H1(2,2));b1=cos((m1*dt*(j-1))+2*pi*thta1(l,1));b2=cos((m2*dt*(j-1))+2*pi*thta2(l,1));c1=a1*b1;c21=a21*b1;c22=a22*b2;d1=(dm).^0.5*c1;d2=(dm).^0.5*(c21+c22);sum1=sum1+d1;sum2=sum2+d2;endsum1=0.8*sum1;sum2=0.8*sum2;v1(j,K)=sum1;v2(j,K)=sum2;endendu1=v1+v5u2=v2+v5t=(0:2*N-1)*dt;subplot(2,2,1);plot(t,v1,'b-');xlabel('t(s)');ylabel('v(t)');axis([-1180-3030]);holdonplot(t,v2,'r-');xlabel('t(s)');ylabel('v(t)');axis([-1180-3030]);Y=fft(v1);%对数值解作傅立叶变换Y(1)=[];%去掉零频量m=length(Y)/2;%计算频率个数；power=150*abs(Y(1:m)).^2/(length(Y).^2);%计算功率谱freq=5*(1:m)/length(Y);%计算频率，因为时间步长为0.125，而不是1，故乘以8subplot(2,2,2);loglog(freq,power,'r-',f,s,'g-');%对数显示，比较plot(t,v2,'r-');axis([-5015-101000]);%画出坐标轴比例axis([xminxmaxyminymax])set(gca,'xtick',[0.1110]);%自动定义刻度set(gca,'ytick',[0.1110]);xlabel('频率');ylabel('功率');