完整程序：
clc
clear
syms w we;
w=0.4:0.05:1.6;mu=90;v=6;%kb=1;kt=1;%航速6m/s，航向90度，即横浪，cos(90)=0
T=3;B=10;Sw=0.785;%船宽10米，吃水3米,水线面系数假设为0.785
we=w+w.^2.*v/9.8;
for i=1:24
delta_we(i)=we(i+1)-we(i);
end
delta_we(25)=delta_we(24);
z=0.194375;
s_kesi=0.778./(w.^5).*exp(-z./(w.^4)); %含风速=15m/s即三一波高=4m在内了 %2.47./(w.^6).*exp(-2.*9.8.^2./24.5.^2./(w.^2));
s_alpha=w.^4.*s_kesi./(9.8^2);%波倾角谱
%以下进行波倾角修正
%按照船宽和吃水定kb,kt,进而求s_alpha_v;
kb=exp(-w.^2.*T./(6.24*pi));
kt=1-sqrt(Sw).*B.^2.*w.^4./(1.56^2*16*pi^4);
s_kesi_v=0.778./(w.^5).*exp(-0.194375./(w.^4))./(1+2.*w./9.8.*v.*cos(mu*pi/180));
s_alpha_v=w.^4.*s_kesi_v./(9.8^2).*kb.^2.*kt.^2;
%figure;hold on;set(gca,'XLim',[0.2,2]);
%plot(w,s_kesi,'r',w,s_alpha,'g--');
s_alpha_ap=sqrt(2.*s_alpha_v.*0.05);%波倾角幅值
s_alpha_v_ap=sqrt(2.*s_alpha_v.*delta_we);
%产生0~2*pi之间均匀分布的随机变量epsilon
y=zeros(25,1);y(1)=1977326743;b=7;M=10^10;
for i=1:25
y(i+1)=mod(7*y(i),M);
x(i)=y(i)/M;
end
epsilon=2*pi*x;
alpha=0;alpha_v=0;
t=0:0.1:30*pi;
for i=1:25
a=s_alpha_ap(i).*sin(w(i)*t+epsilon(i));
alpha=alpha+a;
aa=s_alpha_v_ap(i).*sin(we(i)*t+epsilon(i));
alpha_v=alpha_v+aa;
end
miu_roll=0.1614;%横摇阻尼，可以变化，0~2；
omg_roll=0.8;%横摇固有频率给定？根据船型计算？
A=we./omg_roll;%遭遇频率和固有频率之比
W_roll_2=1./((1-A.^2).^2+4.*(miu_roll^2).*(A.^2));%幅频响应函数
s_roll=W_roll_2.*s_alpha_v;%横摇角能谱,用s_alpha_v代替s_alpha
%plot(w,s_roll)
s_roll_ap=s_alpha_v_ap.*sqrt(W_roll_2);
epsilon_roll=atan(2.*miu_roll.*A./(1-A.^2));
hyj=0;
for i=1:25
r=s_roll_ap(i).*sin(we(i)*t-epsilon_roll(i));
hyj=hyj+r;
end
hold on;
figure(1),plot(t,10.*alpha,'r',t,10*alpha_v,'b');xlabel('t (s)');ylabel('\alpha'); %绘制波倾角时域图
legend('波倾角','有效波倾角');
%plot(t,20.*alpha_v);xlabel('t (s)');ylabel('\alpha_v');
figure(2),plot(t,10.*hyj);xlabel('t (s)');ylabel('横摇角'); %绘制横摇角时域图
%subplot(3,1,3);
%plot(t,20.*hyj);xlabel('t (s)');ylabel('横摇角'); %绘制横摇角时域图