信号处理基础之噪声与降噪(三) | EMD降噪与VMD降噪及python代码实现
接上期信号处理之噪声与降噪(二)|时域降噪方法(平滑降噪、SVD降噪)python代码实现，本期为大家介绍EMD降噪和VMD降噪，并给出python代码。后续将会介绍EMD族和小波族的降噪应用，敬请关注。
目录
1 EMD降噪
1.1 EMD的基本原理
1.2 EMD降噪的实现过程
1.3 EMD的不足
1.4 EMD降噪的python实现
2.1 VMD的基本原理
2.2 VMD降噪的实现过程
2.3 VMD的优点
2.4 VMD降噪的python实现
3.1 测试用例
3.2 降噪结果
3.3 降噪指标
1 EMD降噪
1.1 EMD的基本原理
经验模态分解（Empirical Mode Decomposition, EMD）是将复杂信号分解为有限个固有模态函数（Intrinsic Mode Function, IMF）和残差信号的方法，如图1所示。
图1 EMD与IMF
1.2 EMD降噪的实现过程
根据EMD分解过程不难看出，在进行信号处理时，只需要对信号分解的IMFs和残差信号进行约束，然后采用分解得到的IMFs对信号重构，即可获得降噪的信号。基于EMD的降噪过程包括以下步骤：
（1）EMD分解：得到一系列IMFs和残差信号；
（2）IMFs筛选：筛选包含有用信息的IMFs；
（3）信号重构：将筛选的IMFs叠加，获得降噪后的信号。
不难发现，基于EMD的降噪过程中，最关键的是IMFs的筛选，目前常用的方法包括频谱筛选法、峭度/峰值因子筛选法、相关系数筛选法和自适应法（该部分暂不进行展开）等。
频谱筛选法：计算各IMF的频谱，根据需要选取频段筛选IMFs，但这种方法基于研究人员清晰知道自己所需的频段；
峭度/峰值因子筛选法：计算各IMF的峭度值或峰值因子，一般认为峭度值越大信号中包含较多的毛刺，该类信号通常与噪声相关，重构时进行去除；
模态混叠：在同一个IMF分量中,存在尺度分布范围很宽却又各不相同的信号；在不同的IMF分量中,存在着尺度相近的信号。模态混叠主要是因为：在求取包络线的过程中，信号中存在异常事件，影响极值点的选取，使极值点分布不均匀，从而导致求取的包络为异常事件的局部包络和真实信号包络的组合。经该包络计算出的均值，再筛选出的IMF分量就包含了信号的固有模式和异常事件或者包含了相邻特征时间尺度的固有模式，从而产生了模态混叠现象。
from PyEMD import EMDimport numpy as npfrom scipy.stats import kurtosisdef denoised_with_EMD(data):imf_select = []order_select = []# EMD分解emd = EMD()imfs = emd.emd(data)for i in range(len(imfs)-1):print(kurtosis(imfs[i]))if kurtosis(imfs[i]) < 0:imf_select.append(imfs[i])order_select.append(i+1)else:continuesignal_denoised = np.sum(imf_select[:len(imf_select)], axis=0)return imf_select, order_select, signal_denoised
2 VMD降噪
变分模态分解（Variational Mode Decomposition, VMD）假设信号由一系列具有特定中心频率、有限带宽的子信号（IMFs）组成，如图2所示。因此VMD问题最终即是求解中心频率与带宽限制的问题，其在经典wiener滤波的基础上，通过变分问题进行求解，最终得到模态函数。VMD的数学原理和推导过程较为复杂，本文不进行详细展开，具体过程可参考文献[1]。
图2 VMD与IMF
实际上，变分问题最终是求解泛函的极值问题，最终可以转换成优化求解的问题，VMD过程可简单描述如下：
以下从优化模型的角度对变分问题进行简单介绍，假设信号可分解为K个有限带宽的模态分量
约束条件：模态和等于输入信号，表示如下：
以上优化模型的求解可采用诸多优化求解器，包括二次惩罚项、拉格朗日算子、增广拉格朗日函数等，各求解过程较为复杂，此处不展开讨论，感兴趣可参考文献[2]。
与EMD降噪过程类似，VMD降噪过程如下：
（1）VMD分解：得到各模态分量IMF；
（2）IMFs筛选：筛选包含有用信息的IMFs；
（1）可自定义模态个数；
（2）通过VMD方法分解出来的IMF都具有独立的中心频率，并且在频域上表现出稀疏性的特征，具备稀疏研究的特质；
（3）在对IMF求解过程中，通过镜像延拓的方式避免了类似EMD分解中出现的端点效应；
import numpy as npfrom scipy.stats import kurtosisfrom vmdpy import VMDdef denoised_with_VMD(data):imf_select = []order_select = []# VMD分解,alpha为惩罚因子，K为分解层数，tol为优化终止条件u, u_hat, omega = VMD(data, alpha=2000, tau=0.0, K=4, DC=0, init=1, tol=1e-7)for i in range(len(u)):if kurtosis(u[i]) < 0:imf_select.append(u[i])order_select.append(i+1)else:continuesignal_denoised = np.sum(imf_select[:len(imf_select)], axis=0)return imf_select, order_select, signal_denoised
3 算法测试
import matplotlib.pyplot as pltimport matplotlibn = 500 # 生成500个点的信号t = np.linspace(0, 30*np.pi, n, endpoint=False)s = np.cos(0.1*np.pi*t) + np.sin(0.3*np.pi*t) + np.cos(0.5*np.pi*t) + np.sin(0.7*np.pi*t) # 原始信号r = 0.5 * np.random.randn(n) # VMD噪声# r = 0.1 * np.random.randn(n) # EMD噪声y = s + r # 加噪声imf_select, order_select, signal_denoised = denoised_with_VMD(y)n = len(imf_select) + 2fig = plt.figure(figsize=(16, 12))plt.subplots_adjust(hspace=0.5)font = {'family': 'Times New Roman', 'size': 16, 'weight': 'normal',}matplotlib.rc('font', **font)plt.subplot(n, 1, 1)plt.plot(t, y, linewidth=1.5, color='b')plt.title('原始模拟信号', fontname='microsoft yahei', fontsize=16)for i in range(2,n):plt.subplot(n, 1, i)plt.plot(t, imf_select[i-2], linewidth=1.5, color='r')plt.title('IMF%d'%order_select[i-2], fontname='microsoft yahei', fontsize=16)plt.subplot(n, 1, n)plt.plot(t, signal_denoised, linewidth=1.5, color='b')plt.title('降噪信号', fontname='microsoft yahei', fontsize=16)
图3 EMD降噪
图4 VMD降噪
4 参考文献
[1]Dominique, Zosso, Konstantin, et al. Variational Mode Decomposition[J]. IEEE Transactions on Signal Processing A Publication of the IEEE Signal Processing Society.
[2]蔺梦雄,焦博森,杨琰,等. 基于VMD降噪的RV减速器故障诊断[J]. 北京联合大学学报,2023,37(3):7-13.
[3]https://zhuanlan.zhihu.com/p/608044525?utm_id=0
