阶次分析变转速机械故障诊断求助

%~~~~~~~~~~~~~导入数据部分~~~~~~~~~~~~~%

clc

clear

fs=10240;

load s1.mat;

load rpm.mat;

S=s1(1:1048575);

S2=rpm(1:1048575);

t=(0:length(S2)-1)/fs;

array_time_amp=S; %导入时域振动信号

pluse=S2; %导入脉冲信号

% %~~~~~~~~~~~~~~~计算脉冲发生时刻部分~~~~~~~~~~~~~~~%

% Num_pluse1=1;

Num_pluse1=[];

Threshold=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)<Threshold&&pluse(temp2+1)>=Threshold)

Num_pluse1=[ Num_pluse1,temp2+1];%Num_pluse1=[1,219,6,7.....temp(2+1)]

end

end

t_pluse=(Num_pluse1-1)/fs;%找到脉冲发生的前沿时刻

%

%~~~~~~~~~~~~~~~等角度时间计算~~~~~~~~~~~~~~~~%

%delt_thet=pi/24;

%t_angle=[];

f_pluse=[];

f_pluse(1)=t_pluse(1);

for temp3=1:length(t_pluse)-1

f_pluse(temp3+1)=[1/(t_pluse(temp3+1)-t_pluse(temp3))];

end

speed=f_pluse*60;

thet=cumtrapz(t_pluse,f_pluse);

angle=thet*2*pi;

delt_thet=pi/1000;

%%%%求取对应的时刻点

maxjd=max(angle);

maxk=fix(maxjd/delt_thet); %让x向0靠近取整 （需要离散多少步）

minjd=min(angle);

mink=ceil(minjd/delt_thet); %朝正无穷大取整

k_thet=(mink:maxk)*delt_thet;

t_angle=spline(angle,t_pluse,k_thet); %三次样条函数插值，根据己知的x，y数据，用三次样条函数插值得到拟合曲线，之后计算出xi处的值yi；

array_angle=[1:length(t_angle)].*delt_thet;

array_angle_amp=interp1(t,array_time_amp,t_angle,'spline'); %角域重采样后的信号

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;

figure(2);

subplot(2,1,1),plot(array_angle,array_angle_amp), title('angle dominant振动仿真信号角域图'), xlabel('角度angle /rad'),

ylabel('幅值amplitude');

grid on;

subplot(2,1,2),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');

xlim([0,30]);

grid on;