手把手教你安装GPU版Pytorch环境及手把手带大家跑通第一个故障诊断案例！！！（上）

本文分享阶次分析代码。阶次基本介绍等时间采样无法对变转速设备进行故障分析，由于转速的变化，导致在等时间采样的情况下，无法采集到同相位的点，所以直接进行FFT转换后，信号往往会出现非常大的泄露，如图1所示。图1 等时间采样频谱阶次分析首先要确定参考轴的旋转频率。假设参考轴的转速为 ，则其旋转频率 。某一振动信号的频率为 (除去固有频率之外，其他频率均会随着转速的变化而等比例变化)，那么该振动信号对应的阶次 ，所以无论转速如何变化，故障阶次始终不会改变。实现阶次谱存在两个关键步骤：(1).插值获取等角度间隔的时间点(2).使用等角度间隔的时间点对振动信号重采样得到等角度间隔的振动信号最后只需对重采样的振动信号FFT变换，即可得到阶次谱。需要注意角采样频率 仍需满足奈奎斯特采样定理，大于等于最大阶次2.56倍，阶次谱中最大阶次则为 。案例数据介绍采用渥太华升速轴承内圈损伤数据，数据为.mat数据，包含两类数据振动信号(Channel_1)和脉冲信号(Channel_2)，采样频率200000，每圈脉冲数1024，内圈故障阶次5.43。实现过程讲解读取原始振动数据与转速脉冲(局部放大)，如图2、图3。图2 原始数据图3 转速脉冲2.根据转速脉冲计算转频时间曲线如图4，对给定时间内求取旋转角度的累积变化量如图5。图4 转频时间曲线图5 角度时间曲线3.使用上述获取到的角度以及时间创建一个插值函数(python中使用interp1d实现)，将角度进行等间隔划分并带入，即可获取到等角度间隔的时间点。使用原始振动数据和时间再次创建一个插值函数，将等角度间隔的时间点带入，即可获取到等角度间隔采样的振动信号，如图6。图6 重采样振动信号4.对重采样的振动信号进行FFT转换，得到阶次谱如图7(进行局部放大)，但是故障阶次不够明显，进行包络解调得到包络阶次谱如图8(进行局部放大)，此时故障阶次5.43非常明显。图7 阶次谱图8 阶次包络谱代码import numpy as npfrom scipy.interpolate import interp1dimport matplotlib.pyplot as pltfrom scipy.fftpack import fftfrom scipy.integrate import cumulative_trapezoid, trapezoidfrom scipy.fftpack import hilbert, fft, ifftfrom matplotlib.pylab import mplfrom scipy import fftpackimport pandas as pdimport scipympl.rcParams['font.sans-serif'] = ['SimHei']mpl.rcParams['axes.unicode_minus'] = Falsedef trend_remove(x, y): # 用2次多项式拟合x，y数组,返回多项式系数 param = np.polyfit(x, y, 2) # 拟合完之后，用生成的多项式系数用来生成多项式函数 y_poly_func = np.poly1d(param) # 生成多项式函数之后，就是获取x在这个多项式处的值 y_poly = y_poly_func(x) # 原始信号减去拟合函数的值即为去除趋势项后的信号 de_trend_y = y - y_poly return de_trend_yfile = r'C:\Users\Deng\Desktop\新建文件夹 (2)\I-A-1.mat'dataset = scipy.io.loadmat(file)vib = np.array([i[0] for i in dataset['Channel_1']])pulse = np.array([i[0] for i in dataset['Channel_2']])fs = 200000t = (1 / fs) * np.arange(0, len(vib))# 振动信号去除趋势项vib = trend_remove(t, vib)# 转速计算：SpeedFrePulseTime = []for i in range(len(pulse)-1): if pulse[i + 1] < 2 < pulse[i]: PulseTime.append(t[i])PulseTime = np.array(PulseTime)SpeedFre = []speedTime = []for i in range(0, len(PulseTime) + 1, 1024): if i == 0: SpeedFre.append(0) else: Fre = 1 / (PulseTime[i] - PulseTime[i - 1024]) SpeedFre.append(Fre) speedTime.append(PulseTime[i])SpeedFre = np.array(SpeedFre)SpeedFre = SpeedFre[1:]PulseTime = PulseTime[1:]# 积分求取转过的总角度AngleAngle = cumulative_trapezoid(SpeedFre*360, speedTime, initial=0)plt.figure(figsize=(11, 4))plt.plot(speedTime, SpeedFre)plt.title('转速')plt.xlabel('时间')plt.ylabel('hz')plt.figure(figsize=(11, 4))plt.plot(speedTime, SpeedFre*360)plt.title('角度')plt.xlabel('时间')plt.tight_layout()plt.show()# 获取等角度间隔的时间点：EqualAngleTimeangleStep = 0.1EqualAngle = np.arange(0, trapezoid(SpeedFre*360, speedTime), angleStep)FunctionA = interp1d(Angle, speedTime)EqualAngleTime = FunctionA(EqualAngle)# 用等角度间隔的时间点对振动信号重新采样，得到了等角度间隔的振动信号：EqualAngleVibFunctionV = interp1d(t, vib)EqualAngleVib = FunctionV(EqualAngleTime)# 对重采样信后进行fft变换VibFft = (np.abs(fft(EqualAngleVib)) * 2) / len(EqualAngleVib)VibFftHalf = VibFft[1:int(np.round(len(EqualAngleVib) / 2))]order = np.arange(0, len(EqualAngleVib), dtype='int64') * (1/angleStep*360) / len(EqualAngleVib)orderHalf = order[1:int(np.round(len(EqualAngleVib) / 2))]plt.figure(figsize=(11, 4))plt.plot(EqualAngle, EqualAngleVib)plt.title('重采样之后的振动信号')plt.xlabel('角度')plt.ylabel('振动幅值')plt.figure(figsize=(11, 4))plt.plot(orderHalf, VibFftHalf)plt.title('阶次谱')plt.xlabel('阶次')plt.ylabel('振动幅值')plt.xlim(0, 100)plt.show()angleFs = (1/angleStep)*360orderMin = 1orderMax = 500# 带通滤波后逆变换signal_fft = fftpack.fft(EqualAngleVib)sig_fre = fftpack.fftfreq(n=len(signal_fft), d=1 / angleFs)signal_fft[(abs(sig_fre) < orderMin) | (abs(sig_fre) > orderMax)] = 0signal_filter = fftpack.ifft(signal_fft).real# 包络变换H = np.abs(hilbert(signal_filter) - np.mean(signal_filter))HP = np.abs(fft(H - np.mean(H))) * 2 / len(EqualAngleVib)HP_half = HP[1:int(np.round(len(EqualAngleVib) / 2))]plt.figure(figsize=(11, 4))plt.plot(orderHalf, HP_half)plt.title('阶次包络谱')plt.xlabel('阶次')plt.ylabel('振动幅值')plt.xlim(0, 100)plt.show() 编辑：海洋校核：李正平、陈凯歌、赵栓栓、曹希铭、赵学功、白亮、任超、陈宇航、Tina、王金、陈莹洁该文资料搜集自网络，仅用作学术分享，不做商业用途，若侵权，后台联系小编进行删除来源：故障诊断与python学习