本文提出了一種改進型的經驗模態分解算法用于心音圖(PCG)信號去噪,結合PCG的規則平均Shannon能量包絡算法,可有效提取PCG中的S1/S2成分。首先,通過小波變換和經驗模態分解結合算法對PCG信號進行濾波預處理;然后,提取預處理后PCG信號的固有模函數(IMF)時域、頻域特性及能量包絡;最后,結合信號的Shannon能量包絡和IMF相關特性準確定位出S1和S2。運用該方法對30例PCG信號進行測試,得到S1/S2成分的綜合識別率達99.75%。實驗結果表明,本文算法運用于S1/S2成分提取具有較好的效果,為進一步研究心音身份識別奠定基礎。
引用本文： 龔敬, 聶生東, 王遠軍. 基于改進經驗模態分解與能量包絡的S1/S2提取. 生物醫學工程學雜志, 2015, 32(5): 970-974. doi: 10.7507/1001-5515.20150173 復制
引言
在胸壁一定部位用聽診器聽取到由心肌收縮、心臟瓣膜關閉和血液撞擊心室壁、大動脈壁等引起的振動所產生的聲音，稱為心音。正常情況下，心音包含S1、S2、S3、S4四種成分，一般情況下，只能聽到S1和S2[1-5] 。S1主要由二尖瓣(M1)及三尖瓣(T1)產生，發生于心室收縮時，由于心室壓略大于心房壓導致房室瓣關閉的過程。S2是由于主動脈瓣及肺動脈瓣關閉所產生，代表心室舒張期開始，主要由主動脈瓣關閉音(A2)與肺動脈瓣關閉音(P2)組成。心音信號頻率主要集中在0～600 Hz，其中，S1的頻率成分主要集中于50～150 Hz，S2主要在50～200 Hz，250～350 Hz范圍內又出現第二個小波峰[6-7]。這兩種成分變化顯著，其他頻率范圍的信號則不太明顯，因此，在用心音信號做身份識別時，需要分離出S1和S2成分作為識別的重要特征。
文獻表明，目前對于心音信號的S1/S2成分分離研究，主要是對心音圖(phonocardiogram,PCG)信號進行預處理后，采用能量包絡的方法進行，而預處理方法往往決定后續對S1/S2成分提取的效果。在對PCG信號預處理方面，Meziani等[1]提出運用小波變換(wavelet transform)方法，成功分離出不同病理狀態下的心音信號；Hamza Cherif等[2]選用希爾伯特黃變換對PCG信號進行預處理，并提取出心音成分；Wang等[8]提出利用小波變換和能量包絡的方法提取出S1，準確率為93.24%；Safara等[9]提取小波包熵和貝葉斯分類器對心雜音進行分類，準確率達到96.94%。在算法上，Li等[10]利用經驗模態分解(empirical mode decomposition,EMD)與小波變換的結合算法EMD-Wavelet進行心電圖信號去噪，得到較好的效果。目前，相關參考文獻顯示希爾伯特黃變換和小波變換都已經成功運用于PCG信號處理，然而二者結合算法以及與能量包絡聯合提取S1/S2成分的研究較少。因此，本文提出利用EMD-Wavelet算法結合能量包絡方法，對S1和S2成分進行提取。
1 心音成分提取算法
1.1 心音預處理
通常情況下，經體表傳感器采集到的PCG信號不可避免地含有較多的噪聲，雖然硬件采集系統日趨完善，已具有較好的去噪功能，但僅依靠硬件設備不能完全解決PCG信號的噪聲干擾問題，需要采用數字濾波技術進行預處理。目前，數字濾波技術預處理去除噪聲的方法有很多，根據PCG信號的時頻特性，采用線性和非線性的數字濾波較為常見，如小波變換等方法已經成功被運用于PCG信號去噪，然而不同方法的去噪效果各有利弊[4, 11-13]。本文研究采用EMD和小波變換結合的方法，濾去PCG信號中的干擾噪聲。
1.1.1 經驗模態分解算法
EMD算法作為一種自適應、高效的信號分解方法，可以將非線性不平穩信號從高頻到低頻分解為若干個固有模函數(intrinsic mode function，IMF)的線性疊加，且每個IMF的頻率分辨率也隨原始信號的改變而改變。由于該分解以局部時間尺度為基礎，因此適用于非線性、非平穩信號。其中，IMF滿足以下兩個條件：①在全部數據區間中，極值點數目與過零點數目相等或最多相差1個；②在任一點處，局部極大值點定義的包絡和局部極小值點定義的包絡的均值為零。
對于信號x(t)，其EMD分解過程如下[10, 14]：首先，求得信號x(t)所有的極大值點和極小值點，并選用三次樣條曲線進行插值擬合，從而獲得數據的上、下包絡線。計算上、下包絡線的均值m1，以及x(t)與m1的差值h1=x(t)-m1。隨后，判斷h1是否為IMF，若不是則對其重復上述計算過程，直到第n次的h1n為IMF，令c1=h1n,則稱c1是從x(t)中提取的第1個IMF，包含著信號的最高頻率成分。然后，求原信號x(t)與c1的差值R1=x(t)-
為保證獲得的IMF分量振幅和頻率有充分的物理意義，本文選取標準差作為終止準則，終止條件SD一般采用如下公式定義：
| $ x\left( t \right) = \sum\limits_{i = 1}^n {{c_n}\left( t \right) + {R_n}} $ |
其中，α表示終止系數，通常設置為0.2～0.3，本文選取0.25。
1.1.2 小波變換去噪
選用小波變換的方法進行去噪，具體可分為以下三個步驟[15]：①選擇合適的小波函數將原始信號進行m層分解，獲取各個尺度上的細節和近似分量。通常情況下，噪聲主要包含在高頻部分，即細節系數中。②選取合適的閾值對各個尺度上的細節分量進行閾值量化處理，獲取新的小波系數。③利用第n層近似分量和經閾值量化處理的各層細節分量進行重構，獲得去噪后信號。
1.1.3 EMD-Wavelet預處理算法
利用EMD能將信號分解到不同頻率范圍的特性，用小波變換的方法對特定頻帶范圍內的IMF進行預處理，最終得到去噪后的心音信號。選取進行去噪處理的IMF主要參照標準為頻帶范圍。流程圖如圖 1所示，具體步驟如下：
圖1
改進EMD-Wavelet算法去噪流程圖
Figure1.
EMD-Wavelet algorithm flow chart
(1)對原始PCG信號進行EMD，得到n層IMF分量；
(2)對n層IMF進行FFT變換，計算每層的頻率范圍；
(3)觀察n層IMF的頻譜圖，判斷其是否為高頻，將主要頻帶分布在50 Hz以上的IMF保存得到IMF2′，剩余IMF保留為IMF1′；
(4)利用小波變換對IMF2′分別進行去噪處理，并將所得各IMF2′分量與IMF1′結合重構PCG信號。
利用EMD對正常PCG信號進行分解得到各層IMF分量及其每層IMF的頻譜圖，利用離散小波變換對其進行去噪處理。通過對比試驗本文選用db7小波基函數進行IMF分量的去噪處理，分解尺度為5[16]。最終結合其他IMF分量重建信號，得到去噪后的PCG信號，如圖 2所示。
圖2
心音信號去噪效果圖
Figure2.
Denoised PCG signal
1.2 S1/S2成分提取
對信號s(i)的規則平均Shannon能量變換主要包括以下步驟[17]：
(1)對信號s(i)進行歸一化處理，公式如下：
| $ {s_{norm}}\left( i \right) = \frac{{s\left( i \right)}}{{\max \left( {\left| {s\left( i \right)} \right|} \right)}} $ |
其中max(｜s(i)｜)表示s(i)絕對值的最大值；
(2)計算平均Shannon能量，計算每連續間隔20 ms的s(i)平均Shannon能量，計算時有10 ms和前一次樣本點計算重復，公式如下：
| $ ES\left( {{n_i}} \right) = - \frac{1}{n}\sum\limits_{j - 1}^m {{s^2}\left( i \right)\log {s^2}\left( i \right)} $ |
其中m表示20 ms信號s(i)的長度；
(3)計算全局信號的規則平均Shannon能量，通過計算20 ms信號的平均Shannon能量后，計算其規則平均Shannon能量。具體公式如下：
| $ {p_e}\left( i \right) = \frac{{ES\left( {{n_i}} \right) - mean\left( {ES\left( {{n_i}} \right)} \right)}}{{{\rm{std}}\left( {ES\left( {{n_i}} \right)} \right)}} $ |
其中mean(ES(ni))為信號平均Shannon能量的平均值，std(ES(ni))表示信號平均Shannon能量標準差。
通過計算全局PCG信號的規則平均Shannon能量，提取出其Shannon能量包絡，可以發現在S1、S2處其能量幅值較大，且S1能量幅值大于S2。因而，可進一步識別并提取出S1和S2，具體步驟如下(如圖 3所示)：
圖3
S1、S2檢測流程圖
Figure3.
Flow chart of detecting S1/S2 components
(1)計算各層IMF與原始信號的互相關系數，找出IMF分量中與原始PCG信號相關系數最大的IMF分量IMFn；
(2)分別計算去噪后IMFn分量和去噪后PCG信號的規則平均Shannon能量；
(3)結合兩路規則平均Shannon能量，分別設定閾值判斷S1、S2的起止位置；
(4)綜合上述判斷，校正S1、S2起始位置，提取S1、S2成分。
我們可以從心音規則平均Shannon能量包絡圖中看出，S1、S2成分作為一系列正向波，其幅度較大；但是在包含較多噪聲的PCG信號的心音規則平均Shannon能量包絡圖，其噪聲表現為幅度較小的正向波。因此可通過設定合適的閾值將噪聲去除，本文選取的閾值為Thr=δ*max(ES(ni))，其中δ為參考系數，通過多次對比試驗取δ=5%,max(ES(ni))表示計算的規則平均Shannon能量的最大值[18]。檢測過程中，將大于閾值Thr的成分置為1，小于Thr的成分置為0。
此外，由于通常S1時限在70～150 ms，S2時限在60～120 ms，我們可以在檢測到的S1與S2峰值前后100 ms的時間范圍內尋找其起始點與終止點。同時，由于在PCG圖中有意義的心音成分的寬度一般都在30 ms以上，因此在心音規則平均Shannon能量包絡圖中，可以選用如下準則作為校正標準：若正向波之間的間隔寬度小于20 ms，則可識別為是同一個心音成分；寬度小于30 ms的正向波一般來說可能是噪聲引起的，故判定為噪聲成分。通過校正檢測，我們成功提取出S1、S2成分，如圖 4所示。
圖4
S1/S2成分提
由上到下分別是：去噪后心音信號、心音信號Shannon能量包絡、S1/S2成分標記
Figure4. Results of extracting S1/S2 componentsEach figure from top to bottom: denoised PCG signal, PCG signal’s Shannon energy, location of S1/S2 components
2 結果與討論
由于S1和S2的頻率主要分布于50～200 Hz，且噪聲一般分布在高頻分量中，因而在運用EMD-Wavelet算法進行去噪時，需要篩選出50 Hz以上的高頻IMF分量，并對這些高頻IMF分量進行小波去噪。同時，在進行EMD分解的過程中終止準則的限定也十分重要，一般選用標準差作為限定條件，通常取值為0.2～0.3，本文選取0.25。
在定位S1和S2成分時，選用IMF分量和原始信號的相關系數作為參考，主要是為了找出原始信號的主要成分分布，并以此作為初步定位的條件和標準具有一定的參考性。在最終提取S1和S2成分時，本文在運用IMF分量初步判定的基礎上結合全局的Shannon能量包絡，能夠比較準確地定位出S1和S2成分。
本文選取eGeneralMedical.com心音數據庫中30例PCG信號進行試驗，并與文獻[5]中方法對比，結果如表 1所示，準確率達到99.75%，證明了利用改進EMD和能量包絡提取S1、S2信號方法的可行性與準確性，為進一步研究PCG信號用于身份識別奠定基礎。
表1
S1、S2信號提取結果
Table1.
Extraction of S1 and S2 components
運用改進EMD-Wavelet算法對PCG信號進行預處理，具有一定的優勢，主要表現在該算法可以明確得到信號分解所在的頻帶范圍，可以有針對性地對信號含噪聲的高頻成分進行處理，將該算法與Shannon能量包絡結合能夠準確定位S1和S2成分，避免一次判定對S1和S2成分的丟失，從而比較準確地提取出S1和S2成分。
引言
在胸壁一定部位用聽診器聽取到由心肌收縮、心臟瓣膜關閉和血液撞擊心室壁、大動脈壁等引起的振動所產生的聲音，稱為心音。正常情況下，心音包含S1、S2、S3、S4四種成分，一般情況下，只能聽到S1和S2[1-5] 。S1主要由二尖瓣(M1)及三尖瓣(T1)產生，發生于心室收縮時，由于心室壓略大于心房壓導致房室瓣關閉的過程。S2是由于主動脈瓣及肺動脈瓣關閉所產生，代表心室舒張期開始，主要由主動脈瓣關閉音(A2)與肺動脈瓣關閉音(P2)組成。心音信號頻率主要集中在0～600 Hz，其中，S1的頻率成分主要集中于50～150 Hz，S2主要在50～200 Hz，250～350 Hz范圍內又出現第二個小波峰[6-7]。這兩種成分變化顯著，其他頻率范圍的信號則不太明顯，因此，在用心音信號做身份識別時，需要分離出S1和S2成分作為識別的重要特征。
文獻表明，目前對于心音信號的S1/S2成分分離研究，主要是對心音圖(phonocardiogram,PCG)信號進行預處理后，采用能量包絡的方法進行，而預處理方法往往決定后續對S1/S2成分提取的效果。在對PCG信號預處理方面，Meziani等[1]提出運用小波變換(wavelet transform)方法，成功分離出不同病理狀態下的心音信號；Hamza Cherif等[2]選用希爾伯特黃變換對PCG信號進行預處理，并提取出心音成分；Wang等[8]提出利用小波變換和能量包絡的方法提取出S1，準確率為93.24%；Safara等[9]提取小波包熵和貝葉斯分類器對心雜音進行分類，準確率達到96.94%。在算法上，Li等[10]利用經驗模態分解(empirical mode decomposition,EMD)與小波變換的結合算法EMD-Wavelet進行心電圖信號去噪，得到較好的效果。目前，相關參考文獻顯示希爾伯特黃變換和小波變換都已經成功運用于PCG信號處理，然而二者結合算法以及與能量包絡聯合提取S1/S2成分的研究較少。因此，本文提出利用EMD-Wavelet算法結合能量包絡方法，對S1和S2成分進行提取。
1 心音成分提取算法
1.1 心音預處理
通常情況下，經體表傳感器采集到的PCG信號不可避免地含有較多的噪聲，雖然硬件采集系統日趨完善，已具有較好的去噪功能，但僅依靠硬件設備不能完全解決PCG信號的噪聲干擾問題，需要采用數字濾波技術進行預處理。目前，數字濾波技術預處理去除噪聲的方法有很多，根據PCG信號的時頻特性，采用線性和非線性的數字濾波較為常見，如小波變換等方法已經成功被運用于PCG信號去噪，然而不同方法的去噪效果各有利弊[4, 11-13]。本文研究采用EMD和小波變換結合的方法，濾去PCG信號中的干擾噪聲。
1.1.1 經驗模態分解算法
EMD算法作為一種自適應、高效的信號分解方法，可以將非線性不平穩信號從高頻到低頻分解為若干個固有模函數(intrinsic mode function，IMF)的線性疊加，且每個IMF的頻率分辨率也隨原始信號的改變而改變。由于該分解以局部時間尺度為基礎，因此適用于非線性、非平穩信號。其中，IMF滿足以下兩個條件：①在全部數據區間中，極值點數目與過零點數目相等或最多相差1個；②在任一點處，局部極大值點定義的包絡和局部極小值點定義的包絡的均值為零。
對于信號x(t)，其EMD分解過程如下[10, 14]：首先，求得信號x(t)所有的極大值點和極小值點，并選用三次樣條曲線進行插值擬合，從而獲得數據的上、下包絡線。計算上、下包絡線的均值m1，以及x(t)與m1的差值h1=x(t)-m1。隨后，判斷h1是否為IMF，若不是則對其重復上述計算過程，直到第n次的h1n為IMF，令c1=h1n,則稱c1是從x(t)中提取的第1個IMF，包含著信號的最高頻率成分。然后，求原信號x(t)與c1的差值R1=x(t)-
為保證獲得的IMF分量振幅和頻率有充分的物理意義，本文選取標準差作為終止準則，終止條件SD一般采用如下公式定義：
| $ x\left( t \right) = \sum\limits_{i = 1}^n {{c_n}\left( t \right) + {R_n}} $ |
其中，α表示終止系數，通常設置為0.2～0.3，本文選取0.25。
1.1.2 小波變換去噪
選用小波變換的方法進行去噪，具體可分為以下三個步驟[15]：①選擇合適的小波函數將原始信號進行m層分解，獲取各個尺度上的細節和近似分量。通常情況下，噪聲主要包含在高頻部分，即細節系數中。②選取合適的閾值對各個尺度上的細節分量進行閾值量化處理，獲取新的小波系數。③利用第n層近似分量和經閾值量化處理的各層細節分量進行重構，獲得去噪后信號。
1.1.3 EMD-Wavelet預處理算法
利用EMD能將信號分解到不同頻率范圍的特性，用小波變換的方法對特定頻帶范圍內的IMF進行預處理，最終得到去噪后的心音信號。選取進行去噪處理的IMF主要參照標準為頻帶范圍。流程圖如圖 1所示，具體步驟如下：
圖1
改進EMD-Wavelet算法去噪流程圖
Figure1.
EMD-Wavelet algorithm flow chart
(1)對原始PCG信號進行EMD，得到n層IMF分量；
(2)對n層IMF進行FFT變換，計算每層的頻率范圍；
(3)觀察n層IMF的頻譜圖，判斷其是否為高頻，將主要頻帶分布在50 Hz以上的IMF保存得到IMF2′，剩余IMF保留為IMF1′；
(4)利用小波變換對IMF2′分別進行去噪處理，并將所得各IMF2′分量與IMF1′結合重構PCG信號。
利用EMD對正常PCG信號進行分解得到各層IMF分量及其每層IMF的頻譜圖，利用離散小波變換對其進行去噪處理。通過對比試驗本文選用db7小波基函數進行IMF分量的去噪處理，分解尺度為5[16]。最終結合其他IMF分量重建信號，得到去噪后的PCG信號，如圖 2所示。
圖2
心音信號去噪效果圖
Figure2.
Denoised PCG signal
1.2 S1/S2成分提取
對信號s(i)的規則平均Shannon能量變換主要包括以下步驟[17]：
(1)對信號s(i)進行歸一化處理，公式如下：
| $ {s_{norm}}\left( i \right) = \frac{{s\left( i \right)}}{{\max \left( {\left| {s\left( i \right)} \right|} \right)}} $ |
其中max(｜s(i)｜)表示s(i)絕對值的最大值；
(2)計算平均Shannon能量，計算每連續間隔20 ms的s(i)平均Shannon能量，計算時有10 ms和前一次樣本點計算重復，公式如下：
| $ ES\left( {{n_i}} \right) = - \frac{1}{n}\sum\limits_{j - 1}^m {{s^2}\left( i \right)\log {s^2}\left( i \right)} $ |
其中m表示20 ms信號s(i)的長度；
(3)計算全局信號的規則平均Shannon能量，通過計算20 ms信號的平均Shannon能量后，計算其規則平均Shannon能量。具體公式如下：
| $ {p_e}\left( i \right) = \frac{{ES\left( {{n_i}} \right) - mean\left( {ES\left( {{n_i}} \right)} \right)}}{{{\rm{std}}\left( {ES\left( {{n_i}} \right)} \right)}} $ |
其中mean(ES(ni))為信號平均Shannon能量的平均值，std(ES(ni))表示信號平均Shannon能量標準差。
通過計算全局PCG信號的規則平均Shannon能量，提取出其Shannon能量包絡，可以發現在S1、S2處其能量幅值較大，且S1能量幅值大于S2。因而，可進一步識別并提取出S1和S2，具體步驟如下(如圖 3所示)：
圖3
S1、S2檢測流程圖
Figure3.
Flow chart of detecting S1/S2 components
(1)計算各層IMF與原始信號的互相關系數，找出IMF分量中與原始PCG信號相關系數最大的IMF分量IMFn；
(2)分別計算去噪后IMFn分量和去噪后PCG信號的規則平均Shannon能量；
(3)結合兩路規則平均Shannon能量，分別設定閾值判斷S1、S2的起止位置；
(4)綜合上述判斷，校正S1、S2起始位置，提取S1、S2成分。
我們可以從心音規則平均Shannon能量包絡圖中看出，S1、S2成分作為一系列正向波，其幅度較大；但是在包含較多噪聲的PCG信號的心音規則平均Shannon能量包絡圖，其噪聲表現為幅度較小的正向波。因此可通過設定合適的閾值將噪聲去除，本文選取的閾值為Thr=δ*max(ES(ni))，其中δ為參考系數，通過多次對比試驗取δ=5%,max(ES(ni))表示計算的規則平均Shannon能量的最大值[18]。檢測過程中，將大于閾值Thr的成分置為1，小于Thr的成分置為0。
此外，由于通常S1時限在70～150 ms，S2時限在60～120 ms，我們可以在檢測到的S1與S2峰值前后100 ms的時間范圍內尋找其起始點與終止點。同時，由于在PCG圖中有意義的心音成分的寬度一般都在30 ms以上，因此在心音規則平均Shannon能量包絡圖中，可以選用如下準則作為校正標準：若正向波之間的間隔寬度小于20 ms，則可識別為是同一個心音成分；寬度小于30 ms的正向波一般來說可能是噪聲引起的，故判定為噪聲成分。通過校正檢測，我們成功提取出S1、S2成分，如圖 4所示。
圖4
S1/S2成分提
由上到下分別是：去噪后心音信號、心音信號Shannon能量包絡、S1/S2成分標記
Figure4. Results of extracting S1/S2 componentsEach figure from top to bottom: denoised PCG signal, PCG signal’s Shannon energy, location of S1/S2 components
2 結果與討論
由于S1和S2的頻率主要分布于50～200 Hz，且噪聲一般分布在高頻分量中，因而在運用EMD-Wavelet算法進行去噪時，需要篩選出50 Hz以上的高頻IMF分量，并對這些高頻IMF分量進行小波去噪。同時，在進行EMD分解的過程中終止準則的限定也十分重要，一般選用標準差作為限定條件，通常取值為0.2～0.3，本文選取0.25。
在定位S1和S2成分時，選用IMF分量和原始信號的相關系數作為參考，主要是為了找出原始信號的主要成分分布，并以此作為初步定位的條件和標準具有一定的參考性。在最終提取S1和S2成分時，本文在運用IMF分量初步判定的基礎上結合全局的Shannon能量包絡，能夠比較準確地定位出S1和S2成分。
本文選取eGeneralMedical.com心音數據庫中30例PCG信號進行試驗，并與文獻[5]中方法對比，結果如表 1所示，準確率達到99.75%，證明了利用改進EMD和能量包絡提取S1、S2信號方法的可行性與準確性，為進一步研究PCG信號用于身份識別奠定基礎。
表1
S1、S2信號提取結果
Table1.
Extraction of S1 and S2 components
運用改進EMD-Wavelet算法對PCG信號進行預處理，具有一定的優勢，主要表現在該算法可以明確得到信號分解所在的頻帶范圍，可以有針對性地對信號含噪聲的高頻成分進行處理，將該算法與Shannon能量包絡結合能夠準確定位S1和S2成分，避免一次判定對S1和S2成分的丟失，從而比較準確地提取出S1和S2成分。
