类EMD的“信号分解方法”及MATLAB实现（第二篇）——CEEMD

缩写为CEEMD的方法其实不止一种，包括互补集合经验模态分解方法[1]（Complementary Ensemble Empirical Mode Decomposition,2010）和完全集合经验模态分解方法[2]（Complete Ensemble Empirical Mode Decomposition,2011）。本文中所探讨的是上述第一种方法。

1. CEEMD（互补集合经验模态分解）的概念

上一篇我们介绍了EMD的一种最常见的衍生方法EEMD，这次要讲到的CEEMD是从EEMD方法进一步优化而来的，既然是优化那就必有所针对，CEEMD针对的就是EEMD的“残余辅助噪声”。

为什么会有残余辅助噪声呢？因为EEMD的前提是认为“多组白噪声的叠加近似等于0”。然而当处理的次数不够多的时候，白噪声往往不能被降低到忽略不计的程度。

反过来讲，如果使用EEMD方法时想要获得残余噪声较小的结果，就需要增加平均处理的次数，这样无疑会增加计算量。

为了解决这个问题，CEEMD的解决思路是：

CEEMD：将一对互为相反数的正负白噪声作为辅助噪声加入源信号当中，以消除原来 EEMD 方法分解后重构信号当中残留的多余的辅助白噪声，同时减少分解时所需的迭代次数，降低计算成本。 [3]

具体的方法可以说是非常简单直接了：与EEMD相比，CEEMD的区别仅仅在于添加白噪声的方式上。EEMD添加的是相互独立的白噪声；CEEMD添加的是成对的、互为相反数的白噪声序列。

为了对比残余噪声，我们分别计算使用EEMD和CEEMD方法对信号分解再重构之后的残余量[1]：

EEMD方法的重构残余量

CEEMD方法的重构残余量

可以看出CEEMD方法的残余辅助噪声比EEMD要低十几个数量级。Yeh展示了在某段信号下两种方法处理后的白噪声残余随叠加次数M的变化趋势（下图），EEMD方法要在将近10000次累加之后才能将残余量降到CEEMD方法的水平，而CEEMD则在个位数的处理次数下就能达到这个水平。

2. CEEMD的编程实现

ceemd函数没有出现在MATLAB的官方库中（截至MATLAB 2020b），不过这个方法的编程并不难，因为它的处理流程与eemd方法基本一致，网上也可以找到从eemd基础上改造而来的程序，不过笔者还是部分地修改了一些。修改后的ceemd代码与eemd相同，需要调试的参数有两个，分别是用于平均处理的次数M、添加的白噪声的幅值（白噪声的幅值通常用“白噪声幅值的标准差与原始信号幅值标准差之比”来表征）。

现在我们再来验证一下CEEMD方法相对于EEMD的优越性，还是采用上一篇文章相同的测试信号。

待分解的测试信号

优越性有两个方面：1.相同的累加次数下，CEEMD的白噪声残余更小。2.CEEMD使用更少的计算资源（即更小的累加次数）即可得到理想的分解结果。第一个优越性上一节已经演示过了，不再赘述。

在累加次数为8，白噪声幅值为0.2时，CEEMD分解的结果为：

CEEMD分解结果

CEEMD分解的IMF1、IMF2含有高频的正弦间歇性信号，IMF2可以看做IMF1很小的能量损失，分析高频信号时，可以将IMF1、2叠加起来作为重构的高频信号，会得到更好的分析效果。IMF4也很好地提取了信号中的低频分量。

而在相同的参数设置（M=8,nstd=0.2）下，EEMD的分解结果为：

EEMD分解结果

可以看到IMF1、IMF2、IMF3中，混入的噪声明显更大。如果想得到像CEEMD类似的处理结果，则累加次数要增加到100以上了。由此可见CEEMD方法可以大大节省计算资源。

上述中进行CEEMD分解的程序如下：

Nstd = 0.2; %Nstd为附加噪声标准差与Y标准差之比
NE = 8;   %NE为对信号的平均次数
imf = pCEEMD(sig,t,Nstd,NE);
% 画信号CEEMD分解图
% 输入：
% y为待分解信号
% FsOrT为采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量。如果未知采样频率，可设置为1
% Nstd为附加噪声标准差与Y标准差之比
% NE为对信号的平均次数
% 输出：
% imf为经CEEMD分解后的各imf分量值
% 例1：（FsOrT为采样频率）
% fs = 100;
% t = 1/fs:1/fs:1;
% data = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pCEEMD(data,fs,0.2,100);
% 例2：（FsOrT为时间向量，需要注意此时FsOrT的长度要与y相同）
% t = 0:0.01:1;
% data = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pCEEMD(data,t,0.2,100);

上述程序中的pCEEMD是笔者经过再次封装的ceemd画图程序。因为本专栏计划将“类EMD”方法统一起来，所以对诸如emd/eemd/ceemd以及后边可能会讲到的vmd等方法全部统一格式，以解决各代码来源不同（比如来自MATLAB官方库、第三方库和一些零散的程序）的问题。上边imf即为ceemd分解后的各分量信号，调用函数时CEEMD分解的图也可以画出来。

对于有些应用场景，还需要对各imf分量的频谱进行分析，就需要如下这样的图：

画这个图也同样封装成了一行代码就可以实现的形式：

imf = pCEEMDandFFT(sig,fs,Nstd,NE);% 画信号CEEMD分解与各IMF分量频谱对照图
% function imf = pCEEMDandFFT(y,FsOrT,Nstd,NE)
% 输入：
% y为待分解信号
% FsOrT为采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量
% Nstd为附加噪声标准差与Y标准差之比
% NE为对信号的平均次数
% 输出：
% imf为经CEEMD分解后的各imf分量值
% 例1：（FsOrT为采样频率）
% fs = 100;
% t = 1/fs:1/fs:1;
% y = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pCEEMDandFFT(y,fs,0.2,100);
% 例2：（FsOrT为时间向量，需要注意此时FsOrT的长度要与y相同）
% t = 0:0.01:1;
% y = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pCEEMDandFFT(y,t,0.2,100);

上边的测试代码和封装函数，包括工具箱都可以在下述链接中获取：

CEEMD画图工具（公开版） | 工具箱文档

EMD、EEMD以及HHT相关的程序也有，编程不易，感谢支持~关于EMD、EEMD和HHT的相关介绍可以看这里：

Mr.括号：这篇文章能让你明白经验模态分解（EMD）——EMD在MATLAB中的实现方法

Mr.括号：类EMD的“信号分解方法”及MATLAB实现（第一篇）——EEMD

Mr.括号：希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位——Hilbert分析衍生方法及MATLAB实现

3. 更多

后续还会逐渐补充CEEMDAN、MEEMD、VMD、LMD以及小波分解、小波包分解、SWT、EWT等等“信号分解方法”，把这一系列做的尽量全面一些。有其他想让博主补充的也可以在评论区留言，合适的话会一起加入该系列豪华大餐哦~

[1] Yeh J R , Shieh J S , Huang N E . Complementary Ensemble Empirical Mode Decomposition: a Novel Noise Enhanced Data Analysis Method.[J]. Advances in Adaptive Data Analysis, 2010, 2(2):-.

[2] Mar´ıa E. Torres ★, Marcelo A. Colominas ★, Gasto´n Schlotthauer ★, et al. A complete ensemble empirical mode decomposition with adaptive noise[C]// IEEE International Conference on Acoustics. IEEE, 2011.

[3] 汪一飞. 基于全矢CEEMD的滚动轴承故障诊断研究[D].郑州大学,2019.