希尔伯特解调（包络谱）python代码实战及详细讲解，在CWRU数据上验证
欢迎关注公众号《故障诊断与python学习》
代码位置：https://github.com/HappyBoy-cmd/fault_diagnosis_signal_processing
参考资料：
机械故障诊断及典型案例解析（第2版，时献江）
会议论文：Bearing Intelligent Fault Diagnosis Under Complex Working Condition Based on SK-ES-CNN，2021 Global Reliability and Prognostics and Health Management (PHM-Nanjing)
1、数据介绍
包括内圈、外圈、滚动体和正常数据，分别为一个。
1730_12k_0.007-InnerRace：内圈故障
1730_12k_0.007-OuterRace3：外圈故障
1730_12k_0.014-Ball：滚动体故障
1730_48k_Normal：正常数据
对CWRU轴承数据集不了解的同学见这里：
CWRU数据集介绍 第1期
CWRU数据集介绍 第2期
CWRU数据集介绍 第3期
CWRU数据集介绍 第3期
2、加载CWRU内圈故障数据
下面先以轴承内圈故障数据进行分析。原始数据为mat文件，是matlab文件，定义一个函数进行数据读取
def data_acquision(FilePath):
"""
fun: 从cwru mat文件读取加速度数据
param file_path: mat文件绝对路径
return accl_data: 加速度数据，array类型
"""
data = scio.loadmat(file_path) # 加载mat数据
data_key_list = list(data.keys()) # mat文件为字典类型，获取字典所有的键并转换为list类型
accl_key = data_key_list[3] # 获取'X108_DE_time'
accl_data = data[accl_key].flatten() # 获取'X108_DE_time'所对应的值，即为振动加速度信号,并将二维数组展成一维数组
return accl_data
data_acquision(FilePath)输入参数FilePath，输出一个1维array数据。下面进行演示
file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_12k_0.007-InnerRace.mat'
xt = data_acquision(file_path)
plt.figure(figsize=(12,3))
plt.plot(xt)
print(xt)
输出结果：
[ 0.22269856 0.09323776 -0.14651649 ... -0.36125573 0.31138814
0.17055689]
3、希尔伯特解调(包络谱)分析
希尔伯特解调法，亦叫包络谱分析。
3.1希尔伯特黄变换
设为一个实时域信号，其Hilbert变换定义为：
则原始信号和它的Hilbert变换信号
可以构建一个新的解析信号
:
# step1: 做希尔伯特变换
ht = fftpack.hilbert(xt)
print(ht)
输出结果：
[-0.02520403 -0.28707983 -0.00610516 ... 0.1100125 0.22821944
-0.11203138]
此时输出的是解析信号
的虚部系数
对取模，得到其幅值
注：即为包络信号，也叫解析信号
3.2获得包络信号
t = np.sqrt(ht**2+xt**2) # at = sqrt(xt^2 + ht^2)
接下来对包络信号做fft即为包络谱
3.3获得包络谱
sampling_rate = 12000
am = np.fft.fft(at) # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am) # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)] # 获取正频率
plt.plot(freq, am)
3.4去直流分量
在0Hz的幅值比较高，使得其它频率幅值较低，不便观察。这种现象叫直流分量，去直流分量方法，y = y-mean(y)
sampling_rate = 12000
at = at - np.mean(at) # 去直流分量
am = np.fft.fft(at) # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am) # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)] # 获取正频率
plt.plot(freq, am)
输出结果：
sampling_rate = 12000
at = at - np.mean(at) # 去直流分量
am = np.fft.fft(at) # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am) # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)] # 获取正频率
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)
输出结果：
4、计算故障特征频率
内圈故障特征频率：
外圈故障特征频率：
滚动体故障特征频率：
: 滚动体个数，
: 轴转速
: 滚珠(子)直径
: 轴承节径
轴承型号为：6205-2RSL JME SKF 深沟球滚珠轴承
=7.94mm,
=39.04mm,
=0,
=9
4.1定义一个轴承故障特征频率计算函数
为了方便，定义了一个轴承故障特征频率计算函数，只需输入参数,
,
,
和
即可
def bearing_fault_freq_cal(n, d, D, alpha, fr=None):
'''
基本描述：
计算滚动轴承的故障特征频率
详细描述：
输入4个参数 n, fr, d, D, alpha
return C_bpfi, C_bpfo, C_bsf, C_ftf, fr
内圈 外圈 滚针 保持架 转速
Parameters
----------
n: integer
The number of roller element
fr: float(r/min)
Rotational speed
d: float(mm)
roller element diameter
D: float(mm)
pitch diameter of bearing
alpha: float(°)
contact angle
fr:：float(r/min)
rotational speed
Returns
-------
BPFI: float(Hz)
Inner race-way fault frequency
BPFO: float(Hz)
Outer race-way fault frequency
BSF: float(Hz)
Ball fault frequency
FTF: float(Hz)
Cage frequency
'''
C_bpfi = n*(1/2)*(1+d/D*np.math.cos(alpha))
C_bpfo = n*(1/2)*(1-(d/D)*np.math.cos(alpha))
C_bsf = D*(1/(2*d))*(1-np.square(d/D*np.math.cos(alpha)))
C_ftf = (1/2)*(1-(d/D)*np.math.cos(alpha))
if fr!=None:
return C_bpfi*fr/60, C_bpfo*fr/60, C_bsf*fr/60, C_ftf*fr/60, fr/60
else:
return C_bpfi, C_bpfo, C_bsf, C_ftf, fr
下面计算CWRU在转速1730rpm时的故障特征频率
bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
print('内圈故障特征频率',bpfi)
print('外圈故障特征频率',bpfo)
print('滚动体故障特征频率',bsf)
print(ftf)
print(fr)
5、理论故障特征频率与实际故障特征频率验证
sampling_rate = 12000
at = at - np.mean(at) # 去直流分量
am = np.fft.fft(at) # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am) # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)] # 获取正频率
plt.figure(figsize=(12,3))
plt.plot(freq, am)
plt.xlim(0,500)
plt.vlines(x=156.13, ymin=0, ymax=0.2, colors='r') # 一倍频
plt.vlines(x=156.13*2, ymin=0, ymax=0.2, colors='r') # 二倍频
输出结果：
6、与fft进行对比分析
那如果不希尔伯特变换解调，直接fft，能够看到故障特征频率吗？下面进行对比分析
sampling_rate = 12000
am = np.fft.fft(xt) # 对希尔伯特变换后的at做fft变换获得幅值
am = np.abs(am) # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)]
freq = np.fft.fftfreq(len(xt), d=1 / sampling_rate) # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)] # 获取正频率
plt.plot(freq, am)
输出结果：
plt.plot(freq, am)
plt.xlim(0, 500)
plt.vlines(x=156.13, ymin=0, ymax=0.05, colors='r') # 一倍频
plt.vlines(x=156.13*2, ymin=0, ymax=0.05, colors='r') # 二倍频
7、封装包络谱函数
为了更加方便使用包络谱，这里封装了一个包络谱函数plt_envelope_spectrum
def plt_envelope_spectrum(data, fs, xlim=None, vline= None):
'''
fun: 绘制包络谱图
param data: 输入数据，1维array
param fs: 采样频率
param xlim: 图片横坐标xlim，default = None
param vline: 图片垂直线，default = None
'''
#----去直流分量----#
data = data - np.mean(data)
#----做希尔伯特变换----#
xt = data
ht = fftpack.hilbert(xt)
at = np.sqrt(xt**2+ht**2) # 获得解析信号at = sqrt(xt^2 + ht^2)
am = np.fft.fft(at) # 对解析信号at做fft变换获得幅值
am = np.abs(am) # 对幅值求绝对值（此时的绝对值很大）
am = am/len(am)*2
am = am[0: int(len(am)/2)] # 取正频率幅值
freq = np.fft.fftfreq(len(at), d=1 / fs) # 获取fft频率，此时包括正频率和负频率
freq = freq[0:int(len(freq)/2)] # 获取正频率
plt.plot(freq, am)
if vline: # 是否绘制垂直线
plt.vlines(x=vline, ymax=0.2, ymin=0, colors='r') # 高度y 0-0.2，颜色红色
if xlim: # 图片横坐标是否设置xlim
plt.xlim(0, xlim)
plt.xlabel('freq(Hz)') # 横坐标标签
plt.ylabel('amp(m/s2)') # 纵坐标标签
7.1外圈故障数据测试
file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_12k_0.007-OuterRace3.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bpfo)
输出结果：
7.2滚动体故障数据测试分析
file_path = r'E:/研究生/pytorch/CSDN代码/fault_diagnosis_signal_processing/第4篇-包络谱/1730_12k_0.014-Ball.mat'
data = data_acquision(file_path)
plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bsf)
可见实际滚动体故障特征频率不明显
8、总结
(1)包络谱能够检测出内圈、外圈故障，滚动体比较困难
(2)直接使用fft难以检测故障特征频率，故障特征频率易被高频掩盖
文章出处登录后可见！
已经登录？立即刷新
