%諧波疊加法模擬風速時程（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））;