%諧波疊加法模擬風(fēng)速時(shí)程（Kaimal譜）
clc
clear

%風(fēng)速時(shí)程參數(shù)設(shè)定
m=10;?????%模擬風(fēng)速個(gè)數(shù)
N=2^3;????%頻率采樣點(diǎn)數(shù)，采樣頻率是指計(jì)算機(jī)每秒鐘采集多少個(gè)風(fēng)速樣本
dt=0.5;???%時(shí)間間隔
omegaup=2*pi;?%上限頻率，pi為圓周率，應(yīng)該不是“上限頻率”
%設(shè)定鳳速譜參數(shù)
L=1000;???%斜拉橋跨度
z=50;?????%風(fēng)速測(cè)點(diǎn)離地面高度
z0=0.03;??%地面粗糙度
Uz=40;????%50m處的平均風(fēng)速
delta=100;?%模擬風(fēng)速測(cè)點(diǎn)間的距離
lambda=10;?%空間相關(guān)函數(shù)中的衰減系數(shù)C，需要單獨(dú)定義！
K=0.4;?????%Kaman常數(shù)0.4
M=2*N;?????%nfft傅里葉變換長(zhǎng)度，取采樣個(gè)數(shù)的2倍

%形成風(fēng)速時(shí)程矩陣
v=zeros（mM*m）;%創(chuàng)建m行，M*m列的時(shí)程矩陣;
t=dt*（0:1:（M*m-1））;%創(chuàng)建時(shí)程橫坐標(biāo)時(shí)間點(diǎn)
domega=omegaup/N;%頻率間距
D=zeros（mmN）;%三向矩陣的意義是？
U=K*Uz/log（z/z0）;%測(cè)點(diǎn)位置的摩阻速度，也可以轉(zhuǎn)換為對(duì)述律表示
disp（U）
%形成目標(biāo)譜
omega1=omegaup/N:domega:omegaup;%形成頻率列表，初始化頻率列表，初始頻率為omegaup/N
Sw1=200*U^2.*z/Uz./（1+50.*omega1.*z./（2*pi*Uz））.^（5/3）;%Kaimal譜密度表達(dá)式
for?j=1:m???%對(duì)模擬點(diǎn)風(fēng)速的循環(huán)
????%rand（‘state‘0）;???%生成均勻分布隨機(jī)矩陣?%rand（‘state‘s）
???????????????????????%表示隨機(jī)產(chǎn)生數(shù)的狀態(tài)state，一般情況下不用指定狀態(tài)。
???????????????????????%但是有的書籍作者為了讓讀者能夠看到和他書本的同樣結(jié)果，
???????????????????????%采用了設(shè)置state，rand（‘state‘0）作用在于如果指定狀態(tài)，產(chǎn)生隨機(jī)結(jié)果就相同了。
????thet=2*pi*rand（jN）;%生成隨機(jī)相位
????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）

%計(jì)算譜數(shù)據(jù)庫(kù)矩陣，功率譜計(jì)算，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））;