基于UIC513标准的乘坐舒适度算法MATLAB代码

function[U]=UIC513（ Matrix ）   %Y.S.J ver.2002.9.27
%------------------------------
%close all;
clc;
Matrix = Matrix‘;

frs=250;

waf1=0.4;
waf2=100;
waq1=0.71;
wag1=tf（[4*pi*pi*waf2*waf200][12*pi*waf1/waq14*pi*pi*waf1*waf1]）;
wag2=tf（[1][12*pi*waf2/waq14*pi*pi*waf2*waf2]）;
wag=wag1*wag2;

%H（a）传函系数

[wabwaa]=tfdata（wag‘v‘）;
wab;
waa;
[wabzwaaz]=bilinear（wabwaafrs）;
wabz;
waaz;

w=logspace（03）;
f=w/（2*pi）;
wah=freqs（wabwaaw）;
wamag=abs（wah）;
waphase=angle（wah）;
%figure（1）
%subplot（211）loglog（fwamag）
%subplot（212）loglog（fwaphase）
%pause;
w=logspace（03）;
f=w/（2*pi）;
wah=freqz（wabzwaazw）;
wamag=abs（wah）;
waphase=angle（wah）;
%figure（2）
%subplot（211）loglog（fwamag）
%subplot（212）loglog（fwaphase）
%pause;

Wbf3=16;
Wbf4=16;
Wbf5=2.5;
Wbf6=4;
Wbq2=0.63;
Wbq3=0.8;
Wbq4=0.8;
Wbk=0.4;

%depict the transfer function of Hb
wbg1=tf（[2*pi*Wbk*Wbf4*Wbf4*Wbf6*Wbf6/Wbf34*pi*pi*Wbk*Wbf4*Wbf4*Wbf6*Wbf6][12*pi*Wbf4/Wbq24*pi*pi*Wbf4*Wbf4]）;
wbg2=tf（[1/（Wbf5*Wbf5）2*pi/（Wbf5*Wbq3）4*pi*pi][12*pi*Wbf6/Wbq44*pi*pi*Wbf6*Wbf6]）;
wbg=wbg1*wbg2;

%H（b）传函系数

[wbbwba]=tfdata（wbg‘v‘）;
[wbbzwbaz]=bilinear（wbbwbafrs）;

[HW]=freqz（wbbzwbaz）;
%loglog（W*243/（2*pi）abs（H））;grid

w1=logspace（03）;
wbh=freqs（wbbwbaw1）;
wbmag=abs（wbh）;
%wbmag=abs（wbh-wah）
wbphase=angle（wbh）;
f1=w1/（2*pi）;
%figure（2）
%loglog（f1wbmag）
%pause;
%loglog（f1wbphase）

%depict the transfer function of Hc
Wcf3=8;
Wcf4=8;
Wcq2=0.63;
Wck=1.0;

%H（c）传函系数

Wcb=[2*pi*Wck*Wcf4*Wcf4/Wcf34*pi*pi*Wck*Wcf4*Wcf4];
Wca=[12*pi*Wcf4/Wcq24*pi*pi*Wcf4*Wcf4];
[wcbzwcaz]=bilinear（WcbWcafrs）;

w3=logspace（03）;
wch=freqs（WcbWcaw3）;
wcmag=abs（wch）;
f3=w3/（2*pi）;
%mag = 20*log10（mag）;
%figure（3）
%subplot（211） loglog（f3wcmag）
%pause;

%depict the transfer function of Hd
Wdf3=2;
Wdf4=2;
Wdq2=0.63;
Wdk=1.0;

%H（d）传函系数

Wdb=[2*pi*Wdk*Wdf4*Wdf4/Wdf34*pi*pi*Wdk*Wdf4*Wdf4];
Wda=[12*pi*Wdf4/Wdq24*pi*pi*Wdf4*Wdf4];
[wdbzwdaz]=bilinear（WdbWdafrs）;

w4=logspace（03）;
wdh=freqs（WdbWdaw4）;
wdmag=abs（wdh）;
f4=w4/（2*pi）;
%figure（4）
% %subplot（211） loglog（f4wdmag）
% fs=243;

% %  40Hz巴特沃斯低通滤波器参数计算
% Wp=40;Ws=45;rp=0.3;rs=30;F=fs;   %fs=243
% % Wp=40;Ws=45;rp=0.3;rs=40;F=frs;   %fs=243--973
% [NWn]=buttord（Wp/（F/2）Ws/（F/2）rprs）;
% [ba]=butter（NWn）;
f