類EMD的「信號分解方法」及MATLAB實現（第二篇）——CEEMD

  • 2021 年 2 月 1 日
  • AI

縮寫為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);

上邊的測試代碼和封裝函數，包括工具箱都可以在公眾號（括號的城堡）中獲取，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.

VirMach 便宜 VPS

QNews

QNews