order_变转速_变转速故障_order在matlab_脉冲转速_matlab中order_

%~~~~~~~~~~~~~导入数据~~~~~~~~~~~~~% clc clear Fs=71680; N=Fs*5; t=(0:N-1)/Fs; Adc = 1; %直流分量幅度 S=sin(2*pi*t.^2- pi/6)+sin(4*pi*t.^2- pi/6)+sin(8*pi*t.^2- pi/6) S2=Adc+ sin(2*pi* t.^2- pi/6);%参考轴的转速为n(t)=60t r/min % S=sin(2*pi*t.^2)+sin(4*pi*t.^2)+sin(8*pi*t.^2); % S2=Adc+ sin(2*pi* t.^2); % S = sin(2*pi*t.^2-pi/6)+sin(8*pi*t.^2+pi/6)+sin(40*pi*t); % S = Adc+ sin(2*pi*t.^2-pi/6) ; % S2 = Adc+ 60 * t ; array_time_amp=S; %导入时域振动信号 pluse=S2; %导入脉冲信号 figure(1); % subplot(2,1,1),plot(t,array_time_amp),title('time dominant振动信号时域图'),xlabel('时间time'),ylabel('幅值amplitude'); subplot(2,1,1),plot(t,array_time_amp),title('变转速数据'),xlabel('时间time'),ylabel('幅值amplitude'); grid on; % subplot(2,1,2),plot(t,pluse),title('keyphasor键相脉冲仿真信号时域图'),xlabel('时间time'),ylabel('幅值amplitude'); subplot(2,1,2),plot(t,pluse),title('转速脉冲'),xlabel('时间time'),ylabel('幅值amplitude'); grid on; % %~~~~~~~~~~~~~~~计算脉冲发生时刻~~~~~~~~~~~~~~~% Num_pluse1=1; Num_pluse1=[]; Threshold=1.5;%设定脉冲阈值 for temp2=1:length(pluse)-1;%设temp2为步长为1的[1,2.....,n-1]n-1维数组；即将所有的采样点编号；length函数功能是返回pluse的数组维数；for循环体循环n-1次，再结束。 % for temp2=1:length(pluse)-1; if(pluse(temp2)<1.5&&pluse(temp2+1)>=1.5) Num_pluse1=[ Num_pluse1,temp2+1];%Num_pluse1=[1,219,6,7.....temp(2+1)] end end Num_pluse1=Num_pluse1(1:length(Num_pluse1)); if length(Num_pluse1)<2,return,end;%如果脉冲信号就一个数则跳出计算脉冲发生时刻的函数 t_pluse=(Num_pluse1-1)/Fs;% % % % % %~~~~~~~~~~~~~~~等角度时间计算~~~~~~~~~~~~~~~~% delt_thet=pi/24; t_angle=[]; for temp3=3:length(t_pluse); b=inv([1,t_pluse(temp3-2),t_pluse(temp3-2)^2;1,t_pluse(temp3-1),t_pluse(temp3-1)^2;1,t_pluse(temp3),t_pluse(temp3)^2])*[0,2*pi,4*pi]'; if temp3==3;%(如果数temp3等于3则执行第一个循环，如果temp3不等于3则执行第二个循环，即在脉冲时刻[tt0,tt3]执行第一个循环，等角度时刻记到[t1，t72]) k=0; while k<1.5*2*pi/delt_thet;%0<k<1.5*2*24=72 tt=(sqrt(4*b(3)*(k*delt_thet-b(1))+b(2)^2)-b(2))/(2*(b(3)+eps)); t_angle=[t_angle,tt]; k=k+1; end else k=pi/delt_thet;%k=24,(如果temp3不等于3则执行第二个循环，即在脉冲时刻[tt1,tt4].......[tt13,tt16]执行第二个循环，等角度时刻记到[t72，t120]，[t120,t168]........[t648,t696])) while k>=pi/delt_thet && k<1.5*2*pi/delt_thet;% 24<k<72时，执行循环。否则循环结束。 tt=(sqrt(4*b(3)*(k*delt_thet-b(1))+b(2)^2)-b(2))/(2*(b(3)+eps)); t_angle=[t_angle,tt]; k=k+1; end end end array_angle=[1:length(t_angle)].*delt_thet; % % ~~~~~~~~~~~~~~~~~~~~~~ 数字跟踪滤波器~~~~~~~~~~~~~~~~~~% i=2; while i<=length(t_pluse); ft=2/(t_pluse(i)-t_pluse(i-1));%[t_pluse(i-1),t_pluse(i)]时间段内的采样频率 [b,a]=ellip(6,3,80,ft*25/(Fs/2));%ft*order/(Fs/2)，归一化频率，其中1对应Fs/2 array_time_amp(Num_pluse1(i-1):(Num_pluse1(i)-1))=filter(b,a,array_time_amp(Num_pluse1(i-1):( Num_pluse1(i)-1)));%对[Num_pluse1(1),Num_pluse1(2)-1]内振动信号进行数字跟踪低通滤波 i=i+1; end %~~~~~~~~~~~~~~~~~~等角度时刻的幅值拟合~~~~~~~~~~~~~~~~% syms array_angle_amp; array_angle_amp=[]; array_angle_amp=interp1(t,array_time_amp,t_angle,'spline');%linear%spline%cubic %~~~~~~~~~~~~评价标准：幅值平方求能量谱~~~~~~~~~~~~~~~% nlfx=(0:length(t_angle)-1)*48/length(t_angle); nlpfz=(abs(fft(array_angle_amp))*2*(1/0.5)/length(array_angle_amp)).^2; pks = findpeaks(nlpfz(1:length(nlpfz)/2),'minpeakheight',0.1,'NPEAK',3); znl=sum(nlpfz); figure(5); plot(nlfx(1:length(nlfx)/2),nlpfz(1:length(nlpfz)/2)),title('time dominant振动信号时域图'); %------------------加窗--------------------------------% %-------------------------------------------------------% w=Chebwin(length(array_angle_amp))';%rectwin 矩形窗 Kaiser 凯瑟窗 Hamming 海明窗 hanning 汉宁窗 bartlett 三角窗 Gausswin 高斯窗 Blackman 布莱克曼窗 Chebwin 切比雪夫窗 array_angle_amp=array_angle_amp.*w; % ---------------------------------------------------------------% %-------------- FFT ---------------------------------------% angle_dom_ffty=abs(fft(array_angle_amp))*2/length(array_angle_amp); delt_order=2*pi/(length(angle_dom_ffty)*delt_thet); angle_dom_fx=(0:length(angle_dom_ffty)-1)*delt_order; FFTy=abs(fft(array_time_amp))*2/length(array_time_amp); fx=(0:length(array_time_amp)-1)*Fs/length(array_time_amp)/1000; %---------------2 维图形显示------------------------------- figure(2); subplot(2,1,1),plot(fx(1:length(fx)/2),FFTy(1:length(fx)/2)),title('Frequency domain频域'); xlabel('频率f'), ylabel('幅值amplitude'); title('振动仿真信号 幅度-频率曲线图'); subplot(2,1,2),plot(array_angle,array_angle_amp), title('angle dominant振动仿真信号角域图'), xlabel('角度angle /rad'), ylabel('幅值amplitude'); grid on; figure(3); subplot(1,1,1),plot(angle_dom_fx(1:length(angle_dom_fx)/2),angle_dom_ffty(1:length(angle_dom_fx)/2)), title('order dominant振动仿真信号的阶次域图'), xlabel('阶次order'),ylabel('幅值amplitude'); grid on;