信号处理基础 | 一文带你明白经验模态分解(EMD)与固有模态函数(IMF)，附python代码！

1年前浏览6907

本期给大家推荐一篇能够通俗易懂的理解经验模态分解(EMD)的文章，后期会继续推荐EMD的各种变种方法。

来源：知乎
作者：Mr.看海
链接：https://zhuanlan.zhihu.com/p/40005057

1 为什么要用EMD

在信号处理方面我们有时域处理方法（如有效值、峭度）、频域处理方法（如频谱、功率谱）以及一些时频域处理方法（如小波分析）。EMD（Empirical Mode Decomposition）作为时频域的处理方法，相对于同样是时频域方法的小波分析有什么好处呢？

EMD最显著的特点，就是其克服了基函数无自适应性的问题。啥意思呢？回忆小波分析部分的内容，我们会知道小波分析是需要选定某一个小波基的，小波基的选择对整个小波分析的结果影响很大，一旦确定了小波基，在整个分析过程中将无法更换，即使该小波基在全局可能是最佳的，但在某些局部可能并不是，所以小波分析的基函数缺乏适应性。

通俗的说，用EMD有什么好处呢？对于一段未知信号，不需要做预先分析与研究，就可以直接开始分解。这个方法会自动按照一些固模式按层次分好，而不需要人为设置和干预。

再通俗一点，EMD就像一台机器，把一堆混在一起的硬币扔进去，他会自动按照1元、5毛、1毛、5分、1分地分成几份。

2 固有模态函数

固有模态函数（Intrinsic Mode Functions, IMF）就是原始信号被EMD分解之后得到的各层信号分量。EMD的提出人黄锷认为，任何信号都可以拆分成若干个固有模态函数之和。而固有模态函数有两个约束条件：

1）在整个数据段内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个。
2）在任意时刻，由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零，即上、下包络线相对于时间轴局部对称。

啥意思？

用不严谨的语言和灵魂画师来解释一下：

1）图线要反复跨越x轴，像这样：

在整个数据段内，极值点的个数和过零点的个数必须相等或相差最多不能超过一个

而不能像这样某次穿过零点后出现多个极点：

极点数目偏多

2）包络线要对称，像这样：

包络线对称

而不能像这样：

包络线不对称

洗洗眼睛，看个正常点的例子吧：

EMD分解

上图由7张图片组成，其中第1张为原始信号，后边依次为EMD分解之后得到的6个分量，分别叫做IMF1~IMF5，最后一张图为残差，每一个IMF分量代表了原始信号中存在的一种固有模态函数。可以看出，每个IMF分量都是满足这两个约束条件的。

3 EMD分解步骤

EMD的分解过程是简单直观的：

1）根据原始信号上下极值点，分别画出上、下包络线。

上、下包络线

2）求上、下包络线的均值，画出均值包络线。

均值包络线

3）原始信号减均值包络线，得到中间信号。

原始信号减均值包络线

4）判断该中间信号是否满足IMF的两个条件，如果满足，该信号就是一个IMF分量；如果不是，以该信号为基础，重新做1）~4）的分析。IMF分量的获取通常需要若干次的迭代。

不满足约束2，需要继续迭代使用上述方法得到第一个IMF后，用原始信号减IMF1，作为新的原始信号，再通过1）~4）的分析，可以得到IMF2，以此类推，完成EMD分解。

迭代分解结果

4 IMF的物理含义

IMF1、IMF2这些分量的含义。这里举个简单的例子解释一下：

首先我们生成一段信号，它是由4Hz的正弦波、10Hz的正弦波和白噪声叠加而成的。如下图：

左侧3个信号叠加生成右侧信号

现在我们将合成后的信号进行EMD分解，结果如下图：

5张图分别为：原始信号、IMF1、IMF2、IMF3、残差

从EMD分解图中可以看出，IMF1为10Hz的分量，IMF2为4Hz的分量。所以简单（且理想化）地说，IMF的各个分量分别代表了原始信号中的各频率分量，并按照从高频到低频的顺序依次排列。这也就是（非常简单的情况下的）IMF的物理含义。

那IMF3等分量又代表了什么呢？在这个例子里，它们没有意义，属于EMD端点效应等带来的副作用，至于端点效应，在后续的文章里会讲到。

通过这个例子的说明，我们可以得到EMD的一大特点：自适应地进行信号主要成分分析（不是PCA）。也就是说，EMD分解信号不需要事先预定或强制给定基函数，而是依赖信号本身特征而自适应地进行分解。

当然，现实中的信号分量（IMF）不会像例子中一样保持完全稳定的频率和振幅，也常常无法从各分量中直接看出信号规律。EMD分解经常被用作信号特征提取的一个预先处理手段，将各IMF分量作为后续分析方法的输入，以完成更加复杂的工作。

5 示例代码

这个示例代码先生成4Hz的正弦波、10Hz的正弦波和白噪声叠加的模拟信号，然后进行EMD分解得到4个IMF分量，然后为每个IMF绘制频谱图。

import numpy as npimport matplotlib.pyplot as pltfrom PyEMD import EMD, Visualisation
# 生成一个示例信号，可以是两不同频率的信号叠加并加入高斯噪声t = np.linspace(0, 2, 100)  # 采样频率为50Hzsignal1 = 1*np.sin(2 * np.pi * 4 * t)  # 第一个信号，频率为4 Hzsignal2 = 3*np.sin(2 * np.pi * 10 * t)  # 第二个信号，频率为10 Hznoise = 0.2 * np.random.randn(len(t))  # 添加高斯噪声
signal = signal1 + signal2 + noise
# 创建EMD对象emd = EMD()
# 执行EMD分解imfs = emd(signal)
# 绘制原始信号和分解后的各个IMFsn_imfs = len(imfs)
plt.figure(figsize=(10, 6))plt.subplot(n_imfs + 1, 1, 1)plt.plot(t, signal, 'r')plt.title("Original Signal")
for i in range(n_imfs):    plt.subplot(n_imfs + 1, 1, i + 2)    plt.plot(t, imfs[i], 'b')    plt.title(f"IMF {i + 1}")
plt.tight_layout()
# 绘制每个IMF的频谱图plt.figure(figsize=(10, 6))plt.title("Frequency Spectra of IMFs")for i in range(n_imfs):    plt.subplot(n_imfs, 1, i + 1)    freq_spectrum = np.fft.fft(imfs[i])    freq_spectrum = np.abs(freq_spectrum)    freq = np.fft.fftfreq(len(t), d=t[1] - t[0])    plt.plot(freq[0:int(len(freq)/2)], freq_spectrum[0:int(len(freq_spectrum)/2)], 'b')    plt.title(f"IMF {i + 1} Frequency Spectrum")
plt.tight_layout()plt.show()

确保你已经安装了PyEMD库，可以使用以下命令来安装：

pip install EMD-signal -i https://pypi.tuna.tsinghua.edu.cn/simple pytest

输出结果：

从结果可以看出，IMF1为10Hz，IMF2为4Hz，证实了EMD可以自适应分解信号，但是需要注意的是这个只适合简单的信号。

编辑：张泽明

校核：李正平、张勇、王畅、赵栓栓、陈凯歌

该文为转载，已取得原作者许可，若侵权，后台联系小编进行删除。

来源：故障诊断与python学习

著作权归作者所有，欢迎分享，未经许可，不得转载

首次发布时间：2023-11-23