針對基于小波分析和經驗模式分解的降噪方法本質上不能追蹤并消除噪聲且容易造成心音失真的問題,本文提出了一種融合改進最小值控制遞歸平均和最優修正對數譜幅度估計的心音降噪方法。該方法采用短時窗平滑動態追蹤、估計噪聲最小值,并將噪聲估計用于獲取最優頻譜增益函數,通過最小化干凈心音與估計的干凈心音的差異來最大限度地抑制噪聲。此外,結合心音時頻圖主觀分析和對正常與異常心音分類系統貢獻性的客觀分析,提出了一種更為嚴格的評價機制。實驗結果表明,本文提出的方法有效地改善了心音信號的時頻特征,并在正常與異常心音分類系統中獲得了更高的評分。本文提出方法能夠幫助醫務工作者提高聽診的準確性,對計算機輔助診斷系統的構建與應用也具有十分重要的參考價值。
引用本文: 許春冬, 周靜, 應冬文, 辛鵬麗. 噪聲動態估計下的心音降噪. 生物醫學工程學雜志, 2020, 37(5): 775-785. doi: 10.7507/1001-5515.202002023 復制
引言
心音信號分析是診斷心血管疾病的重要手段,被廣泛應用于臨床診斷和體檢當中[1-2]。心音信號的分析主要包括醫生經驗聽診和計算機輔助診斷[2-3],現今多借助于電子聽診器進行采集[4]。然而,在采集過程中,心音信號不可避免地會受到環境噪聲、采集設備噪聲以及人體其它器官噪聲等因素的干擾[3, 5],使得原本就十分弱小的心音信號的可分析性極大降低。因此,在對心音信號分析之前,需要先做降噪處理。
近幾十年來,研究和應用較廣的心音降噪方法當屬基于小波分析的方法[3, 6]和基于經驗模式分解(empirical mode decomposition,EMD)的方法[5]。其中,基于小波分析的方法是采用小波分解獲取小波域內不同尺度的小波系數,根據噪聲對應的小波系數小于心音對應的小波系數的假設,建立適當的閾值函數對各尺度的小波系數進行線性濾波處理,抑制噪聲的小波系數,達到降噪的目的[3, 6];此外,文獻[3, 5]也指出根據心音信號的頻帶特點,將部分高頻細節系數和低頻近似系數置零,能進一步去除不必要的噪聲。基于 EMD 的方法是將信號分解為多個模態固有函數(intrinsic mode function,IMF)分量,通過對模態混疊臨界點的判斷,將臨界點前的 IMF 分量視為噪聲并直接丟棄,同時丟棄心音成分少、噪聲含量大的高頻 IMF 分量[5, 7],并建立合適的濾波機制,將心音信號從保留的混疊 IMF 分量中分離出來,如互補總體 EMD(complementary ensemble EMD,CEEMD)[5]。然而,實際采集環境下的干擾噪聲具備隨機、非平穩、混沌等特性[8],且心音具備多樣性[9],小波系數的大小難以精準區分心音信號和噪聲(如類似心音的噪聲)[10],這使得獲取精準的自適應閾值函數更為困難,因此要想通過小波或 EMD 的空間變換和線性濾波處理將噪聲和心音完全分離開來十分困難,且容易出現心音失真情況[10-11]。同時,近些年的研究表明[4, 12-13],高頻部分也蘊含著有用的特征,尤其是在一些異常心音(如病變心雜音)中,故不可簡單濾波或者直接丟棄。此外,心音降噪還涉及到心雜音等由心臟異常引起的額外音,且許多額外音與噪聲特征相近,如三尖瓣關閉不全時在第一心音(S1)與第二心音(S2)間產生的心雜音[14-15]。作為病理特征的心雜音對心音分類有著重要的作用,不可直接濾除[13],這也進一步加大了心音降噪的難度。在心音降噪處理中,應以心音不失真為前提,最大程度地抑制噪聲[3, 5]。
實際中無法直接獲得完全干凈的心音信號,這就使得各類算法的降噪結果沒有標準的參考。目前常用的評價方法主要包括以相對干凈的心音信號為參考計算評價指標,如信噪比(signal-to-noise ratio,SNR)、均方根誤差(root mean square error,RMSE)[3, 5]、主觀譜圖分析[5]、平均意見得分(mean opinion score,MOS)[16]。其中,以相對干凈的心音信號為參考計算評價指標的分析并不嚴格,其客觀性并不強。
針對心音降噪面臨的問題,本文提出了一種融合改進最小值控制遞歸平均(improved minimum control recursive average,IMCRA)和最優修正對數譜幅度(optimally modified log-spectral amplitude,OMLSA)估計的心音降噪方法(后文簡稱:IMCRA-OMLSA)。該方法能夠從不同于小波分析和 EMD 分解降噪方法的角度,通過最優增益函數譜估計,最小化目標估計心音與干凈心音的差異,有效地抑制噪聲,降低心音失真的風險。在對算法的評價中,由于完全干凈的心音不可知,故提出譜圖主觀分析和降噪算法對分類系統貢獻度的客觀分析相結合的評價機制,確保對算法評價的客觀性及有效性。本文提出算法在心音噪聲中具備良好的性能,該算法對降噪聽診器的設計、心血管疾病計算機輔助診斷系統的構建具有十分重要意義,能夠幫助醫務工作者有效提高聽診與病況分析的準確性。
1 算法框架
IMCRA 和 OMLSA 算法在語音信號增強處理中已經獲得了較好的應用[17-19],本文將兩類算法相結合并首次應用于心音信號降噪處理中,提出的算法框架如圖1 所示。算法中,采用 IMCRA 通過兩次迭代追蹤短時窗內頻譜最小值。在較大窗平滑、較小窗跟蹤下能夠較好地追蹤噪聲變化情況,使得對噪聲譜估計的更新延遲更低,提高了追蹤的準確性;以最小值為噪聲估計且不進行補償,極大地減小過估計情況的出現。采用 OMLSA 算法,通過使目標估計心音與干凈心音間的差異最小化,較好地估計出干凈心音;在計算條件概率時,以 IMCRA 算法跟蹤計算的噪聲譜估計來計算后驗信噪比等參數,獲取最優頻譜增益函數。通過增益函數和帶噪心音的幅度譜即可估計出干凈心音幅度譜,經傅里葉逆變換和幀合成即可得到最終估計結果。

2 基于 IMCRA-OMLSA 的心音降噪
2.1 基于 OMLSA 的心音降噪
OMLSA 算法是一種最小化干凈信號和目標信號差異的譜估計方法[17, 19],不同于譜減法等經典譜估計,OMLSA 通過帶噪信號的頻譜直接估計目標信號的頻譜[20],能夠適應多種噪聲環境[21-22];同時,OMLSA 能夠避免音樂噪聲殘留,能夠有效保護較弱的聲學成分[23]。理論上,OMLSA 算法適用于心音降噪處理。
OMLSA 算法最早由 Cohen 等[24-25]提出。該算法基于心音與噪聲滿足獨立統計和高斯分布兩個假設提出[24, 26],在心音信號降噪中的流程如下:
(1)根據心音與噪聲獨立統計的假設,則帶噪心音如式(1)所示:
![]() |
其中,x(n) 為帶噪心音,s(n) 為干凈心音,d(n) 為噪聲。
對帶噪心音分幀加窗,并做傅里葉變換,其結果如式(2)所示:
![]() |
其中,k代表頻點,l代表幀數,和
分別為s(n) 和d(n) 分幀傅里葉變換的結果。
可表示為分幀的形式,如式(3)所示,相位譜
如式(4)所示:
![]() |
![]() |
上式中,win表示幀長,inc表示幀移,Lm 表示取虛部,Re 表示取實部,h(n)為窗函數。其中,窗函數一般取漢明窗,而不是矩形窗,這是因為漢明窗在卷積過程中旁瓣較小可減小振鈴效應,頻域更為平滑;心音具備短時平穩性,幀長win取 30 ms,為了減小計算,幀移inc取 15 ms,為幀長的一半。
(2)設干凈心音頻譜幅度為,目標估計心音頻譜幅度為
,按照最小化干凈心音與目標估計心音的差異[26],可得到頻譜幅度差異Delta,如式(5)、式(6)所示:
![]() |
![]() |
其中,min 表示最小值。
根據式(5)和文獻[26-27]的結論,可獲得估計的目標信號,如式(7)所示:
![]() |
則可表示出干凈心音頻譜幅度估計和帶噪信號頻譜幅度
、增益函數G和心音存在條件概率
的關系[21, 23],如式(8)、式(9)所示:
![]() |
![]() |
其中,為增益函數最小值,為一常數;
為心音存在條件下的增益函數。其中,增益函數、條件概率都需要通過一定方法來求取。
(3)估計增益函數[21, 27],如式(10)、式(11)所示:
![]() |
![]() |
其中,為先驗信噪比,
為后驗信噪比。
(4)心音存在的條件概率的估計如式 (12) 所示[26]:
![]() |
其中,為心音不存在的先驗概率。
(5)估計先驗信噪比[23]。
設為前一幀中心音存在條件下的增益函數,
為前一心音幀的后驗信噪比,則可計算當前心音幀的先驗信噪比
,如式(13)、式(14)所示:
![]() |
![]() |
其中,為噪聲功率譜。為了平衡噪聲殘留和心音失真的性能,這里引入了權重因子α,則可計算當前幀的先驗信噪比
,如式(15)所示:
![]() |
(6)估計心音不存在的先驗概率。
對的估計較為復雜,現有的研究中多使用軟判決(soft-decision)的方法進行估計[21]。這里直接給出
的估計方程式,如式(16)所示:
![]() |
其中,和
分別表示局部帶寬和整體帶寬上心音信號存在的似然概率[21-22],
為消除噪聲引入的參數。在求
和
之前,先定義先驗信噪比的遞歸平均
和其滑動平均
,如式(17)所示:
![]() |
其中,β 為平滑權重參數,? Wτ 和 Wτ 為卷積區間上限和下限,hτ(.)為歸一化漢明窗。則可計算心音信號存在的似然概率 ,如式(18)所示:
![]() |
其中,下標τ為“local”時表示局部帶寬似然概率,為“global”時表示整體帶寬似然概率[21],和
為最大、最小的判定閾值。
是通過先驗信噪比對心音幀內所有頻點平均得到的,如式(19)、式(20)所示:
![]() |
![]() |
其中,u(l) 的取值參考文獻[28]。
(7)噪聲估計。
噪聲估計在 OMLSA 降噪過程中十分重要,對噪聲估計的好壞往往決定了整個算法的降噪性能[20, 29]。對噪聲估計過高,則會造成心音失真;對噪聲估計過低,則會有較多的噪聲殘留。噪聲估計的方法主要包括基于最小值統計(minimum statistics,MS)的估計方法、最小值控制遞歸平均法(minimum control recursive average,MCRA)、IMCRA,其中屬 IMCRA 算法效果較好[21]。
2.2 基于 IMCRA 的噪聲估計
IMCRA 是在 MCRA 的基礎上進行的改進,將最小值控制遞歸過程由一次變為了兩次迭代[20-22]。噪聲估計是基于心音信號存在與否的假設,如式(21)所示[21]。
![]() |
心音信號存在與否的情況如圖2 所示,與時域波形圖相對應,心音存在情況即為 S1、S2 等心臟泵血過程中產生的信號,心音不存在情況
即為收縮期、舒張期成分中的平靜部分(不包含心雜音部分),圖2 還顯示了對數功率譜密度(power spectral density,PSD)的變化。

在 IMCRA 中,心音信號存在的條件概率可表示為式(22)、式(23)所示[29]:
![]() |
![]() |
其中,為心音存在的先驗概率,
和
分別為后驗信噪比和先驗信噪比。則可估計噪聲功率譜
,如式(24)、式(25)所示[20]:
![]() |
![]() |
因此,如何求取心音不存在時的先驗概率成為了關鍵所在。
的值與帶噪心音的平滑功率譜的最小值有關,在 IMCRA 中通過采用兩次迭代操作來跟蹤最小值[21]。第一次迭代的作用是對幀內頻點進行簡單的活動性檢測,第二次迭代則進行平滑操作去除心音存在時能量較強的部分(尤其是基礎心音存在的部分),從而使得可以通過較短的時窗來跟蹤功率譜的最小值。較短的時窗下可以使得對心音中噪聲的更新速度更快,這將更有利于算法去跟蹤非平穩噪聲[20, 26]。
第一次迭代過程中,根據判決準則粗略估計各頻點上心音的活動性,如式(26)所示:
![]() |
其中,和
為閾值參數,
為 1 時表示心音存在,為 0 時表示心音不存在。
及
的求取方式如式(27)、(28)所示。
![]() |
![]() |
其中,為最小噪聲估計偏差,
為長度為D的觀察窗內平滑后的帶噪心音功率譜
的最小值,求取方式如式 (29) 所示,平滑過程請參考文獻[21, 30]。
![]() |
在第二次迭代中,基于第一次迭代中的頻譜活動性檢測結果,除去功率譜較大的心音存在部分,對剩余的心音不存在的頻點進行平滑,平滑結果為,如式(30)所示[30]:
![]() |
其中,h(i) 為窗函數,WL為窗函數長度。在獲得后,還需要對其在時域上進行平滑,并在一個有限窗內跟蹤其最小值[26, 30],如式(31)~(32)所示:
![]() |
![]() |
其中,為平滑系數,b為最小值追蹤窗長,
為時域一階遞歸平滑結果,
為追蹤窗內最小值,則可求得心音信號不存在的先驗概率
[30],如式(33)~(35)所示:
![]() |
![]() |
![]() |
其中,為常數閾值。
在獲得心音不存在的先驗概率以后,就可以按式(22)所示求出心音存在的條件概率,按式(24)所示可求出對噪聲功率譜的估計
,按式(14)所示即可獲得后驗信噪比
,繼而獲得最優頻譜增益函數G。算法中一些參數的取值如表1 所示。其中,
和
為用于約束的條件概率下限和先驗概率上限。需要說明的是本文未對噪聲估計結果進行補償,即默認補償因子取為 1,因為這樣可以防止噪聲譜估計過大,造成心音失真。

3 降噪結果的主觀分析
實驗中測試了三類數據,第一類是直接通過相對干凈的心音直接加噪構建含噪心音,第二類是選用 2016 年心臟病學中物理網絡與計算挑戰賽的數據集(PhysioNet/Computing in Cardiology Challenge 2016,CinC 2016)(網址:
第一類心音選用了南京郵電大學成謝鋒教授課題組提供的心音,選用其中 30 例健康人群心音和 20 例心臟病患者人群心音(本文已獲權使用),數據的詳細描述可參考文獻[33-34]。以其中一條三尖瓣關閉不全的心音為例,實驗中將白噪聲以 10 dB 的信噪比加入原始心音中,降噪結果如圖3 所示。由圖3 中原始心音時域圖和頻譜圖不難發現,三尖瓣關閉不全時在 S1 和 S2 之間產生了較強的心雜音,主要由右心室收縮后血液反流至右心房產生[13-14],因此出現在收縮期的雜音不可被濾除。從 CEEMD 降噪結果的時域圖和頻譜圖不難發現,CEEMD 降噪中噪聲殘留現象較為嚴重,且對三尖瓣關閉不全的病理心雜音抑制性較大;從小波降噪結果的時域圖和頻譜圖可以看出,小波分析降噪同 CEEMD 降噪類似,能夠在一定程度上保留 S1 和 S2,但會抑制病理心雜音,且比 CEEMD 抑制的效果更明顯;而采用 IMCRA-OMLSA 算法的降噪結果的時域圖和頻譜圖上均顯示,噪聲得到了有效地抑制,且較好地保留了病理心雜音。

第二類心音選用 CinC 2016 中的心音,以其中 a0001 心音記錄為例,降噪結果如圖4 所示。從圖4 可以看出,CEEMD 降噪和小波分析降噪會殘留許多噪聲,頻譜可分析性較低;此外,從頻譜分析可看出,對于一些靠近心音的干擾噪聲,CEEMD 降噪和小波分析降噪并未做到有效抑制,如圖4 中橢圓標注處所示。從 IMCRA-OMLSA 算法的降噪結果來看,雖然在防止噪聲過估計的同時保留了許多細節波動,但整體上較好地增強了心音特征,且能夠抑制靠近心音的強干擾噪聲,頻譜的可分析性明顯提高。

第三類心音為自主采集心音,采集設備為蘿卜醫生電子智能聽診器(型號為 LBTZ-01),10 名知情同意受試者均為課題組成員,共采集 30 條心音記錄。以其中一人的心音為例,獲得降噪結果如圖5 所示。圖5 顯示,原始采集的心音中含有大量的低頻噪聲,主要來源于實驗室空調噪聲和設備自身電源的干擾噪聲。采用 CEEMD 和小波分析降噪的時域圖和頻域圖結果顯示,心音中同樣存在噪聲殘留。采用 IMCRA-OMLSA 算法的降噪結果顯示,頻帶內的噪聲得到了較好地抑制,頻譜圖中心音信號特征獲得了明顯突出,變得更為干凈。

因此,從算法降噪處理的結果進行主觀分析,可以得出以下結論:
(1)EMD 分解和小波分析等方法通過線性空間變換將心音信號分解為多個不同尺度下的分量,通過建立閾值函數進行濾波處理,抑制各尺度下的噪聲分量。在心音處理中的優點是能夠將頻帶劃分為多個尺度進行單獨處理,具備較高分辨率[3, 5];缺點是線性閾值判斷并不能夠準確區分心音與噪聲,將判斷為不是心音的部分大幅度抑制,對判斷為心音存在的部分完整保留,且不能很好地抑制一些非心音的高頻噪聲和靠近心音的強干擾噪聲。
(2)IMCRA-OMLSA 算法通過估計心音存在的條件概率、增益函數、噪聲譜等信息估計干凈心音,相比于空間變換和閾值濾波的方法,更為有效地利用了心音自身的信息,同時在時域和頻域角度分別進行了平滑估計運算,能夠從局部到整體快速追蹤噪聲,且不會出現噪聲嚴重的過估計現象,能夠有效地保護心雜音等弱小的異常特征,缺點是在防止噪聲譜過估計時不能夠很好地抑制非平穩噪聲。
4 降噪結果的客觀評估
心音信號處理中,研究人員均難以有效采集到完全干凈的心音數據,因此不可通過直接計算 SNR、RMSE 等客觀指標對算法進行嚴格地評價。故本文在主觀分析的基礎上,通過對正常與異常心音分類的貢獻性來客觀評估 IMCRA-OMLSA 算法降噪的有效性。正常與異常心音的分類過程主要涉及到訓練集和測試集的構建、特征提取以及分類器構建等內容[2, 31],降噪算法對分類系統貢獻性可直接通過分類結果評價指標直觀反映。
4.1 數據集與特征提取
實驗選用了 CinC 2016 數據集[31-32]。整個數據集由 a、b、c、d、e、f 這 6 個子數據庫構成,分別采集于 1 000 多位不同年齡、不同性別、不同身體狀況的受試者,共計 3 251 條心音記錄。
文獻[35]提出通過提取心動周期的多種特征(324 種)并進行融合來構建訓練集和測試集。本文采用此方法來提取心音的特征集,其步驟如下:
(1)采用 IMCRA-OMLSA 算法對心音降噪;
(2)采用基于邏輯回歸的隱半馬爾可夫模型(logistic regression hidden-semi-Markov model,LR-HSMM)[36]將各心音記錄的心動周期分割為 S1、收縮期、S2、舒張期;
(3)根據分割結果提取各心音記錄的 324 維特征;
(4)構建深度神經網絡,將特征集按照 9∶1 左右的比例隨機分為訓練集與測試集,訓練網絡并測試分類結果。
4.2 分類器構建
實驗中構建了具備兩個隱含層的深度神經網絡模型,如圖6 所示。為了確定最佳網絡結構,對隱含層的神經元數目進行了測試。可以確定的是輸入層神經元為 324,輸出層神經元為 2,隱含層數目為 2,隱含層神經元數量在 10~40 之間搜索,通過訓練結果確定。

4.3 評估結果及分析
將 2 950 條心音記錄的特征依次輸入分類模型中進行網絡訓練,迭代訓練至模型基本收斂后,在 301 條心音記錄的特征測試集下獲得測試結果。分類結果用到了靈敏度(sensitivity,Se)(以符號Se表示)、特異性(specificity,Sp)(以符號Sp表示)、總體得分(overall,Ov)(以符號Ov表示)三類評價指標,其定義如式(36)~式(38)所示:
![]() |
![]() |
![]() |
其中,樣本標識為異常(abnormal)狀態的心音記錄用A表示,樣本標識為正常(normal)狀態的心音記錄用N表示,樣本被分類為異常心音用a表示,樣本被分類為正常心音用n表示,樣本被分類為不確定(uncertain)類別的心音用u表示,可分析性較好的心音記錄用 1 標識,可分析性較差的心音用 2 標識,各類型的比例權重(wight)用w表示。wa1、wa2、wn1、wn2 為權重因子,wa1、wa2 的值為異常心音信號中可分析性較好和可分析性較差的心音記錄所占的百分比,wn1、wn2 的值為正常心音信號中可分析性較好和可分析性較差的心音記錄所占的百分比,取值如表2 所示[37]。Aa1、Aa2 分別為可分析性較好和可分析較差的異常心音被分類為異常心音的數量,An1、An2 分別為可分析性較好和可分析較差的異常心音被分類為正常心音的數量,Au1、Au2 分別為可分析性較好和可分析較差的異常心音被分類為不確定類別心音的數量;Na1、Na2 分別為可分析性較好和可分析較差的正常心音被分類為異常心音的數量,Nn1、Nn2 分別為可分析性較好和可分析較差的正常心音被分類為正常心音的數量,Nu1、Nu2 分別為可分析性較好和可分析較差的正常心音被分類為不確定類別心音的數量。

在帶通濾波[35]、小波分析[3]、CEEMD[5]以及 IMCRA-OMLSA 算法處理下,提取的特征集經網絡訓練后的測試結果如表3 所示。帶通濾波降噪處理下在隱含層神經元為(10, 10)時獲得最佳結果,Se = 0.802,Sp = 0.823,Ov = 0.812;小波分析降噪處理下在隱含層神經元為(40, 10)時獲得最佳結果,Se = 0.825,Sp = 0.837,Ov = 0.831;CEEMD 降噪處理下在隱含層神經元為(30, 10)時獲得最佳結果,Se = 0.824,Sp = 0.839,Ov = 0.831;IMCRA-OMLSA 降噪處理下在隱含層神經元為(30, 10)時獲得最佳結果,Se = 0.826,Sp = 0.843,Ov = 0.834。

此外,在最佳網絡狀態下,將四種降噪算法進行了對比,表3 結果表明,帶通濾波降噪處理后測試集在分類系統中得分最低,這是因為帶通濾波只是將低于 20 Hz、高于 120 Hz 的頻帶成分直接濾除了,丟掉了高于 120 Hz 的特征,且對 20~120 Hz 之間的部分并未進行去噪;CEEMD 降噪與小波分析降噪方法效果相近,但 CEEMD 效果較小波略好,原因在于 CEEMD 算法能夠較好地將信號分解為各模態函數分量,而小波分解的效果除了與信號本身有關,還與選取的小波基函數有關;IMCRA-OMLSA 降噪算法處理后,測試集在分類器中獲得了較好的效果,這是因為算法在保證心音特征不丟失的情況下很大限度地抑制了噪聲,能更為精確地提取心音的特征信息。
因此,通過降噪算法對分類的貢獻性進行客觀分析,即對分類結果客觀評價,可以得出結論:IMCRA-OMLSA 降噪算法相較于帶通濾波降噪、CEEMD 降噪、小波分析降噪等方法效果更好,能更有效地提高心音的可分析性,促進心音分類系統獲得更好的性能。
5 總結與討論
心音降噪可以提高心音信號的可分析性,在輔助診斷系統中起著重要的作用。針對經驗模式分解和小波分析降噪的局限性,本文提出了一種基于 IMCRA-OMLSA 的心音降噪方法,能夠在短時間內有效追蹤噪聲,動態估計長時采集的心音中的噪聲,并通過最佳頻譜增益函數和條件概率來逼近真實的干凈心音信號。實驗結果也進一步驗證,基于 IMCRA-OMLSA 的心音降噪方法能夠有效抑制采集環境下的噪聲,尤其是靠近基礎心音的噪聲,保護了反映生理特征的心雜音等成分,較 CEEMD 和小波分析的降噪算法效果更優,且能夠有效地提升分類系統的準確性。本文提出算法在心音降噪中具備良好的表現,為心音信號前端處理提供了一種新的方法,對降噪聽診器的設計、心血管疾病計算機輔助診斷系統的構建具有十分重要意義,能夠進一步促進心音信號在疾病診斷等方面的應用。
實驗算法中也存在兩個問題需進一步解決:① 在 IMCRA 算法的第一次迭代中,當活動性檢測效果不理想時,可能導致平滑處理結果偏大,因此可進一步增加約束條件來提升算法魯棒性;② 當心音信號不存在時,干凈心音的估計結果不為 0,而是最小增益與幅值的乘積,因此還可進一步考慮在 OMLSA 算法中增加心音不存在時的約束條件,使估計誤差較小且頻譜特征更為清晰。
利益沖突聲明:本文全體作者均聲明不存在利益沖突。
引言
心音信號分析是診斷心血管疾病的重要手段,被廣泛應用于臨床診斷和體檢當中[1-2]。心音信號的分析主要包括醫生經驗聽診和計算機輔助診斷[2-3],現今多借助于電子聽診器進行采集[4]。然而,在采集過程中,心音信號不可避免地會受到環境噪聲、采集設備噪聲以及人體其它器官噪聲等因素的干擾[3, 5],使得原本就十分弱小的心音信號的可分析性極大降低。因此,在對心音信號分析之前,需要先做降噪處理。
近幾十年來,研究和應用較廣的心音降噪方法當屬基于小波分析的方法[3, 6]和基于經驗模式分解(empirical mode decomposition,EMD)的方法[5]。其中,基于小波分析的方法是采用小波分解獲取小波域內不同尺度的小波系數,根據噪聲對應的小波系數小于心音對應的小波系數的假設,建立適當的閾值函數對各尺度的小波系數進行線性濾波處理,抑制噪聲的小波系數,達到降噪的目的[3, 6];此外,文獻[3, 5]也指出根據心音信號的頻帶特點,將部分高頻細節系數和低頻近似系數置零,能進一步去除不必要的噪聲。基于 EMD 的方法是將信號分解為多個模態固有函數(intrinsic mode function,IMF)分量,通過對模態混疊臨界點的判斷,將臨界點前的 IMF 分量視為噪聲并直接丟棄,同時丟棄心音成分少、噪聲含量大的高頻 IMF 分量[5, 7],并建立合適的濾波機制,將心音信號從保留的混疊 IMF 分量中分離出來,如互補總體 EMD(complementary ensemble EMD,CEEMD)[5]。然而,實際采集環境下的干擾噪聲具備隨機、非平穩、混沌等特性[8],且心音具備多樣性[9],小波系數的大小難以精準區分心音信號和噪聲(如類似心音的噪聲)[10],這使得獲取精準的自適應閾值函數更為困難,因此要想通過小波或 EMD 的空間變換和線性濾波處理將噪聲和心音完全分離開來十分困難,且容易出現心音失真情況[10-11]。同時,近些年的研究表明[4, 12-13],高頻部分也蘊含著有用的特征,尤其是在一些異常心音(如病變心雜音)中,故不可簡單濾波或者直接丟棄。此外,心音降噪還涉及到心雜音等由心臟異常引起的額外音,且許多額外音與噪聲特征相近,如三尖瓣關閉不全時在第一心音(S1)與第二心音(S2)間產生的心雜音[14-15]。作為病理特征的心雜音對心音分類有著重要的作用,不可直接濾除[13],這也進一步加大了心音降噪的難度。在心音降噪處理中,應以心音不失真為前提,最大程度地抑制噪聲[3, 5]。
實際中無法直接獲得完全干凈的心音信號,這就使得各類算法的降噪結果沒有標準的參考。目前常用的評價方法主要包括以相對干凈的心音信號為參考計算評價指標,如信噪比(signal-to-noise ratio,SNR)、均方根誤差(root mean square error,RMSE)[3, 5]、主觀譜圖分析[5]、平均意見得分(mean opinion score,MOS)[16]。其中,以相對干凈的心音信號為參考計算評價指標的分析并不嚴格,其客觀性并不強。
針對心音降噪面臨的問題,本文提出了一種融合改進最小值控制遞歸平均(improved minimum control recursive average,IMCRA)和最優修正對數譜幅度(optimally modified log-spectral amplitude,OMLSA)估計的心音降噪方法(后文簡稱:IMCRA-OMLSA)。該方法能夠從不同于小波分析和 EMD 分解降噪方法的角度,通過最優增益函數譜估計,最小化目標估計心音與干凈心音的差異,有效地抑制噪聲,降低心音失真的風險。在對算法的評價中,由于完全干凈的心音不可知,故提出譜圖主觀分析和降噪算法對分類系統貢獻度的客觀分析相結合的評價機制,確保對算法評價的客觀性及有效性。本文提出算法在心音噪聲中具備良好的性能,該算法對降噪聽診器的設計、心血管疾病計算機輔助診斷系統的構建具有十分重要意義,能夠幫助醫務工作者有效提高聽診與病況分析的準確性。
1 算法框架
IMCRA 和 OMLSA 算法在語音信號增強處理中已經獲得了較好的應用[17-19],本文將兩類算法相結合并首次應用于心音信號降噪處理中,提出的算法框架如圖1 所示。算法中,采用 IMCRA 通過兩次迭代追蹤短時窗內頻譜最小值。在較大窗平滑、較小窗跟蹤下能夠較好地追蹤噪聲變化情況,使得對噪聲譜估計的更新延遲更低,提高了追蹤的準確性;以最小值為噪聲估計且不進行補償,極大地減小過估計情況的出現。采用 OMLSA 算法,通過使目標估計心音與干凈心音間的差異最小化,較好地估計出干凈心音;在計算條件概率時,以 IMCRA 算法跟蹤計算的噪聲譜估計來計算后驗信噪比等參數,獲取最優頻譜增益函數。通過增益函數和帶噪心音的幅度譜即可估計出干凈心音幅度譜,經傅里葉逆變換和幀合成即可得到最終估計結果。

2 基于 IMCRA-OMLSA 的心音降噪
2.1 基于 OMLSA 的心音降噪
OMLSA 算法是一種最小化干凈信號和目標信號差異的譜估計方法[17, 19],不同于譜減法等經典譜估計,OMLSA 通過帶噪信號的頻譜直接估計目標信號的頻譜[20],能夠適應多種噪聲環境[21-22];同時,OMLSA 能夠避免音樂噪聲殘留,能夠有效保護較弱的聲學成分[23]。理論上,OMLSA 算法適用于心音降噪處理。
OMLSA 算法最早由 Cohen 等[24-25]提出。該算法基于心音與噪聲滿足獨立統計和高斯分布兩個假設提出[24, 26],在心音信號降噪中的流程如下:
(1)根據心音與噪聲獨立統計的假設,則帶噪心音如式(1)所示:
![]() |
其中,x(n) 為帶噪心音,s(n) 為干凈心音,d(n) 為噪聲。
對帶噪心音分幀加窗,并做傅里葉變換,其結果如式(2)所示:
![]() |
其中,k代表頻點,l代表幀數,和
分別為s(n) 和d(n) 分幀傅里葉變換的結果。
可表示為分幀的形式,如式(3)所示,相位譜
如式(4)所示:
![]() |
![]() |
上式中,win表示幀長,inc表示幀移,Lm 表示取虛部,Re 表示取實部,h(n)為窗函數。其中,窗函數一般取漢明窗,而不是矩形窗,這是因為漢明窗在卷積過程中旁瓣較小可減小振鈴效應,頻域更為平滑;心音具備短時平穩性,幀長win取 30 ms,為了減小計算,幀移inc取 15 ms,為幀長的一半。
(2)設干凈心音頻譜幅度為,目標估計心音頻譜幅度為
,按照最小化干凈心音與目標估計心音的差異[26],可得到頻譜幅度差異Delta,如式(5)、式(6)所示:
![]() |
![]() |
其中,min 表示最小值。
根據式(5)和文獻[26-27]的結論,可獲得估計的目標信號,如式(7)所示:
![]() |
則可表示出干凈心音頻譜幅度估計和帶噪信號頻譜幅度
、增益函數G和心音存在條件概率
的關系[21, 23],如式(8)、式(9)所示:
![]() |
![]() |
其中,為增益函數最小值,為一常數;
為心音存在條件下的增益函數。其中,增益函數、條件概率都需要通過一定方法來求取。
(3)估計增益函數[21, 27],如式(10)、式(11)所示:
![]() |
![]() |
其中,為先驗信噪比,
為后驗信噪比。
(4)心音存在的條件概率的估計如式 (12) 所示[26]:
![]() |
其中,為心音不存在的先驗概率。
(5)估計先驗信噪比[23]。
設為前一幀中心音存在條件下的增益函數,
為前一心音幀的后驗信噪比,則可計算當前心音幀的先驗信噪比
,如式(13)、式(14)所示:
![]() |
![]() |
其中,為噪聲功率譜。為了平衡噪聲殘留和心音失真的性能,這里引入了權重因子α,則可計算當前幀的先驗信噪比
,如式(15)所示:
![]() |
(6)估計心音不存在的先驗概率。
對的估計較為復雜,現有的研究中多使用軟判決(soft-decision)的方法進行估計[21]。這里直接給出
的估計方程式,如式(16)所示:
![]() |
其中,和
分別表示局部帶寬和整體帶寬上心音信號存在的似然概率[21-22],
為消除噪聲引入的參數。在求
和
之前,先定義先驗信噪比的遞歸平均
和其滑動平均
,如式(17)所示:
![]() |
其中,β 為平滑權重參數,? Wτ 和 Wτ 為卷積區間上限和下限,hτ(.)為歸一化漢明窗。則可計算心音信號存在的似然概率 ,如式(18)所示:
![]() |
其中,下標τ為“local”時表示局部帶寬似然概率,為“global”時表示整體帶寬似然概率[21],和
為最大、最小的判定閾值。
是通過先驗信噪比對心音幀內所有頻點平均得到的,如式(19)、式(20)所示:
![]() |
![]() |
其中,u(l) 的取值參考文獻[28]。
(7)噪聲估計。
噪聲估計在 OMLSA 降噪過程中十分重要,對噪聲估計的好壞往往決定了整個算法的降噪性能[20, 29]。對噪聲估計過高,則會造成心音失真;對噪聲估計過低,則會有較多的噪聲殘留。噪聲估計的方法主要包括基于最小值統計(minimum statistics,MS)的估計方法、最小值控制遞歸平均法(minimum control recursive average,MCRA)、IMCRA,其中屬 IMCRA 算法效果較好[21]。
2.2 基于 IMCRA 的噪聲估計
IMCRA 是在 MCRA 的基礎上進行的改進,將最小值控制遞歸過程由一次變為了兩次迭代[20-22]。噪聲估計是基于心音信號存在與否的假設,如式(21)所示[21]。
![]() |
心音信號存在與否的情況如圖2 所示,與時域波形圖相對應,心音存在情況即為 S1、S2 等心臟泵血過程中產生的信號,心音不存在情況
即為收縮期、舒張期成分中的平靜部分(不包含心雜音部分),圖2 還顯示了對數功率譜密度(power spectral density,PSD)的變化。

在 IMCRA 中,心音信號存在的條件概率可表示為式(22)、式(23)所示[29]:
![]() |
![]() |
其中,為心音存在的先驗概率,
和
分別為后驗信噪比和先驗信噪比。則可估計噪聲功率譜
,如式(24)、式(25)所示[20]:
![]() |
![]() |
因此,如何求取心音不存在時的先驗概率成為了關鍵所在。
的值與帶噪心音的平滑功率譜的最小值有關,在 IMCRA 中通過采用兩次迭代操作來跟蹤最小值[21]。第一次迭代的作用是對幀內頻點進行簡單的活動性檢測,第二次迭代則進行平滑操作去除心音存在時能量較強的部分(尤其是基礎心音存在的部分),從而使得可以通過較短的時窗來跟蹤功率譜的最小值。較短的時窗下可以使得對心音中噪聲的更新速度更快,這將更有利于算法去跟蹤非平穩噪聲[20, 26]。
第一次迭代過程中,根據判決準則粗略估計各頻點上心音的活動性,如式(26)所示:
![]() |
其中,和
為閾值參數,
為 1 時表示心音存在,為 0 時表示心音不存在。
及
的求取方式如式(27)、(28)所示。
![]() |
![]() |
其中,為最小噪聲估計偏差,
為長度為D的觀察窗內平滑后的帶噪心音功率譜
的最小值,求取方式如式 (29) 所示,平滑過程請參考文獻[21, 30]。
![]() |
在第二次迭代中,基于第一次迭代中的頻譜活動性檢測結果,除去功率譜較大的心音存在部分,對剩余的心音不存在的頻點進行平滑,平滑結果為,如式(30)所示[30]:
![]() |
其中,h(i) 為窗函數,WL為窗函數長度。在獲得后,還需要對其在時域上進行平滑,并在一個有限窗內跟蹤其最小值[26, 30],如式(31)~(32)所示:
![]() |
![]() |
其中,為平滑系數,b為最小值追蹤窗長,
為時域一階遞歸平滑結果,
為追蹤窗內最小值,則可求得心音信號不存在的先驗概率
[30],如式(33)~(35)所示:
![]() |
![]() |
![]() |
其中,為常數閾值。
在獲得心音不存在的先驗概率以后,就可以按式(22)所示求出心音存在的條件概率,按式(24)所示可求出對噪聲功率譜的估計
,按式(14)所示即可獲得后驗信噪比
,繼而獲得最優頻譜增益函數G。算法中一些參數的取值如表1 所示。其中,
和
為用于約束的條件概率下限和先驗概率上限。需要說明的是本文未對噪聲估計結果進行補償,即默認補償因子取為 1,因為這樣可以防止噪聲譜估計過大,造成心音失真。

3 降噪結果的主觀分析
實驗中測試了三類數據,第一類是直接通過相對干凈的心音直接加噪構建含噪心音,第二類是選用 2016 年心臟病學中物理網絡與計算挑戰賽的數據集(PhysioNet/Computing in Cardiology Challenge 2016,CinC 2016)(網址:
第一類心音選用了南京郵電大學成謝鋒教授課題組提供的心音,選用其中 30 例健康人群心音和 20 例心臟病患者人群心音(本文已獲權使用),數據的詳細描述可參考文獻[33-34]。以其中一條三尖瓣關閉不全的心音為例,實驗中將白噪聲以 10 dB 的信噪比加入原始心音中,降噪結果如圖3 所示。由圖3 中原始心音時域圖和頻譜圖不難發現,三尖瓣關閉不全時在 S1 和 S2 之間產生了較強的心雜音,主要由右心室收縮后血液反流至右心房產生[13-14],因此出現在收縮期的雜音不可被濾除。從 CEEMD 降噪結果的時域圖和頻譜圖不難發現,CEEMD 降噪中噪聲殘留現象較為嚴重,且對三尖瓣關閉不全的病理心雜音抑制性較大;從小波降噪結果的時域圖和頻譜圖可以看出,小波分析降噪同 CEEMD 降噪類似,能夠在一定程度上保留 S1 和 S2,但會抑制病理心雜音,且比 CEEMD 抑制的效果更明顯;而采用 IMCRA-OMLSA 算法的降噪結果的時域圖和頻譜圖上均顯示,噪聲得到了有效地抑制,且較好地保留了病理心雜音。

第二類心音選用 CinC 2016 中的心音,以其中 a0001 心音記錄為例,降噪結果如圖4 所示。從圖4 可以看出,CEEMD 降噪和小波分析降噪會殘留許多噪聲,頻譜可分析性較低;此外,從頻譜分析可看出,對于一些靠近心音的干擾噪聲,CEEMD 降噪和小波分析降噪并未做到有效抑制,如圖4 中橢圓標注處所示。從 IMCRA-OMLSA 算法的降噪結果來看,雖然在防止噪聲過估計的同時保留了許多細節波動,但整體上較好地增強了心音特征,且能夠抑制靠近心音的強干擾噪聲,頻譜的可分析性明顯提高。

第三類心音為自主采集心音,采集設備為蘿卜醫生電子智能聽診器(型號為 LBTZ-01),10 名知情同意受試者均為課題組成員,共采集 30 條心音記錄。以其中一人的心音為例,獲得降噪結果如圖5 所示。圖5 顯示,原始采集的心音中含有大量的低頻噪聲,主要來源于實驗室空調噪聲和設備自身電源的干擾噪聲。采用 CEEMD 和小波分析降噪的時域圖和頻域圖結果顯示,心音中同樣存在噪聲殘留。采用 IMCRA-OMLSA 算法的降噪結果顯示,頻帶內的噪聲得到了較好地抑制,頻譜圖中心音信號特征獲得了明顯突出,變得更為干凈。

因此,從算法降噪處理的結果進行主觀分析,可以得出以下結論:
(1)EMD 分解和小波分析等方法通過線性空間變換將心音信號分解為多個不同尺度下的分量,通過建立閾值函數進行濾波處理,抑制各尺度下的噪聲分量。在心音處理中的優點是能夠將頻帶劃分為多個尺度進行單獨處理,具備較高分辨率[3, 5];缺點是線性閾值判斷并不能夠準確區分心音與噪聲,將判斷為不是心音的部分大幅度抑制,對判斷為心音存在的部分完整保留,且不能很好地抑制一些非心音的高頻噪聲和靠近心音的強干擾噪聲。
(2)IMCRA-OMLSA 算法通過估計心音存在的條件概率、增益函數、噪聲譜等信息估計干凈心音,相比于空間變換和閾值濾波的方法,更為有效地利用了心音自身的信息,同時在時域和頻域角度分別進行了平滑估計運算,能夠從局部到整體快速追蹤噪聲,且不會出現噪聲嚴重的過估計現象,能夠有效地保護心雜音等弱小的異常特征,缺點是在防止噪聲譜過估計時不能夠很好地抑制非平穩噪聲。
4 降噪結果的客觀評估
心音信號處理中,研究人員均難以有效采集到完全干凈的心音數據,因此不可通過直接計算 SNR、RMSE 等客觀指標對算法進行嚴格地評價。故本文在主觀分析的基礎上,通過對正常與異常心音分類的貢獻性來客觀評估 IMCRA-OMLSA 算法降噪的有效性。正常與異常心音的分類過程主要涉及到訓練集和測試集的構建、特征提取以及分類器構建等內容[2, 31],降噪算法對分類系統貢獻性可直接通過分類結果評價指標直觀反映。
4.1 數據集與特征提取
實驗選用了 CinC 2016 數據集[31-32]。整個數據集由 a、b、c、d、e、f 這 6 個子數據庫構成,分別采集于 1 000 多位不同年齡、不同性別、不同身體狀況的受試者,共計 3 251 條心音記錄。
文獻[35]提出通過提取心動周期的多種特征(324 種)并進行融合來構建訓練集和測試集。本文采用此方法來提取心音的特征集,其步驟如下:
(1)采用 IMCRA-OMLSA 算法對心音降噪;
(2)采用基于邏輯回歸的隱半馬爾可夫模型(logistic regression hidden-semi-Markov model,LR-HSMM)[36]將各心音記錄的心動周期分割為 S1、收縮期、S2、舒張期;
(3)根據分割結果提取各心音記錄的 324 維特征;
(4)構建深度神經網絡,將特征集按照 9∶1 左右的比例隨機分為訓練集與測試集,訓練網絡并測試分類結果。
4.2 分類器構建
實驗中構建了具備兩個隱含層的深度神經網絡模型,如圖6 所示。為了確定最佳網絡結構,對隱含層的神經元數目進行了測試。可以確定的是輸入層神經元為 324,輸出層神經元為 2,隱含層數目為 2,隱含層神經元數量在 10~40 之間搜索,通過訓練結果確定。

4.3 評估結果及分析
將 2 950 條心音記錄的特征依次輸入分類模型中進行網絡訓練,迭代訓練至模型基本收斂后,在 301 條心音記錄的特征測試集下獲得測試結果。分類結果用到了靈敏度(sensitivity,Se)(以符號Se表示)、特異性(specificity,Sp)(以符號Sp表示)、總體得分(overall,Ov)(以符號Ov表示)三類評價指標,其定義如式(36)~式(38)所示:
![]() |
![]() |
![]() |
其中,樣本標識為異常(abnormal)狀態的心音記錄用A表示,樣本標識為正常(normal)狀態的心音記錄用N表示,樣本被分類為異常心音用a表示,樣本被分類為正常心音用n表示,樣本被分類為不確定(uncertain)類別的心音用u表示,可分析性較好的心音記錄用 1 標識,可分析性較差的心音用 2 標識,各類型的比例權重(wight)用w表示。wa1、wa2、wn1、wn2 為權重因子,wa1、wa2 的值為異常心音信號中可分析性較好和可分析性較差的心音記錄所占的百分比,wn1、wn2 的值為正常心音信號中可分析性較好和可分析性較差的心音記錄所占的百分比,取值如表2 所示[37]。Aa1、Aa2 分別為可分析性較好和可分析較差的異常心音被分類為異常心音的數量,An1、An2 分別為可分析性較好和可分析較差的異常心音被分類為正常心音的數量,Au1、Au2 分別為可分析性較好和可分析較差的異常心音被分類為不確定類別心音的數量;Na1、Na2 分別為可分析性較好和可分析較差的正常心音被分類為異常心音的數量,Nn1、Nn2 分別為可分析性較好和可分析較差的正常心音被分類為正常心音的數量,Nu1、Nu2 分別為可分析性較好和可分析較差的正常心音被分類為不確定類別心音的數量。

在帶通濾波[35]、小波分析[3]、CEEMD[5]以及 IMCRA-OMLSA 算法處理下,提取的特征集經網絡訓練后的測試結果如表3 所示。帶通濾波降噪處理下在隱含層神經元為(10, 10)時獲得最佳結果,Se = 0.802,Sp = 0.823,Ov = 0.812;小波分析降噪處理下在隱含層神經元為(40, 10)時獲得最佳結果,Se = 0.825,Sp = 0.837,Ov = 0.831;CEEMD 降噪處理下在隱含層神經元為(30, 10)時獲得最佳結果,Se = 0.824,Sp = 0.839,Ov = 0.831;IMCRA-OMLSA 降噪處理下在隱含層神經元為(30, 10)時獲得最佳結果,Se = 0.826,Sp = 0.843,Ov = 0.834。

此外,在最佳網絡狀態下,將四種降噪算法進行了對比,表3 結果表明,帶通濾波降噪處理后測試集在分類系統中得分最低,這是因為帶通濾波只是將低于 20 Hz、高于 120 Hz 的頻帶成分直接濾除了,丟掉了高于 120 Hz 的特征,且對 20~120 Hz 之間的部分并未進行去噪;CEEMD 降噪與小波分析降噪方法效果相近,但 CEEMD 效果較小波略好,原因在于 CEEMD 算法能夠較好地將信號分解為各模態函數分量,而小波分解的效果除了與信號本身有關,還與選取的小波基函數有關;IMCRA-OMLSA 降噪算法處理后,測試集在分類器中獲得了較好的效果,這是因為算法在保證心音特征不丟失的情況下很大限度地抑制了噪聲,能更為精確地提取心音的特征信息。
因此,通過降噪算法對分類的貢獻性進行客觀分析,即對分類結果客觀評價,可以得出結論:IMCRA-OMLSA 降噪算法相較于帶通濾波降噪、CEEMD 降噪、小波分析降噪等方法效果更好,能更有效地提高心音的可分析性,促進心音分類系統獲得更好的性能。
5 總結與討論
心音降噪可以提高心音信號的可分析性,在輔助診斷系統中起著重要的作用。針對經驗模式分解和小波分析降噪的局限性,本文提出了一種基于 IMCRA-OMLSA 的心音降噪方法,能夠在短時間內有效追蹤噪聲,動態估計長時采集的心音中的噪聲,并通過最佳頻譜增益函數和條件概率來逼近真實的干凈心音信號。實驗結果也進一步驗證,基于 IMCRA-OMLSA 的心音降噪方法能夠有效抑制采集環境下的噪聲,尤其是靠近基礎心音的噪聲,保護了反映生理特征的心雜音等成分,較 CEEMD 和小波分析的降噪算法效果更優,且能夠有效地提升分類系統的準確性。本文提出算法在心音降噪中具備良好的表現,為心音信號前端處理提供了一種新的方法,對降噪聽診器的設計、心血管疾病計算機輔助診斷系統的構建具有十分重要意義,能夠進一步促進心音信號在疾病診斷等方面的應用。
實驗算法中也存在兩個問題需進一步解決:① 在 IMCRA 算法的第一次迭代中,當活動性檢測效果不理想時,可能導致平滑處理結果偏大,因此可進一步增加約束條件來提升算法魯棒性;② 當心音信號不存在時,干凈心音的估計結果不為 0,而是最小增益與幅值的乘積,因此還可進一步考慮在 OMLSA 算法中增加心音不存在時的約束條件,使估計誤差較小且頻譜特征更為清晰。
利益沖突聲明:本文全體作者均聲明不存在利益沖突。