基于MATLAB实现的uic513标准的算法振动，分垂横纵3个方向，输入，经MATLAB算法得出振动舒适度的指标+使用说明文档

function[ride]=yallmax %Y.S.J ver.2002.9.27 %------------------------------ %close all; clc; frs=2500; waf1=0.4; waf2=100; waq1=0.71; wag1=tf([4*pi*pi*waf2*waf2,0,0],[1,2*pi*waf1/waq1,4*pi*pi*waf1*waf1]); wag2=tf([1],[1,2*pi*waf2/waq1,4*pi*pi*waf2*waf2]); wag=wag1*wag2; %H(a)传函系数 [wab,waa]=tfdata(wag,'v'); wab; waa; [wabz,waaz]=bilinear(wab,waa,frs); wabz; waaz; w=logspace(0,3); f=w/(2*pi); wah=freqs(wab,waa,w); wamag=abs(wah); waphase=angle(wah); %figure(1) %subplot(2,1,1),loglog(f,wamag) %subplot(2,1,2),loglog(f,waphase) %pause; w=logspace(0,3); f=w/(2*pi); wah=freqz(wabz,waaz,w); wamag=abs(wah); waphase=angle(wah); %figure(2) %subplot(2,1,1),loglog(f,wamag) %subplot(2,1,2),loglog(f,waphase) %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/Wbf3,4*pi*pi*Wbk*Wbf4*Wbf4*Wbf6*Wbf6],[1,2*pi*Wbf4/Wbq2,4*pi*pi*Wbf4*Wbf4]); wbg2=tf([1/(Wbf5*Wbf5),2*pi/(Wbf5*Wbq3),4*pi*pi],[1,2*pi*Wbf6/Wbq4,4*pi*pi*Wbf6*Wbf6]); wbg=wbg1*wbg2; %H(b)传函系数 [wbb,wba]=tfdata(wbg,'v'); [wbbz,wbaz]=bilinear(wbb,wba,frs); [H,W]=freqz(wbbz,wbaz); %loglog(W*243/(2*pi),abs(H));grid w1=logspace(0,3); wbh=freqs(wbb,wba,w1); wbmag=abs(wbh); %wbmag=abs(wbh-wah) wbphase=angle(wbh); f1=w1/(2*pi); %figure(2) %loglog(f1,wbmag) %pause; %loglog(f1,wbphase) %depict the transfer function of Hc Wcf3=8; Wcf4=8; Wcq2=0.63; Wck=1.0; %H(c)传函系数 Wcb=[2*pi*Wck*Wcf4*Wcf4/Wcf3,4*pi*pi*Wck*Wcf4*Wcf4]; Wca=[1,2*pi*Wcf4/Wcq2,4*pi*pi*Wcf4*Wcf4]; [wcbz,wcaz]=bilinear(Wcb,Wca,frs); w3=logspace(0,3); wch=freqs(Wcb,Wca,w3); wcmag=abs(wch); f3=w3/(2*pi); %mag = 20*log10(mag); %figure(3) %subplot(2,1,1), loglog(f3,wcmag) %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/Wdf3,4*pi*pi*Wdk*Wdf4*Wdf4]; Wda=[1,2*pi*Wdf4/Wdq2,4*pi*pi*Wdf4*Wdf4]; [wdbz,wdaz]=bilinear(Wdb,Wda,frs); w4=logspace(0,3); wdh=freqs(Wdb,Wda,w4); wdmag=abs(wdh); f4=w4/(2*pi); %figure(4) % %subplot(2,1,1), loglog(f4,wdmag) % 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 % [N,Wn]=buttord(Wp/(F/2),Ws/(F/2),rp,rs); % [b,a]=butter(N,Wn); fs=2500; wp=2*80*pi/fs; wdelta=wp; N=12; wn=wp/pi; b=fir1(N,wn,hanning(N+1)); disp(' '); disp(' '); disp('------------------------------------------------------'); disp('| 舒 适 度 计 算 程 序 |'); disp('| 2004.10 V2.0 |'); disp('| |'); disp('| SWJTU_TPL_YSJ |'); disp('------------------------------------------------------'); disp(' '); disp(' 计算每段5秒数据的三向振动舒适度 '); disp(''); for k=1:3 [filename,pathname,fileIndex]=uigetfile('*.txt','pick a file'); if fileIndex~=0 datafile=[pathname,filename]; if exist(datafile,'file'); fid=fopen(datafile,'r'); sy=fscanf(fid,'%g'); fclose(fid); end end [m n]=size(sy); sy=sy*9.81; sy=sy-mean(sy); data=sy; %data(:,1)=sy; %for j=1:2 %[filename,pathname,fileIndex]=uigetfile('*.txt','pick a file'); %if fileIndex~=0 % datafile=[pathname,filename]; % if exist(datafile,'file'); % fid=fopen(datafile,'r'); % sy=fscanf(fid,'%g'); % fclose(fid); % end %end %m1=length(sy); m1=length(data); if m1<m i=m-m1; for ii=1:i sy=[sy;0]; end; else if m1>m i=m1-m; sy=sy(1:i); end; end; sy=sy-mean(sy); %data(:,j+1)=sy; data=sy; % for k=1:3 data_1=fftfilt(b,data); % 40Hz巴特沃斯低通滤波器 滤波 data_a=filter(wabz,waaz,data_1); %Wa 计权 if k==1 data_end=filter(wbbz,wbaz,data_a); % Wb 垂向滤波 else if k==2 data_end=filter(wdbz,wdaz,data_a);% Wd 横向滤波 else data_end=filter(wdbz,wdaz,data_a);%Wc 纵向滤波 end; end; tfield=fix(m/(5*fs)); data_end=data_end(1:tfield*5*fs); datam=reshape(data_end,fs*5,tfield); %512*65_ [x,y]=size(datam); for j=1:y datam(:,j)=datam(:,j)-mean(datam(:,j)); end; mean_sq_r_v=sqrt(diag(cov(datam))); data_v(:,k)=mean_sq_r_v; yy(:,k)=data_v(:,k); % yy(k)=prctile(data_v,95); %计算各列均方根值的95％值 end; Azp95=diag(yy(:,1)*yy(:,1)'); Ayp95=diag(yy(:,2)*yy(:,2)'); Axp95=diag(yy(:,3)*yy(:,3)'); % leq_z=20*log10(yy(1)/10^(-6)); % leq_y=20*log10(yy(2)/10^(-6)); % leq_x=20*log10(yy(3)/10^(-6)); N=6*sqrt(Azp95+Ayp95+Axp95); %计算综合舒适度N=6*sqrt(x^2+y^2+z^2) % yy1=prctile(data_v,50); % ty=(5.6/3.15/yy1(lateral_acc))^2*0.167; %y(横向)向疲劳时间(h) % tz=(5.6/3.15/yy1(vertical_acc))^2*0.167; %z(垂向)向疲劳时间(h) % t=(5.6/3.15/(sqrt(1.4^2*yy1(lateral_acc)^2+yy(vertical_acc)^2))) % ^2*0.167; %综合疲劳时间(h) m=length(filename); createfile=strcat(pathname,'\'); createfile=char(createfile); if exist(createfile,'file')==0 createfile=strcat(pathname); mkdir(createfile); end createfile=strcat(pathname,'\'); savename=strcat(createfile,'\',filename(1:m-4),'_s','.txt'); fid=fopen(savename,'wt'); fprintf(fid,'%g\n',N); fclose(fid); clear data_v; % disp(' '); % ask_c=input('继续计算吗(是=y，否=n)？ ','s'); % % if ask_c=='y' sqrtd2 % end; % disp(' '); % disp('程序运行结束，谢谢使用！');