脉动风时程matlab程序
根据风的记录，脉动风可作为高斯平稳过程来考
虑。
观察
个具有零均值的平稳高斯过程，其谱密度函数矩阵为：
（9）
将
进行Cholesky分解，得有效方法。
（10）
其中，
（11）
为
的共轭转置。
根据文献[8]，对于功率谱密度函数矩阵为
的多维随机过程向量，模拟风速具有如下形式：
（12）
其中，风谱在频率范围内划分成
个相同部分，
为频率增量，
为上述下三角矩阵的模，
为两个不同作用点之间的相位角，
为介于
和
之间均匀分布的随机数，
是频域的递增变量。
文中模拟开孔处的来流风，因而只作单点模拟。
即式（4）可简化为：
（13）
本文采用Davenport水平脉动风速谱：
（14）
式中，
脉动风速功率谱；
脉动风频率(Hz)；
地面粗糙度系数；
标准高度为10m处的风速(m/s)。
Matlab程序:
N=10;
d=0.001;
n=d:d:N;%%频率区间（0.01～10）
v10=16;
k=0.005;
x=1200*n/v10;
s1=4*k*v10^2*x.^2./n./(1+x.^2).^(4/3);%%Davenport谱subplot(2,2,1)
loglog(n,s1)%%画谱图
axis([-100 15 -100 1000])
xlabel('freq');
ylabel('S');
for i=1:1:N/d
H(i)=chol(s1(i));%%Cholesky分解
end
thta=2*pi*rand(N/d,1000);%%介于0和2pi之间均匀分布的随机数t=1:1:1000;%%时间区间（0.1～100s）
for j=1:1:1000
a=abs(H);
b=cos((n*j/10)+thta(:,j)');
c=sum(a.*b);
v(j)=(2*d).^(1/2)*c;%%风荷载模拟
end
subplot(2,2,2)
plot(t/10,v)%%显示风荷载
xlabel('t(s)');
ylabel('v(t)');
Y=fft(v);%%对数值解作傅立叶变换
Y(1)=[];%%去掉零频量
m=length(Y)/2;%%计算频率个数；
power=abs(Y(1:m)).^2/(length(Y).^2);%%计算功率谱
freq=10*(1:m)/length(Y);%%计算频率，因为步长为0.1，而不是1，故乘以10
subplot(2,2,3)
loglog(freq,power,'r',n,s1,'b')%%比较
axis([-100 15 -100 1000])
xlabel('freq');
ylabel('S');
对源程序的修改：
z=xcorr(v);
Y=fft(z);%%对数值解作傅立叶变换
Y(1)=[];%%去掉零频量
m=length(Y)/2;%%计算频率个数；
power=abs(Y(1:m)).^2/(length(Y).^2);%%计算功率谱
freq=10*(1:m)/length(Y);%%计算频率，因为步长为0.1，而不是1，故乘以10
subplot(2,2,3)
loglog(freq,power,'r',n,s1,'b')%%比较
axis([-100 15 -100 1000])
xlabel('freq');
ylabel('S');
楼主的修改使模拟得到的功率谱与源谱的数量级对上了，但是吻合不是太好。
但是好像这样做是不对的。
求信号x(t)的功率谱有两种方法，一是对X(t)做傅立叶变换，再平方
S=abs(fft(x))^2
一是先对X(t)求相关系数，再进行傅立叶变换:
S=fft(xcorr(X))
楼主的方法好像是这两个方法的混合。
欢迎大家拍砖^_^
N=1024;
d=0.01;
n=d:d:N/100;
v10=16;
k=0.03;
x=1200*n/v10;
s1=4*k*v10^2*x.^2./n./(1+x.^2).^(4/3);%%Davenport谱for i=1:1:N
H(i)=chol(s1(i));%%Cholesky分解
end
rand('state',0);
thta=2*pi*rand(1,N);%%介于0和2pi之间均匀分布的随机数i=sqrt(-1);
B=H.*exp(i*thta);
G=fft(B,2*N);
for p=1:1:2*N
v(p)=2*sqrt(d)*real(G(p)*exp(i*2*pi*d*p*0.1/3));
end
[power,freq]=psd(v,1024*2,10,boxcar(1024),0,'mean'); power = power * 2 *0.1;%规一化修正
loglog(freq,power,'r',n,s1,'b')。
(完整版)脉动风时程matlab程序
根据风的记录，脉动风可作为高斯平稳过程来考虑。
观察n 个具有零均值的平稳高斯过程，其谱密度函数矩阵为：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)(...)()(............)(...)()()(...)()()(212222111211ωωωωωωωωωωnn n n n n s s s s s s s s s S （9）将)(ωS 进行Cholesky 分解，得有效方法。
T H H S )()()(*ωωω⋅= （10）其中，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)(...)()(............0...)()(0...0)()(21222111ωωωωωωωnn n n H H H H H H H （11） T H )(*ω为)(ωH 的共轭转置。
根据文献[8]，对于功率谱密度函数矩阵为)(ωS 的多维随机过程向量，模拟风速具有如下形式：[]∑∑==++⋅∆⋅=j m N l ml l jm l l jm j t H t v 11)(cos 2)()(θωψωωω n j ...,3,2,1= （12）其中，风谱在频率范围内划分成N 个相同部分，N ωω=∆为频率增量，)(l jm H ω为上述下三角矩阵的模，)(l jm ωψ为两个不同作用点之间的相位角，ml θ为介于0和π2之间均匀分布的随机数，ωω∆⋅=l l 是频域的递增变量。
文中模拟开孔处的来流风，因而只作单点模拟。
即式（4）可简化为：[]∑=+⋅∆⋅=Nl l l l t H t v 1cos 2)()(θωωω （13）本文采用Davenport 水平脉动风速谱：3/422210)1(4)(x n kx v n S v += （14） 式中，－－)(n S v 脉动风速功率谱；－－n 脉动风频率(Hz)；－－k 地面粗糙度系数；;120010v n x－－10v 标准高度为10m 处的风速(m/s)。
Matlab 程序:N=10;d=0.001;n=d:d:N;%%频率区间（0.01～10）v10=16;k=0.005;x=1200*n/v10;s1=4*k*v10^2*x.^2./n./(1+x.^2).^(4/3);%%Davenport 谱subplot(2,2,1)loglog(n,s1)%%画谱图axis([-100 15 -100 1000])xlabel('freq');ylabel('S');for i=1:1:N/dH(i)=chol(s1(i));%%Cholesky 分解endthta=2*pi*rand(N/d,1000);%%介于0和2pi 之间均匀分布的随机数t=1:1:1000;%%时间区间（0.1～100s ）for j=1:1:1000a=abs(H);b=cos((n*j/10)+thta(:,j)');c=sum(a.*b);v(j)=(2*d).^(1/2)*c;%%风荷载模拟endsubplot(2,2,2)plot(t/10,v)%%显示风荷载xlabel('t(s)');ylabel('v(t)');Y=fft(v);%%对数值解作傅立叶变换Y(1)=[];%%去掉零频量m=length(Y)/2;%%计算频率个数；power=abs(Y(1:m)).^2/(length(Y).^2);%%计算功率谱freq=10*(1:m)/length(Y);%%计算频率，因为步长为0.1，而不是1，故乘以10subplot(2,2,3)loglog(freq,power,'r',n,s1,'b')%%比较axis([-100 15 -100 1000])xlabel('freq');ylabel('S');1010100102freq S 050100-20-1001020t(s)v (t )10-2100100freq S对源程序的修改：z=xcorr(v);Y=fft(z);%%对数值解作傅立叶变换Y(1)=[];%%去掉零频量m=length(Y)/2;%%计算频率个数；power=abs(Y(1:m)).^2/(length(Y).^2);%%计算功率谱freq=10*(1:m)/length(Y);%%计算频率，因为步长为0.1，而不是1，故乘以10subplot(2,2,3)loglog(freq,power,'r',n,s1,'b')%%比较axis([-100 15 -100 1000])xlabel('freq');ylabel('S');楼主的修改使模拟得到的功率谱与源谱的数量级对上了，但是吻合不是太好。
风时程生成程序技术说明.
风时程⽣成程序技术说明.⽬录1程序原理 (3)1.1风荷载动⼒分析⽅法简介 (3)1.2风速时程模拟的AR法 (4)1.2.1AR模型 (4)1.2.2AR模型模拟风速时程的基本过程 (5)1.3风时程⽣成程序实现 (7)1.4风时程⽣成程序特点 (9)1.5风时程⽣成程序局限性说明 (10)2参数说明 (11)2.1顺向脉动风速功率谱密度函数()S n (11)v2.2脉动风空间相⼲函数r (13)ij2.3地⾯粗糙系数k（紊流度） (14)2.4平均风速v (14)F x y z t (16)2.5风压⼒时程(,,,)w2.6数值计算的参数 (17)3操作说明 (18)3.1制作空间点信息表格（*.csv） (18)3.2导⼊表格及输⼊参数 (19)3.3计算风时程 (20)3.4显⽰计算结果 (20)3.5输出时程结果及分析代码 (21)3.6接⼒SAP2000进⾏时程分析 (21)3.7接⼒ETABS进⾏时程分析 (22)3.8SAP2000与ETABS的分析代码例⼦ (23)3.8.1ETABS分析代码 (23)3.8.2SAP02000分析代码： (24)4计算实例 (25)4.2.1风速时程结果 (29)4.2.2风振分析计算结果与按现⾏《荷载规范》得出的结果对⽐ (31)4.2.3风振分析的顶点加速度计算与按《⾼钢规》⼿算结果对⽐ (32)5关于风振时程分析的若⼲建议 (34)5.1分析参数设置 (34)5.2输出结果处理 (34)6参考⽂献 (36)程序原理风荷载动⼒分析⽅法简介风荷载是作⽤在结构上的重要动⼒荷载之⼀，尤其对于⾼层、⾼耸及⼤跨结构来说，设计中必须考虑风荷载的作⽤。
计算⾼层、⼤跨、悬索桥以及塔架结构的动⼒风振相应的⼀个有效⽅法是Monte Carlo法。
即根据某些既定的统计参数产⽣⼀系列的时程样本，再对每个样本函数进⾏线性或⾮线性的结构分析。
通过对结构不同单元在样本函数下的时程响应的统计分析，计算整个结构是否安全。
基于线性滤波法的脉动风速模拟
收稿 日期 ： ０６一 ２ １ ２０ １ 一 ６
线 归滤波器法。 性回 研究表明‘ ， 〔 对大型工程 ］ ｓ ，
结构而言， 其自由度是非常大的。特别是大型空 间结构， 对风荷载的三维分布都比较敏感， 必须精 确模拟各点的风谱。Ｃ ｗ 法与 ｗ Ｗ ＡＳ Ａ Ａ法计算 量巨大， 所产生的风速过程不能考虑时间相关性； ｏｒ Ｓａ 提出的线性回归滤波器法很容易求出模型 ｌ 参数， 但模型精度受风谱变化的影响， 风谱的 差异 越大， 精度越低； ａ ｉ ｗ ａ Ｉｔ 提出的线性回归滤波器 ｎ 法具有较好的普适性， 但模型参数一般采用迭代、 递推的方法求解， 容易产生并累积误差， 导致模型
１ 引
言
目 国内外对风速时程的模拟方法主要是 前， Ｃ Ｗ （ ｏｓｎ Ａ ｐｔ ｅｗ ｖ ｕｅ ｏｔｎ Ａ Ｓ Ｃｎａｔ ｍｆ ｄ ａ ｔ ｉ ｕ ｅｓｐ印 ｓｏ ） ｉ ｉ 法、 Ａ Ａ Ｗ ｅ ｗｈｗ ｉｔ ｍｌ ｄ） Ｗ Ｗ （ ａ ｓ ｉ ｅｈ ｄＡ ｐｔ ｅ法及 ｖ ｔ ｇｅ ｉ ｕ
｝（） （） ：ｔ｝＝［ 」ｕｔ｝ Ｃ ｛（） １ 式中，ｕｔ｝ ｛ ）为互不相关的高斯过程；创为互相 （ 仁 关矩阵， 可以由后面的公式求得。 因而， 在某一时刻ｔ 、 高度： 风速 侧： 可 处， ， ） ｔ 以看作高度 ： 处的平均风速 ：：与脉动风速 ｕ （） （， 之和。 ２ｔ ）
的精度不够。
万方数据
Ｓ ｃ ｒ Ｅ ｇ ｅｒ Ｖ ｌ ３ Ｎ ． ｔ ｔ ａ ｎｉ ｅ ｏ ２ ，ｏ４ ｕ ｒ ｕｌ ｎ ｓ ．
Ｅｒ ｑａｅ ｎ ｎ ｅ ｓ ｎｅ ａｈｕｋ ａｄｗｉ ｔ ｄＲ ｓｔ ｃ ｉａ
２ 风的基本特性
风对结构的作用可以看成由平均风作用和脉 动风作用两部分组成， 其中平均风在一定的时间 间隔内， 风的大小和方向不随时间变化； 经过实测 风时程记录可知， 平均风剖面沿结构高度往往按 指数或者对数规律变化； 而脉动风荷载是随机荷 载， 是风力中的动力部分， 它使结构产生随机振
基于Matlab的大型兆瓦级风电机脉动风速时程数值模拟
第 ４期
曹玉 生等
基于 Ｍ ａ ｔ ｌ ａ ｂ的大型兆瓦级风电机 脉动风 速时程数值模拟
２ ７９
作 用在 结构 上 的 自然风 可分 为顺 风 向风 力 、 横 风 向风力 和 垂直 向风 力 ， 而在 三 者 中顺 风 向风 力 又起 了决
定性作用 ， 垂直 向与横风 向风力对于高耸结构实际影响可忽略不计。对顺风向风力 ， 其可分解为周期在 １ ０ ｍ ｉ ｎ以上 的长周 期 部分 ， 即平 均 风 ； 还 有周 期在 几秒 至几 十 秒 区间 内的短周 期 部分 ， 即脉 动风 。在 进行
型风 电机得 到 了大规 模 的应 用 。 因此 ， 在 内陆 地 区普 及 大 型 兆 瓦 级 风 电 机组 是 风 能发 展 的必 然 趋 势 。 由于对大 型风 电机组 进行 实地 的数 据测 量及 周边环 境 分析工 作 的开展 是 极其 困难 的 ， 另外 ， 已有 强风 作
matlab单自由度的时程分析程序
clear;clc;% 结构模型初始参数---------------------------------------------------------- m=3e3; %质量(单位：kg)k=1e6; %刚度（(单位：N/m)）kesai=0.05; %阻尼比取0.05c=2*kesai*sqrt(k*m); %阻尼系数% 读取地震波数据------------------------------------------------------------ acc=textread('D:\处理后的smc文件\51WCW_90_chnua370295.smc_090501.a','%f','headerlines',56); PGA_Max=max(abs(acc)) %最大地面加速度绝对值% Newmark-beta法的基本参数--------------------------------------------------beta=1/6; gama=0.5; %按线性加速度法计算更接近真实结果，故取此组参数dt=0.02; %地震加速度时程波记录时间间隔b1=1/(beta*dt^2); b2=1/(beta*dt); b3=1-1/(2*beta); %计算参数b4=gama/(beta*dt); b5=gama/beta-1; b6=(1-gama/(2*beta)) *dt;ke=k+m*b1+c*b4; %等效刚度% 设定结构初始状态为零,生成向量空间存储计算值---------------------------------u=zeros(100/dt,1); v=zeros(100/dt,1); a=zeros(100/dt,1);% Newmark-beta法的主计算程序------------------------------------------------for n=2:100/dtfe=-m*acc(n)+[b1*u(n-1)+b2*v(n-1)-b3*a(n-1)]*m+[b4*u(n-1) +b5*v(n-1)-b6*a(n-1)]*c; %等效荷载u(n)=fe/ke;a(n)=b1*[u(n)-u(n-1)]-b2*v(n-1)+b3*a(n-1);v(n)=b4*[u(n)-u(n-1)]-b5*v(n-1)+b6*a(n-1);end% 绘制结构在地震作用下的位移、速度、加速度时程曲线-----------------------------subplot(3,1,1)t=(0:length(a)-1)*dt;plot(t,a) %加速度时程曲线Acc_Max=max(abs(a))title('Earthquake Response Curve Of Station 51WCW-90','fontsize',15) ylabel('Acc(cm/s^2)','fontsize',12)subplot(3,1,2)plot(t,v) %速度时程曲线Vel_Max=max(abs(v))ylabel('Vel(cm/s)','fontsize',12)subplot(3,1,3)plot(t,u) %位移时程曲线Dis_Max=max(abs(u))xlabel('Time/s','fontsize',12)ylabel('Dis/cm','fontsize',12)% End---程序结束-------------小弟初次发贴，恳请达人们帮分析一下，不胜感激！其中的循环部分是根据结构动力学书上的写的，感觉问题就出在那部分了，请高人们指点一下线性加速度法是直接数值积分法求解地震反应的方法之一，本文所采用的线性加速度法参考大崎顺彦的《地震动的谱分析入门》第二版。
Matlab 时程分析
Matlab 时程分析0 动力平衡方程及相关参数取值波浪、风载作用下的单桩动力反应计算（把结构简化为质点剪切型）wave wind []{}[]{}[]{}M u C u K u F F ∙∙∙++=+式中：u 为结构水平位移；[]M 桩-结构集中质量矩阵；[][][]P S C C C =+体系的阻尼矩阵；体系的阻尼矩阵由结构和土体的阻尼矩阵集成，其中结构的阻尼按瑞雷阻尼理论，土体阻尼由材料阻尼和辐射阻尼组成。
[][][]P S K K K =+体系的刚度矩阵；体系的刚度矩阵由结构和土体的刚度矩阵集成，土体刚度由动力P-Y 曲线对Y 求导得到。
地震、波浪、风载作用下的单桩动力反应计算1 自由场地震分析（远离桩，取单位面积土柱）[]{}[]{}[]{}[]{}f f f f f f f g M u C u K u M E u ∙∙∙∙∙++=- 土的刚度矩阵：；/f i i i k G h =土的阻尼矩阵：；离桩较远，可采用刚度比例阻尼2 桩土相互作用地震分析[]{}[]{}[]{}[]{}[]{}[]{}s s f f g wave wind M u C u K u M E u F F C u K u ∙∙∙∙∙∙++=-++++也可写为[]{}[]{}[]{}[]{}[]{}[]{}P S P S f f g wave wind M u C u C u u K u K u u M E u F F ∙∙∙∙∙∙∙++-++-=-++[]K 桩-结构集中质量矩阵；[]P C 结构的阻尼矩阵；[]S C 土体的阻尼矩阵；[]P K 结构的刚度矩阵；[]S K 土体的刚度矩阵；g u ∙∙基岩加速度。
自由场反应{}{0()}f T T f u u ∙∙=；{}{0()}fT T f u u = 1 结构动力响应数值方法威尔逊θ法，隐式积分格式，θ一般取1.4时无条件稳定。
θ=1即为常规线性加速度法。
风力机MATLAB设计程序
makedata%根据profili导出到翼型性能数据Excel表格生成翼型的结构体clear airfoil;for n=2:100 %sheet nairfoil(n-1).Re=22;%%%%找到sheet n 截面翼型雷诺数所在行数nRe%%%%try[num,txt,~] = xlsread('yxdata.xlsx',n,'A1:I5000');catchbreak;endnstr=strfind(txt,'Re = ');%找到每一个雷诺数翼型的起始记录位置k=1; %nRe的变量for i=1:length(nstr)%行数从第一行到最后一行开始判断j=nstr{i}; %将第i行的值赋给临时变量jif j %如果j存在则将行数给nReairfoil(n-1).nRe(k)=i;k=k+1;endairfoil(n-1).nRe(k)=length(num)+5;end%%%%%%%%%%%%%%%%%%找到翼型相近的雷诺数下的性能数据和name%%%%%%%%%%%%%%%%%%%k=length(airfoil(n-1).nRe);airfoil(n-1).Re(k)=0;airfoil(n-1).name=txt{airfoil(n-1).nRe(1)}(1:(nstr{airfoil(n-1).nRe(1)}-4));%将第i行的值赋给临时变量jfor i=1:length(airfoil(n-1).nRe) %行数从第一行到最后一行开始判断airfoil(n-1).Re(i)=str2double(txt{airfoil(n-1).nRe(i)}((nstr{airfoil(n-1).nRe(i) }+5):length(txt{airfoil(n-1).nRe(i)}))); %将第i行的值赋给临时变量jend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%try[wnum,~,~] = xlsread('yxdata.xlsx',n,'H1:I500');lth=length(wnum);airfoil(n-1).x(1:lth,1)=wnum(1:lth,1);airfoil(n-1).y(1:lth,1)=wnum(1:lth,2);catchend%%%%读入截面翼型拟合各雷诺数性能曲线和其它数据%%%%for i=1:(length(airfoil(n-1).nRe)-1) %Retemp=(airfoil(n-1).nRe(i):(airfoil(n-1).nRe(i+1)-5));lth=length(temp);airfoil(n-1).Alf(1:lth,i)=num(temp,1);airfoil(n-1).Cl(1:lth,i)=num(temp,2);airfoil(n-1).Cd(1:lth,i)=num(temp,3);airfoil(n-1).ClCd(1:lth,i)=num(temp,4);tempn=find(airfoil(n-1).ClCd(:,i)==max(airfoil(n-1).ClCd(:,i)));airfoil(n-1).zAlf(i)=airfoil(n-1).Alf(tempn,i);airfoil(n-1).zCl(i)=airfoil(n-1).Cl(tempn,i);airfoil(n-1).zCd(i)=airfoil(n-1).Cd(tempn,i);[airfoil(n-1).xCl(:,i) airfoil(n-1).SxCl(:,i) ]= polyfit(airfoil(n-1).Alf(:,i),airfoil(n-1).Cl(:,i),6);[airfoil(n-1).xCd(:,i) airfoil(n-1).SxCd(:,i)] = polyfit(airfoil(n-1).Alf(:,i),airfoil(n-1).Cd(:,i),6);endendsave airfoildataqdclc;clear;filename='name';load(filename)%load xcload airfoilData airfoilpi=3.141592653;qR=287.64;k=1.4;fq=0.12;u=1.698e-05;%pi 圆周率；qR气体常数；k 等商指数；fq 风切指数；u 动力粘度；Pr=1200000;Ve=8.5;Pa=85.8;T=15;B=3;DJ_eta=0.95;CD_eta=0.95;%Pr 额定功率；Vr 额定风速；Pa 风场平均压强；T 平均气温；B 叶片数%DJ_eta 电机效率；CD_eta传动效率；%tempV=70;Cp=0.43;n=30;namR=9;BL1=0.15;BL2=0.05;%Cp 风能利用系数；n 等分段数；namR ;叶尖速比;BL1 叶根园比例；BL1 轮毂园比例；min_n=900;max_n=1950;e_n=1620;%发电机的转速范围iname=1; %%%%%%%%%%%%%%%%%%%%%%%%%开始迭代计算轮毂高度%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Hhub=95;temp=0;while abs(Hhub-temp)>2Vr=Ve*(Hhub/10)^fq;rou=Pa*1000/((273+T)*qR);%Vr 设计风速；rou 空气密度D=(8*Pr/(Cp*DJ_eta*CD_eta*rou*Vr^3*pi))^0.5;D1=floor(D);%取比圆整风轮直径向上取对Cp和功率的大小又决定性作用。
(完整版)脉动风时程matlab程序
根据风的记录，脉动风可作为高斯平稳过程来考虑。
观察n 个具有零均值的平稳高斯过程，其谱密度函数矩阵为：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)(...)()(............)(...)()()(...)()()(212222111211ωωωωωωωωωωnn n n n n s s s s s s s s s S （9）将)(ωS 进行Cholesky 分解，得有效方法。
T H H S )()()(*ωωω⋅= （10）其中，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=)(...)()(............0...)()(0...0)()(21222111ωωωωωωωnn n n H H H H H H H （11） T H )(*ω为)(ωH 的共轭转置。
根据文献[8]，对于功率谱密度函数矩阵为)(ωS 的多维随机过程向量，模拟风速具有如下形式：[]∑∑==++⋅∆⋅=j m N l ml l jm l l jm j t H t v 11)(cos 2)()(θωψωωω n j ...,3,2,1= （12）其中，风谱在频率范围内划分成N 个相同部分，N ωω=∆为频率增量，)(l jm H ω为上述下三角矩阵的模，)(l jm ωψ为两个不同作用点之间的相位角，ml θ为介于0和π2之间均匀分布的随机数，ωω∆⋅=l l 是频域的递增变量。
文中模拟开孔处的来流风，因而只作单点模拟。
即式（4）可简化为：[]∑=+⋅∆⋅=Nl l l l t H t v 1cos 2)()(θωωω （13）本文采用Davenport 水平脉动风速谱：3/422210)1(4)(x n kx v n S v += （14） 式中，－－)(n S v 脉动风速功率谱；－－n 脉动风频率(Hz)；－－k 地面粗糙度系数；;120010v n x－－10v 标准高度为10m 处的风速(m/s)。
Matlab 程序:N=10;d=0.001;n=d:d:N;%%频率区间（0.01～10）v10=16;k=0.005;x=1200*n/v10;s1=4*k*v10^2*x.^2./n./(1+x.^2).^(4/3);%%Davenport 谱subplot(2,2,1)loglog(n,s1)%%画谱图axis([-100 15 -100 1000])xlabel('freq');ylabel('S');for i=1:1:N/dH(i)=chol(s1(i));%%Cholesky 分解endthta=2*pi*rand(N/d,1000);%%介于0和2pi 之间均匀分布的随机数t=1:1:1000;%%时间区间（0.1～100s ）for j=1:1:1000a=abs(H);b=cos((n*j/10)+thta(:,j)');c=sum(a.*b);v(j)=(2*d).^(1/2)*c;%%风荷载模拟endsubplot(2,2,2)plot(t/10,v)%%显示风荷载xlabel('t(s)');ylabel('v(t)');Y=fft(v);%%对数值解作傅立叶变换Y(1)=[];%%去掉零频量m=length(Y)/2;%%计算频率个数；power=abs(Y(1:m)).^2/(length(Y).^2);%%计算功率谱freq=10*(1:m)/length(Y);%%计算频率，因为步长为0.1，而不是1，故乘以10subplot(2,2,3)loglog(freq,power,'r',n,s1,'b')%%比较axis([-100 15 -100 1000])xlabel('freq');ylabel('S');1010100102freq S 050100-20-1001020t(s)v (t )10-2100100freq S对源程序的修改：z=xcorr(v);Y=fft(z);%%对数值解作傅立叶变换Y(1)=[];%%去掉零频量m=length(Y)/2;%%计算频率个数；power=abs(Y(1:m)).^2/(length(Y).^2);%%计算功率谱freq=10*(1:m)/length(Y);%%计算频率，因为步长为0.1，而不是1，故乘以10subplot(2,2,3)loglog(freq,power,'r',n,s1,'b')%%比较axis([-100 15 -100 1000])xlabel('freq');ylabel('S');楼主的修改使模拟得到的功率谱与源谱的数量级对上了，但是吻合不是太好。
基于MATLAB-的脉搏信号处理软件系统
基于MATLAB 的脉搏信号处理软件系统摘要: 本文根据在实验室里测得的脉搏数据，基于MATLBA设计一个脉搏信号的GUI处理界面，并利用MATLAB强大数字信号处理功能复原脉搏波形，并对波形的特征信息进行提取及存储。
原始信号进行了去除基线漂移、通过巴特沃斯带通滤波器以及二阶切比雪夫滤波器去除50HZ工频干扰，并且能计算实时的脉率并更新，显示脉率变化趋势曲线，进行频谱分析和输出文档。
此软件有两个GUI界面，第一个为密码登陆界面，第二个为脉搏信号处理系统GUI界面。
第二个GUI界面主要分为五大模块:1.打开与退出模块包括打开数据和退出系统;2.信号回放模块包括对原信号和滤波信号的回放、暂停回放、继续回放、关闭窗口;3.信号放大与缩小模块包括对信号的X轴和Y轴的放大、缩小处理;信号快进退模块包括对信号的快进、慢进、快退、慢退处理;4.脉率实时处理模块包括输出脉率曲线、暂停回放、输出脉搏信息、脉搏频谱分析、清除波形、输出文档;5.脉率信号输出模块包括输出实时的脉率更新、以及脉搏数据的信息，诸如脉搏采样频率、采样时间、最大脉率值、最小脉率等。
关键词：脉搏;脉率;Matlab ;GUI ;1 引言人体内部各个生理系统之间(如循环系统、呼吸系统等)是相互耦合的。
反映人身体健康状态相对最重要、最全面的是心脏血液循环系统,因此通过采集脉搏波进而分析心脏循环系统功能,能从一个方面较全面反映人体的健康情况。
从脉搏波中提取人体的生理病理信息作为临床诊断和治疗的依据，历来都受到中外医学界的重视。
几乎世界上所有的民族都用过“摸脉”作为诊断疾病的手段。
脉搏波所呈现出的形态(波形)、强度(波幅)、速率(波速)和节律(周期)等方面的综合信息，在很大程度上反映出人体心血管系统中许多生理病理的血流特征,因此对脉搏波采集和处理具有很高的医学价值和应用前景。
目前脉搏信息的研究已经应用于以下几个方面:(1)中医脉象信息的检测与识别;(2)血压的临床检测;(3)心率稳定性的一种简便估计方法;(4)心输出量的一种测量方法;(5)血管功能的一种早期、无创检测方法。
风电机组塔架的脉动风速时程模拟
风电机组塔架的脉动风速时程模拟叶赟;宫兆宇【摘要】In this paper, an autoregressive model(AR model) is used to simulate wind speed time series. e spectrum of simulated wind speed time series is found in agreement with the target spectrum, Davenport wind speed spectrum. Samples of the uctuating wind load on the nodes of a structure are obtained. Using the WAWS, the article builds an AR model to calculate the model order and edit a simulation program. rough the analysis on some wind turbines tower, the feasibility and e ciency of this simulation model is veri ed.% 本文简述了谐波合成法中的自回归模型(AR)模拟出给定风速功率谱的风速时程序列，并验证其与目标谱(Davenport谱)的一致性，从而得到作用在各节点的脉动风荷载时程样本的方法。
本文采用谐波合成法，建立了脉动风速时程的 AR 模型，编辑出脉动风速时程模拟程序，并对某风电机组塔架进行脉动时程分析，验证了该脉动风速时程模拟的可行性与有效性。
【期刊名称】《风能》【年(卷),期】2013(000)001【总页数】6页(P72-77)【关键词】脉动风;数值模拟;自回归模型;风电机组塔架【作者】叶赟;宫兆宇【作者单位】内蒙古科技大学建筑与土木工程学院，包头 014010;内蒙古科技大学建筑与土木工程学院，包头 014010【正文语种】中文【中图分类】TM614高耸结构风荷载是结构设计时必须要考虑的一类重要的随机荷载，风振响应成为控制结构设计的重要因素。