信号处理-基于希尔伯特解调（包络谱）的轴承故障诊断实战，通过python代码实现超详细讲解

一个包括内圈、外圈、滚动体和正常数据。
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期

以下是轴承内圈故障数据的分析。原始数据为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第四篇-包络谱/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] 

希尔伯特解调法，又称包络谱分析。

设$x(t)$为一个实时域信号，其Hilbert变换定义为：
$h(t)=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }\frac{x(\tau )}{t-\tau }\mathrm{d}\tau =x(t)\ast \frac{1}{\pi t}$
则原始信号$x(t)$和它的Hilbert变换信号$h(t)$可以构建一个新的解析信号$z(t)$:
$z(t)=x(t)+jh(t)=a(t)e^{j\phi t}$

# step1: 做希尔伯特变换
ht = fftpack.hilbert(xt)
print(ht)

输出结果：

[-0.02520403 -0.28707983 -0.00610516 ...  0.1100125   0.22821944
 -0.11203138]

此时输出的$h(t)$是解析信号$a(t)$的虚部系数

对$z(t)$取模，得到其幅值$a(t):$

$a(t)=|z(t)|=\sqrt{x^2(t)+h^2(t)}$

注：$a(t)$即为包络信号，也叫解析信号

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)

观察发现：
(1)在频率为0Hz的地方幅值比较大
(2)在低频部分貌似看到一倍频和二倍频

在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)

输出结果：

内圈故障特征频率：$F_{\mathrm{BPFI}}=\frac{n{f}_r}{2}\left(1+\frac{d}{D}\cos \alpha \right)$

外圈故障特征频率：$F_{\mathrm{BPFO}}=\frac{n{f}_r}{2}\left(1-\frac{d}{D}\cos \alpha \right)$

滚动体故障特征频率：$F_{\mathrm{BSF}}=\frac{D{f}_r}{2d}\left[1-{\left(\frac{d}{D}\cos \alpha \right)}^2\right]$

$n$: 滚动体个数，$f_r$: 轴转速 $d$: 滚珠(子)直径 $D$: 轴承节径

轴承型号为：6205-2RSL JME SKF 深沟球滚珠轴承
$d$=7.94mm, $D$=39.04mm, $\alpha$=0, $n$=9

为了方便，定义了一个轴承故障特征频率计算函数，只需输入参数$d$, $D$, $\alpha$, $n$和$f_r$即可

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')  # 二倍频

输出结