EMD和VMD学习记录

EMD

来源：浅谈经验模态分解-EMD | MuYi's Blog (muyi110.github.io)、EMD经验模态分解介绍_from pyemd import emd-CSDN博客

经验模态分解(Empirical Mode Decomposition, EMD)是依据数据自身的时间尺度特征来进行信号分解，无需预先设定任何基函数，是一种时频域信号处理方式。EMD在处理非平稳及非线性数据上具有明显的优势，适合分析非线性非平稳的信号序列，具有较高的信噪比。 平稳信号：分布参数（均值，方差，协方差等）或者分布规律随着时间不发生变化 非平稳信号：分布参数（均值，方差，协方差等）或者分布规律随着时间发生变化 该方法关键是经验模式分解，使复杂的信号分解为有限个本征模函数(Intrinsic Mode Function, IMF)，分解出来的各个IMF分量包含了原信号的不同时间尺度的局部特征信息。通俗理解就是，EMD是一台机器，将一堆混杂的硬币投进去，其会自动按照1角、5毛、1元等分好。

理论知识

原理

在物理上，如果瞬时频率有意义，那么函数必须是对称的，局部均值为零，并且具有相同的过零点和极值点数目。任何信号都是由若干本征模函数组成，一个本征模函数必须满足以下两个条件：

1. 函数在整个时间范围内，局部极值点和过零点的数目必须相等或最多相差一个

2. 在任意时刻点，局部最大值的包络（上包络线）和局部最小值的包络（下包络线） 平均必须为零

对于上述第二条说明：它把经典的全局性要求修改为局部性要求，使瞬时频率不再受不对称波形所形成的不必要的波动所影响, 实际上，这个条件应为“数据的局部均值是零”。但是对于非平稳数据来说，计算局部均值涉及到“局部时间尺度”的概念，而这是很难定义的。因此，在第二个条件中使用了局部极大值包络和局部极小值包络的平均为零来代替，使信号的波形局部对称

EMD将输入信号分解为几个本征模函数和一个残差组成，即由下列公式组成

$$
I(n)=\sum_{m=1}^M \operatorname{IMF}_m(n)+\operatorname{Res}_M(n)
$$

其中 I(n) 表示输入信号， \operatorname{IMF}_m(n) 表示 \mathrm{M}^{t h} 的本征模函数， \operatorname{Res}_M(n) 表示残差。

基本步骤

见解

与其他方法对比

变换类型	傅里叶 (Fourier)	小波 (Wavelet)	希尔伯特 (Hilbert)
基础	先验的	先验的	自适应的
频率	卷积：全局，不确定性	卷积：区域，不确定性	微分：局部，确定性
表示	能量-频率	能量-时间-频率	能量-时间-频率
非线性	否	否	是
非平稳性	否	是	是
特征提取	否	离散：否，连续：是	是
理论基础	理论完整	理论完整	经验的

无论是傅里叶变换还是小波分解，他们分解的信号都是从-inf~inf,这一点在实际中无法实现。 经验模态分解没有基函数，小波和傅里叶分解都有基函数，从这个角度来讲，经验模态分解是一种自适应分解方法。导致数据分解分解出来的若干个信号非常有限。 正常情况下，一般分解最多可以到15个IMF,正常情况下是分解到6个IMF。 经验模态分解没有严格的数学理论作为基础支撑。

缺点

EMD的主要缺点有如下

• 末端效应

  末端效应发生在信号的开始和结尾，因为在在信号开始之前和结尾之后没有样本点被考虑。大多数情况下，末端点并不是信号的极值，但是在执行EMD的过程中，极值包络线会在末断点发散(diverge)，导致错误，进而扭曲了IMF在末端点的波形，而且这种错误在EMD分解过程中会累积。

• 模态混叠问题

  模态混叠问题发生在EMD的执行过程中。出现下列情况之一就称为模态混合

  1. 在同一个IMF分量中,存在尺度分布范围很宽却又各不相同的信号

  2. 在不同的IMF分量中,存在着尺度相近的信号

  模态混叠问题使得特征提取、模型训练、模式识别变得困难，IMF失去了单一特征尺度的特征。集成经验模态分解(Ensemble empirical mode decomposition, EEMD)被提出用来解决模态混叠问题

改进方法

EEMD

为了改善测量的准确性，集成平均是有效的方法（也就是多次测量取平均值）。信号极值点影响IMF,若分布不均匀时会出现模态混叠，白噪声的频谱均匀分布，白噪声使得信号会自动分布到合适的参考尺度上。由于零均值噪声的特性，噪音经过多次的平均计算后会相互抵消，这样集成均值的计算结果就可以直接视作最终结果。集成均值的计算结果与原始信号的差值随着集成平均的次数增加而减少

EMD会产生模态混叠现象，即一个MF的信号频带会很宽，涵盖了低频和高频信号，实际中，我们希望某一个MF信号的频带尽可能集中。 EEMD(ensemble empirical mode decomposition)由于白噪声的尺度特征在整个时频域上均匀分布，进行EMD分解时，各种尺度特征中都会含有白噪声的成份。如果给信号添加白噪声后进行分解，那么信号含有的某种尺度特征，会被分解到与白噪声相对应的尺度特征上去。虽然噪声的不确定性，使得其在单次分解时，可能影响到信号的分解结果，但根据白噪声的统计特性，如果添加不同的白噪声，对大量分解结果求均值，各尺度特征中的噪声成份会相互抵消，最后只剩下待分解信号。

1. 给原始信号添加白噪声序列

2. 分解带有白噪声序列的输入信号，得到IMFs（第一次测量得到一系列的IMFs）

3. 重复第一和第二步，每次添加不同的白噪声序列（执行多次测量）

4. 获取相关IMFs集成的均值作为为最后的结果（理解为多次测量取平均值

添加的白噪声的两个属性

• 添加的白噪声导致所有时间尺度上极值分布的相对均匀分布(The added white noise leads to relatively even distribution of extrema distribution on all timescales)

• 通过集成平均，添加的噪声会被移除（噪声的均值为0）

CCEMD

CCEMD（complementary ensemble empirical mode decomposition），添加的高斯白噪声都是互补的，即满足：w^i[n]+w^{i+1}[n]=0​，一对一对的对噪声进行添加。

VMD

变分模态分解（VMD）简介 |作者：U Vamsi Krishna |中等的 --- Introduction to Variational Mode Decomposition(VMD) | by U Vamsi Krishna | Medium

变分模态分解（VMD）_variational mode decomposition-CSDN博客

浅谈：变分模态分解VMD | Feng Zhiheng's Blog

理论知识

原理

作者认为信号都可以分解为本征模态函数的叠加。

$$
u_k(t)=A_k(t) \cos \left(\phi_k(t)\right)
$$

基本步骤

1. 构造变分问题

   假设原始信号f被分解为k个分量，保证分解序列为具有中心频率的有限带宽的模态分量，同时各模态的估计带宽之和最小，约束条件为所有模态之和与原始信号 相等，则相应约束变分表达式为

   

2. 求解式（1），引入Lagrange乘法算子λ，将约束变分问题转变为非约束变分问题，得到增广Lagrange表达式为

   $$
   \begin{aligned} & L\left(\left\{u_k\right\},\left\{\omega_k\right\}, \lambda\right)= \\ & \alpha \sum_k\left\|\partial_t\left[(\delta(t)+j / \pi t) * u_k(t)\right] e^{-j \omega_k t}\right\|_2^2+ \\ & \left\|f(t)-\sum_k u_k(t)\right\|_2^2+\left\langle\lambda(t), f(t)-\sum_k u_k(t)\right\rangle \end{aligned}
   $$

    

   式中： \alpha 为二次征罚因子，作用是降低高斯噪声的干扰。

3. 利用交替方向乘子 (ADMM) 迭代算法结合Parseval/Plancherel、傅里叶等距变换，优化得到各模态分量和中心频率，并搜寻增广Lagrange函数的鞍点，交替寻优迭代后的 u_k ， \omega_k 和 \lambda 的表达式如下，详细过程:

   $$
   \begin{gathered} \widehat{u}_k^{n+1}(\omega) \leftarrow \frac{\widehat{f}(\omega)-\sum_{i / k} \widehat{u_i}(\omega)+\widehat{\lambda}(\omega) / 2}{1+2 \alpha\left(\omega-\omega_k\right)^2} \\ \omega_k^{n+1} \leftarrow \frac{\int_0^{\infty} \omega\left|\hat{u}_k^{n+1}(\omega)\right|^2 d \omega}{\int_0^{\infty}\left|\hat{u}_k^{n+1}(\omega)\right|^2 d \omega} \\ \widehat{\lambda}^{n+1}(\omega) \leftarrow \widehat{\lambda}^n(\omega)+\gamma\left(\widehat{f}(\omega)-\sum_k \widehat{u}_k^{n+1}(\omega)\right) \end{gathered}
   $$

    

   式中: \gamma 为噪声容忍度，满足信号分解的保真度要求， \widehat{u}_k^{n+1}(\omega) 、 \widehat{u} i(\omega) 、 \widehat{f}(\omega) 和 \widehat{\lambda}(\omega) 分别对应 u_k^{n+1}(t) 、 u_i(t) 、 f(t) 和 \lambda(t)​ 的傅里叶变换。

4. VMD主要迭代求解过程如下：

   S1: 初始化 \widehat{u}_k^1 、 \omega_k^1 、 \lambda^1 和最大迭代次数 N, n \leftarrow 0 ； S2: 利用公式 (3) 和（4）更新 \widehat{u}_k 和 \omega_k ； S3: 利用公式 (5) 更新 \widehat{\lambda}; S4: 精度收玫判据 \varepsilon>0 ，若不满足 \sum_k\left\|\widehat{u}_k^{n+1}-\widehat{u}_k^n\right\|_2^2 /\left\|\widehat{u}_k^n\right\|_2^2<\varepsilon 且 n<N ，则返回第二步，否则完成迭代，输出最终的 \widehat{u}_k 和 \omega_k​ 。

见解

浅谈：变分模态分解VMD | Feng Zhiheng's Blog

变分模态分解其实是经验模态分解EMD的一种扩展。

虽然EMD能够结合信号本身的特点分解出多个模态函数IMF，但是它也存在一些缺点：

首先：EMD对噪声敏感，怎么理解呢？在进行EMD分解时，我们需要利用信号的极大值和极小值来计算包络，而信号的极值很有可能收到噪声干扰，因此EMD对噪声比较敏感。

其次，EMD分解存在端点效应，什么是端点效应呢？端点效应指的是在信号的端点处，由于没有前后帧的信息，求得的包络面可能不准确，在后期分解的过程中，可能会得到不准确的IMF分量。

最后，EMD是经验模态分解，它很多时候依赖的是经验，数学基础还不完善。

EMD和VMD python实现

• EMD

import numpy as np
from PyEMD import  EMD
import pylab as plt
t = np.linspace(0,1,200)
s = np.cos(11*2*np.pi*np.pi*t*t) + 6*t*t
#EMD分解
IMF = EMD().emd(s,t)
N = IMF.shape[0]+1
#画图
plt.subplot(N,1,1)
plt.plot(t,s,"r")
plt.title("Input signal: $S(t)=cos(22/pi t^2 + 6t^2)$")
for n,imf in enumerate(IMF):
    plt.subplot(N,1,n+2)
    plt.plot(t,imf,"g")
    plt.title("IMF" + str(n+1))
plt.tight_layout()
plt.show()

• EEMD

import numpy as np
from PyEMD import  EEMD
import pylab as plt
t = np.linspace(0,1,200)
s = np.cos(11*2*np.pi*np.pi*t*t) + 6*t*t
#EEMD分解
eemd = EEMD()
emd = eemd.EMD
emd.extrema_detection = "parabol"
IMF = eemd.eemd(s,t)
N = IMF.shape[0]+1
#画图
plt.subplot(N,1,1)
plt.plot(t,s,"r")
plt.title("Input signal: $S(t)=cos(22/pi t^2 + 6t^2)$")
for n,imf in enumerate(IMF):
    plt.subplot(N,1,n+2)
    plt.plot(t,imf,"g")
    plt.title("IMF" + str(n+1))
plt.tight_layout()
plt.show()

• VMD

#. some sample parameters for VMD  
alpha = 5000      # moderate bandwidth constraint  
tau = 0           # noise-tolerance (no strict fidelity enforcement)  
K = 3              # 3 modes  
DC = 0             # no DC part imposed  
init = 1           # initialize omegas uniformly  
tol = 1e-7
​
#. Run actual VMD code  
u, u_hat, omega = VMD(v, alpha, tau, K, DC, init, tol)