谐波合成法-matlab程序文件

m=10;                   %模拟风速个数
N=2^8;                   %频率采样个数？？？？？？？为了在计算中使用FFT技术，一般取值为2^α，α为正整数
dt=0.5 ;                   %时间间隔
omegaup=2*pi;                %上限频率
%风速谱参数设定
L=1000;                         %跨度
z=50;                           %模拟点离地面高度
z0=0.03;                          %地面粗糙度  默认取值0.03
Uz=40;                          %50米处的平均风速
delta=100;                      %模拟点间距（m）
lambda=10;                      %空间相关函数中的系数     水平向衰减因子
K=0.4;                          %Kaman常数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为摩阻速速  气流的剪切速度
%形成目标谱
omegal=omegaup/N:omegaup/N:  omegaup;    %△ω，2△ω，3△ω......N△ω
Sw1=200*U^2.*z/Uz./(1+50.*omegal.*z/(2*pi*Uz)).^(5/3);
for j=1:m                             %对m条风速循环
   rand('state',0) ;                 %模拟随机风序列
   Thet=2*pi*rand(j,N) ;                 %生成随机相位角
   for l=1:N                            %生成频率序列
   omega(l)=(l-1)*domega+j/m*domega;    %ωlk双索引频率
   end
end
Sw=200*U^2.*z/Uz./(1+50.*omega.*z/(2*pi*Uz)).^(5/3);
%计算谱数据库矩阵，功率谱计算，Kaimal谱