读研之掉进故障检测（四）—EMD经验模态分解（CWRU轴承数据集）_轴承故障诊断传统emd分解

读研之掉进故障检测（四）—EMD经验模态分解（CWRU轴承数据集）

---

文章目录

• 读研之掉进故障检测（四）—EMD经验模态分解（CWRU轴承数据集）
• • @[TOC](文章目录)
• 前言
• 一、EMD经验模态分解
• • 1.1EMD概念
  • 1.2EMD算法
• 二、实现代码
• 三、不依赖库的EMD（未能实现）

前言

读研之掉进故障检测（四）—EMD经验模态分解（CWRU轴承数据集）
读研之掉进故障检测（三）为论文笔记，目前暂不发布。
但是该博文，将对论文中的EMD经验模态分解做相应的概念阐述，以及代码实现
参考论文：[1] Magrin-Chagnolleau I , Baraniuk R G .Empirical mode decomposition based frequency attributes[J].ann.internat.mtg.soc.exp.geophys.expanded abstracts, 1999.DOI:10.1190/1.1820932.
[2] Ne. H , Sr. L , Mlc. W ,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society. Mathematical, physical and engineering sciences, 1998(1971):454.

一、EMD经验模态分解

1.1EMD概念

------经验模态分解（EMD）是Huang[2] 等人提出的一种新的模态分解方法。该技术自适应地将信号分解为振荡分量。通过该方法，任何复杂的数据集都可以分解为有限的且通常是少量的“固有模式函数”，这些函数允许良好的希尔伯特变换。不同的分量与信号本身非常匹配，并且从分解中提取的时频表示（称为Hilbert幅度谱）对时频平面的平铺方式没有任何约束。

1.2EMD算法

------该技术的原理是将信号x(t)分解为函数之和。该函数具有两个特点（1）具有相同的过零点和极值数目、（2）关于局部平均值对称。第一条件类似于平稳高斯过程的窄带要求。第二个条件将全局要求修改为局部要求，并且是必要的，以确保瞬时频率不会有不希望的波动，如由不对称波形引起的。这些函数被称为本征模式函数（IMF），并表示为$im{f}_i(t)$。而具体获得本征函数的方法，可以由以下步骤所示。$x(t)$是原信号
1、初始化：$r_o(t)=x(t),i=1$
2、提取第i个IMF：
（a）初始化：$h_o(t)=r_i(t),j=1$
（b）提取$h_{j-1}(t)$的局部极小值和极大值。
（c）用三次样条插值局部最大值和局部最小值，以形成$h_{j-1}(t)$的的上下包络。
（d）计算上下包络线的平均值$m_{\mathrm{j}-1}(t)$
（e）$h_j(t)=h_{j-1}(t)-m_{j-1}(t)$
（f）如果满足停止标准，则设置$im{f}_i(t)=h_j(t)$，否则转到（b），$j=j+1$
3、$r_i(t)=r_{i-1}(t)-im{f}_i(t)$
4、如果$r_i(t)$仍然具有至少2个极值，则在$i=i+1$的情况下转到2，否则分解完成，并且$r_i(t)$是余下的。

则在算法完成后，我们能获得 其中$r_n(t)$是分解时的剩余部分。
$x(t)=\sum _{i=1}^{n}im{f}_i(t)+r_n(t)$

------除了以上的算法解释外，另一种解释经验模式分解如何工作的方法是，它挑选出信号中保留的最高频率振荡。因此，局部地，每个IMF包含比之前提取的IMF更低的频率振荡。这个属性对于提取频率变化非常有用，因为变化在IMF的水平上会表现得更清楚。

二、实现代码

数据为CWRU轴承数据集中12k Drive End Bearing Fault Data的B007_0.mat
频率为12K

具体代码如下：

import scipy.io
import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD

mat = scipy.io.loadmat('F:/pycharm_dataset/CWRU_offcial/12k Drive End Bearing Fault Data/B007_0.mat')

signal = mat['X118_DE_time'].flatten()  
rpm = mat['X118RPM'][0][0]  


fs = 12000  
t = np.arange(len(signal)) / fs  


start_time = 0  
end_time = 0.5  
mask = (t >= start_time) & (t <= end_time)
selected_signal = signal[mask]
selected_t = t[mask]


emd = EMD()
imfs = emd(selected_signal)


plt.figure(figsize=(10, 6))
for i in range(5):
    plt.subplot(5, 1, i+1)
    plt.plot(selected_t, imfs[i])
    plt.title(f'IMF {i+1}')

plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

---

结果为：

三、不依赖库的EMD（未能实现）

这里给出一个写了半个下午，但是一直错的代码，错误的原因在stopping_criterion（）处理h时总是出错，导致一直没有imf分量获得。
stopping_criterion（）的原理是对数据序列中除了最后两个元素和第一个元素之外的所有元素进行相邻元素相乘的操作，并判断结果是否小于零。如果某一相邻元素的乘积小于零，则表示这两个元素之间发生了符号变化，即可能存在极值点。通过统计满足这个条件的相邻元素对的数量，可以得到一个极值点的估计值。
可能的错误原因是平坦区域或连续出现多个极值点的情况

import scipy.io
import numpy as np
import matplotlib.pyplot as plt

mat = scipy.io.loadmat('F:/pycharm_dataset/CWRU_offcial/12k Drive End Bearing Fault Data/B007_0.mat')

signal = mat['X118_DE_time'].flatten()
rpm = mat['X118RPM'][0][0]

fs = 12000
t = np.arange(len(signal)) / fs

start_time = 0
end_time = 0.5
mask = (t >= start_time) & (t <= end_time)
selected_signal = signal[mask]
selected_t = t[mask]


# Empirical Mode Decomposition (EMD)
def emd(signal, n):
    r = signal.copy()
    imfs = []

    def stopping_criterion(data):
        return len(np.argwhere(data[:-2] * data[1:-1] < 0)) <= 1

    while not stopping_criterion(r) and len(imfs) < n:
        h = r.copy()
        extrema = []

        while True:
            maxima = np.array([i for i in range(1, len(h) - 1) if (h[i - 1] < h[i] and h[i] > h[i + 1])])
            minima = np.array([i for i in range(1, len(h) - 1) if (h[i - 1] > h[i] and h[i] < h[i + 1])])

            maxima = np.concatenate(([0], maxima, [len(h) - 1]))
            minima = np.concatenate(([0], minima, [len(h) - 1]))

            max_env = np.interp(range(len(h)), maxima, h[maxima])
            min_env = np.interp(range(len(h)), minima, h[minima])

            mean_env = (max_env + min_env) / 2

            h = h - mean_env

            if stopping_criterion(h):
                imfs.append(h)
                break

        r = r - h

    imfs.append(r)
    return imfs

imfs = emd(selected_signal,5)

# Plot the first five IMF components
plt.figure(figsize=(10, 6))
for i in range(5):
    plt.subplot(5, 1, i + 1)
    plt.plot(selected_t, imfs[i])
    plt.title(f'IMF {i + 1}')

plt.tight_layout()
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66