风速时程模拟

%谐波叠加法模拟风速时程（Kaimal谱）
clc
clear

%风速时程参数设定
m=10;     %模拟风速个数
N=2^3;    %频率采样点数，采样频率是指计算机每秒钟采集多少个风速样本
dt=0.5;   %时间间隔
omegaup=2*pi; %上限频率，pi为圆周率，应该不是“上限频率”
%设定凤速谱参数
L=1000;   %斜拉桥跨度
z=50;     %风速测点离地面高度
z0=0.03;  %地面粗糙度
Uz=40;    %50m处的平均风速
delta=100; %模拟风速测点间的距离
lambda=10; %空间相关函数中的衰减系数C，需要单独定义！
K=0.4;     %Kaman常数0.4
M=2*N;     %nfft傅里叶变换长度，取采样个数的2倍

%形成风速时程矩阵
v=zeros（mM*m）;%创建m行，M*m列的时程矩阵;
t=dt*（0:1:（M*m-1））;%创建时程横坐标时间点
domega=omegaup/N;%频率间距
D=zeros（mmN）;%三向矩阵的意义是？
U=K*Uz/log（z/z0）;%测点位置的摩阻速度，也可以转换为对述律表示
disp（U）
%形成目标谱
omega1=omegaup/N:domega:omegaup;%形成频率列表，初始化频率列表，初始频率为omegaup/N
Sw1=200*U^2.*z/Uz./（1+50.*omega1.*z./（2*pi*Uz））.^（5/3）;%Kaimal谱密度表达式
for j=1:m   %对模拟点风速的循环
    %rand（‘state‘0）;   %生成均匀分布随机矩阵 %rand（‘state‘s）
                       %表示随机产生数的状态state，一般情况下不用指定状态。
                       %但是有的书籍作者为了让读者能够看到和他书本的同样结果，
                       %采用了设置state，rand（‘state‘0）作用在于如果指定状态，产生随机结果就相同了。
    thet=2*pi*rand（jN）;%生成随机相位
    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）;
  % disp（Sw）

%计算谱数据库矩阵，功率谱计算，Kaimal谱
    for j1=1:m
        for l=1:m
            for k=1:N
                Coh（j1lk）=（exp（-lambda*omega（k）*delta/（2*pi*Uz）））^（abs（j1-l））;