线性叠加法

（大家反映程序没法运行，原来是之前没有放子程序的缘故，这里统一放上~）
海浪可看做一系列不同周期不同初相位的线性波叠加而成的：
η(t)=∑i=1Maicos⁡(kix−ωit+ϵi)\eta(t)=\sum\limits_{i=1}^{M}a_i\cos(k_ix-\omega_it+\epsilon_i)η(t)=i=1∑Maicos(kix−ωit+ϵi),
aia_iai为第iii个组成波的振幅，ki和ωik_i和\omega_iki和ωi为第iii个组成波的波数和圆频率。ϵi\epsilon_iϵi为(0,2π)(0,2\pi)(0,2π)之间的随机数，代表随机相位。假设靶谱的能量大多分布在区间[ωLωH][\omega_L\quad\omega_H][ωLωH],其他部分可忽略不计。将该区间平分为M个子区间，其间距为Δωi=ωi−ωi−1\Delta\omega_i=\omega_i-\omega_{i-1}Δωi=ωi−ωi−1，取
ωi^=(ωi−1+ωi)/2,ai=2Sηη(ωi^)Δωi\hat{\omega_i}=(\omega_{i-1}+\omega_i)/2,a_i=\sqrt{2S_{\eta\eta}(\hat{\omega_i})\Delta\omega_i}ωi^=(ωi−1+ωi)/2,ai=2Sηη(ωi^)Δωi,
则将代表M个区间波能的M个线性波叠加起来，可得到海浪波面：
η(t)=∑i=1M2Sηη(ωi^)Δωicos⁡(ωi~t+ϵi)\eta(t)=\sum\limits_{i=1}^{M}\sqrt{2S_{\eta\eta}(\hat{\omega_i})\Delta\omega_i}\cos(\tilde{\omega_i}t+\epsilon_i)η(t)=i=1∑M2Sηη(ωi^)Δωicos(ωi~t+ϵi),
式中ωi~\tilde{\omega_i}ωi~为第i个组成波的代表频率。

程序

主程序

% Randam wave simulation
% Designed by: JN-Cui
% Modified on 12/09/2019
%% DEFINITIONS
% alpha - energy scale factor; gama - spectral peak elevation factor;
% omega_m - spectral peak circular frequency; f_m - spectral peak frequency;
% U - wind speed at 10 m above sea surface; H_s - significant wave height;
% g - gravity acceleration;
%% FOR AVERAGE JONSWAP SPECTRAL
% gama=3.3; k=83.7; sigma_a=0.07; sigma_b=0.09;
% alpha=0.076*(X_bar)^(-0.22);
% X_bar=10^(-1)~20^(5); omega_m=22(g/U)*(X_bar)^(-0.33);
% f_m=3.5(g/U)(X_bar)^(-0.33);
%% FUNCTIONfunction [W_s,t,x,Omega,S]=random_wave_v1_0(H_s,T_s,dm,dt,T,dx,X)
[S,Omega,omega_p,T_m]=Improved_Jonswap_spectral(H_s,T_s,dm);
d=2000;
g=9.8;
N=200;
M=T/dt;
e=rand(1,N)*2*pi;
w=0:max(Omega)/(N-1):max(Omega);
ww(N)=max(w);
LT1=length(w);
LT2=M;
LT3=M;
for i=1:N-1ww(i)=unifrnd(w(i),w(i+1));
end
ww(N)=w(N);
Wait=waitbar(0,'程序计算中，请稍后', 'CreateCancelBtn','setappdata(gcbf,''canceling'',1)');
setappdata(Wait,'canceling',0);
SS_o=zeros(1,length(w));
f=zeros(1,length(w));
T_w=zeros(1,length(w));
K=zeros(1,length(w));
k=zeros(1,length(w));
a=zeros(1,length(w));
L=zeros(1,length(w));
% for i=1:length(w)
for i=1:Nwaitbar(i/(LT1+LT2+LT3));if getappdata(Wait,'canceling')breakendSS_o(i)=S_Improved_Jonswap_spectral(ww(i),H_s,T_s,dm);f(i)=ww(i)/2/pi;T_w(i)=1/f(i);K(i)=g*T_w(i)^2/(2*pi);fun=@(L1) L1/tanh(2*pi/L1*d)-K(i);L(i)=Division(fun,0.0001,0,K(i));k(i)=(2*pi/L(i));a(i)=sqrt(2*SS_o(i)*(w(2)-w(1)));
end
% Time steps
t=zeros(1,M);
for i=1:Mt(i)=(i-1)*dt;
end
eta=zeros(N,1);
NN=floor(X/dx)+1;
x=0:dx:X;
Eta=zeros(M,NN-1);
for j=1:Mwaitbar((j+LT1)/(LT1+LT2+LT3));if getappdata(Wait,'canceling')breakendfor ii=1:NNEta(j,ii)=sum(a.*cos(k.*x(ii)-ww.*t(j)+e));end
end
%%
% Wave Surface
W_s=Eta;
delete(Wait);
end

子程序1

%% DEFINITIONS
% alpha - energy scale factor; gama - spectral peak elevation factor;
% omega_m - spectral peak circular frequency; f_m - spectral peak frequency;
% U - wind speed at 10 m above sea surface; H_s - significant wave height;
% g - gravity acceleration;
%% FOR AVERAGE JONSWAP SPECTRAL
% gama=3.3; k=83.7; sigma_a=0.07; sigma_b=0.09;
% alpha=0.076*(X_bar)^(-0.22);
% X_bar=10^(-1)~20^(5); omega_m=22(g/U)*(X_bar)^(-0.33);
% f_m=3.5(g/U)(X_bar)^(-0.33);
%% IMPUT PARAMETERS
% H_s - significant wave height; T_s -  wave period at 1/3 wave height
% dm - calculation interval of omega
%% FUNCTION
function [S,Omega,omega_p,T_p]=Improved_Jonswap_spectral(H_s,T_p,dm)%the input of dm??
gama=3.3; sigma_a=0.07;sigma_b=0.09;
beta_j=0.06238/(0.23+0.0336*gama-0.185*(1.9+gama)^(-1))*(1.094-0.01915*log(gama));
T_s=T_p*(1-0.132*(gama+0.2)^(-0.559));
f_p=1/T_p;
omega_p=f_p*2*pi;
i=1;
df=dm/2/pi;
S_o=zeros(1,length(0:dm:1/T_s*2*pi*4));
Omega1=zeros(1,length(0:dm:1/T_s*2*pi*4));
for omega=0:dm:1/T_s*2*pi*4 %?显著波高对应频率的四倍  包含绝大部分谱的能量if omega<omega_p sigma=sigma_a;S_o(i)=beta_j*H_s^2*T_p^(-4)*(omega/2/pi)^(-5)*exp(-5/4*((omega_p/omega))^(4))...*gama^(exp(-((omega)/omega_p-1)^2/(2*sigma^2)))/(2*pi);elsesigma=sigma_b;S_o(i)=beta_j*H_s^2*T_p^(-4)*(omega/2/pi)^(-5)*exp(-5/4*((omega_p/omega))^(4))...*gama^(exp(-((omega)/omega_p-1)^2/(2*sigma^2)))/(2*pi);endOmega1(i)=omega; i=i+1;
end
Omega=Omega1(2:end);
S=S_o(2:end);
end

子程序2

function [S_omega]=S_Improved_Jonswap_spectral(omega,H_s,T_s,dm)
gama=3.3; sigma_a=0.07;sigma_b=0.09;
beta_j=0.06238/(0.23+0.0336*gama-0.185*(1.9+gama)^(-1))*(1.094-0.01915*log(gama));
T_p=T_s/(1-0.132*(gama+0.2)^(-0.559));
% T_p=T_bar/(1-0.532*(gama+2.5)^(-0.569));
f_p=1/T_p;
omega_p=f_p*2*pi;
i=1;
df=dm/2/pi;
S_o=zeros(1,length(0:dm:1/T_s*2*pi*4));
Omega1=zeros(1,length(0:dm:1/T_s*2*pi*4));if omega<omega_psigma=sigma_a;S_omega=beta_j*H_s^2*T_p^(-4)*(omega/2/pi)^(-5)*exp(-5/4*((omega_p/omega))^(4))...*gama^(exp(-((omega)/omega_p-1)^2/(2*sigma^2)))/(2*pi);elsesigma=sigma_b;S_omega=beta_j*H_s^2*T_p^(-4)*(omega/2/pi)^(-5)*exp(-5/4*((omega_p/omega))^(4))...*gama^(exp(-((omega)/omega_p-1)^2/(2*sigma^2)))/(2*pi);end
end

子程序3

function output=Division(fun,x,a,b)
p=-1;
while (fun(a)*fun(b) <=0) && (abs(a-b)>x) c=(a+b)/2; if fun(c)*fun(b)<=0 a=c; p=p+1; elsep=p+1;b=c;end
end
output=(a+b)/2;
end

等分频率法模拟随机波列（线性波叠加原理）相关推荐

matlab谐波合成法模拟风速时程,基于谐波合成法的输电塔线体系风致响应分析
0引言塔线的耦合作用使动力风荷载作用下大跨.高柔的输电塔线体系表现出很强的非线性特征,导线与输电塔之间的风动力耦联作用不容忽略[1-2].动力风致响应分析的正确性和精确度将直接影响到输电塔线结构的安全 ...
《算法设计编程实验：大学程序设计课程与竞赛训练教材》——2.3　构造法模拟的实验范例...
2.3 构造法模拟的实验范例构造法模拟需要完整.精确地构造出反映问题本质的数学模型,根据该模型设计状态变化的参数,计算模拟结果.由于数学模型建立了客观事物间准确的运算关系,因此其效率一般比较高. 构 ...
matlab小波脊线,小波脊线提取，模极大值法。运行的结果不太对，代码有些地方我也没完全看懂...
本帖最后由 1393107100 于 2019-5-1 11:10 编辑 clear,clc close all %%%%%% 小波变换 %%%%%%%%%%%%%%%% fs=1024; t=1/f ...
matlab中拟合函数中的gian值,如何在Matlab中优化基本周期图法对随机信号进行的功率谱估计...
首都师范大学学报(自然科学版)第27卷第5期2006年10月 Journal of Capital N ormal University (Natural Science Edition ) V o ...
matlab实现鬼波信号压制算法代码实践--第二篇频率-波数域鬼波压制
第二篇 matlab实现频率-波数域鬼波压制方法本篇用来介绍频率-波数域鬼波压制的实现思路和压制效果算法实现思路见2.3节,除了文中代码外,需配置鬼波压制算法工具包(https://downlo ...
根据geojson文件模拟随机点
一.判断一个点是否在某一个封闭多边形的范围内可以参考另一篇博客:判断一个点是否在封闭多边形内我们选取其中一个python方式建立文件 is_point_in_polygon.py "&q ...
【数字信号处理】基本序列 ( 正弦序列 | 数字角频率 ω | 模拟角频率 Ω | 数字频率 f | 模拟频率 f0 | 采样频率 Fs | 采样周期 T )
文章目录一.正弦序列 ( 数字信号 ) 二.模拟角频率与数字角频率关系三.模拟信号四.数字角频率 ω 与模拟角频率 Ω 与模拟频率 f 的关系五.数字频率 f 与模拟频率 f0 的 ...
信号处理中数字频率和模拟频率简明讲解
数字频率和模拟频率在数字信号处理的学习中,很多刚入门朋友常常为模拟频率.数字频率及其相互之间的关系所迷惑,甚至是一些已经对数字信号处理有所了解的朋友也为这个问题所困惑. 我们通常所说的频率,在没有特 ...
SQLIO 模拟随机或者顺序的方式来测试磁盘IO的性能
SQLIO 功能:磁盘IO压力测试工具,SQLIO主要是模拟随机或者顺序的方式来测试磁盘IO的性能. SQLIO Disk Subsystem Benchmark Tool工具下载地址: http:/ ...
文献记录(part68)--K- 近邻分类器鲁棒性验证:从约束放松法到随机平滑法
学习笔记,仅供参考,有错必纠关键词:监督学习 , 对抗机器学习 , 对抗鲁棒性 , 鲁棒性验证 , K- 近邻分类器 K- 近邻分类器鲁棒性验证:从约束放松法到随机平滑法摘要本文研究 K- 近邻 ...

等分频率法模拟随机波列（线性波叠加原理）

线性叠加法

程序

主程序

子程序1

子程序2

子程序3

等分频率法模拟随机波列（线性波叠加原理）相关推荐

最新文章

热门文章