基于活動性檢測動態(tài)估計噪聲的心音降噪算法

許春冬，辛鵬麗，閔 源+，應冬文，2，周 靜，李海兵

(1.江西理工大學 信息工程學院，江西 贛州 341000；2.中國科學院大學 電子電氣與通信工程學院，北京 100049)

0 引 言

心音信號蘊含著大量心臟和周邊血管功能狀態(tài)的生理與病理信息[1，2]，是心血管疾病臨床診斷的重要依據(jù)之一[3]。然而，在進行心音信號的采集與處理時，由于心音信號能量較弱[4]，致使其極易受到來自外部環(huán)境噪聲、采集設備自噪聲及人體其它生理器官雜音等噪聲的干擾[1，2，5]。這些噪聲的存在往往會導致心音信號的質(zhì)量降低，進而對后續(xù)的分析與診斷結果造成嚴重的干擾[1，6]。

近十幾年來，國內(nèi)外學者提出了許多心音降噪的方法，主要可分為3類：基于小波的心音降噪方法、基于經(jīng)驗模態(tài)分解的心音降噪方法、基于噪聲功率估計的心音降噪方法[7-10]。雖然上述算法均取得了一定的降噪效果，但仍存在一定的不足，例如基于小波的心音降噪算法[7，8]在處理異常心音時易造成心音丟失；基于經(jīng)驗模態(tài)分解的降噪方法[9，10]容易使降噪后的心音信號出現(xiàn)失真或低、中頻噪聲殘留?；谠肼晞討B(tài)估計的方法[11]存在低頻噪聲殘留和部分心音失真的問題。

針對上述心音降噪算法所存在的不足，本文提出了一種基于心音活動性檢測動態(tài)估計噪聲的心音降噪算法(后文中均稱為HSAD-IMCRA算法)。該方法通過設計的HSAD方法來判斷心音幀是否為基礎心音，并對基礎心音和非基礎心音中的噪聲采用不同的方法進行估計，從而克服了傳統(tǒng)噪聲估計的不足，能夠在保證心音信號不失真的情況下，實現(xiàn)最大程度的噪聲抑制。

1 算法框架

提出算法的流程如圖1所示。首先，該算法將時域含噪心音信號經(jīng)短時傅里葉變換(short-time Fourier transform，STFT)變換到頻域，并根據(jù)心音信號的結構特點[12]，設計適合的HSAD方法，判別基礎心音(fundamental heart sound，F(xiàn)HS)與非基礎心音(non-FHS)成分；其次，分別采用IMCRA和遞歸平滑的方法進行噪聲估計；然后，通過估計的噪聲功率譜與含噪心音的前后幀信息來估計非因果先驗信噪比(signal-to-noise ratio，SNR)，求取干凈心音估計器的增益函數(shù)GLSA[13]；最后，結合原始相位信息，通過逆STFT(ISTFT)重構時域降噪心音。其中，采用非因果SNR估計方式代替因果SNR的估計方式，更準確地估計了FHS出現(xiàn)和截止邊緣部分及一些異常突變心音信號的先驗SNR[14，15]，在一定程度上減少了FHS的失真；采用HSAD輔助噪聲估計，能更準確地跟蹤噪聲功率的變化情況，從而減少了噪聲的殘留。

圖1 心音降噪算法流程

2 算法原理

2.1 MMSE-LSA算法原理

由Ephraim等[16]提出的最小均方誤差對數(shù)幅度譜(minimum mean square error of log-spectral amplitude，MMSE-LSA)估計算法通過對聲學信號進行最優(yōu)幅度譜估計[17]來恢復干凈信號，相較于一些經(jīng)典的降噪算法，能更有效降低信號失真[18]，提高信號的可分析性，故適用于心音信號的降噪處理。假設心音與噪聲不相關[13]，設s(n) 為干凈心音，d(n) 為噪聲，則含噪心音y(n) 可表示為

y(n)=s(n)+d(n)

(1)

對y(n) 進行短時傅里葉變換[12，13]，可得

Y(k，l)=S(k，l)+D(k，l)

(2)

將式(2)表示為極坐標形式[13]

R(k，l)ejφ(k，l)=A(k，l)ejφ(k，l)+N(k，l)ej?(k，l)

(3)

其中，l和k分別代表幀數(shù)和頻點數(shù)的索引，Y(k，l)、S(k，l)和D(k，l) 分別是y(n)、s(n) 和d(n) 進行STFT后的結果，R(k，l) 和φ(k，l) 分別是Y(k，l) 的頻譜幅度和相位；A(k，l) 和φ(k，l) 分別是S(k，l) 的頻譜幅度和相位；N(k，l) 和?(k，l) 分別是D(k，l) 的頻譜幅度和相位。

(4)

(5)

(6)

其中，ζ1(k，l)=λS(k，l)/λD(k，l) 為先驗信噪比，γ1(k，l)=|Y(k，l)|2/λD(k，l) 為后驗信噪比，λS(k，l) 為心音功率譜，λD(k，l) 為噪聲功率譜。

2.2 非因果先驗信噪比估計

在增益函數(shù)GLSA(k，l) 的求解過程中，用非因果信噪比[15]代替因果信噪比，能更準確估計存在突變心音信號或脈沖干擾噪聲時[13]的噪聲功率。在非因果信噪比的估計中，利用第0至l+M幀含噪心音來估計第l幀的先驗信噪比[15]，而不是局限于使用前一幀信息。

(7)

(8)

(9)

2.3 噪聲估計

2.3.1 IMCRA噪聲估計算法

(10)

基于最小均方誤差(minimum mean square error，MMSE)的準則下[13]，噪聲功率譜的估計如式(11)所示[13，20]

(11)

(12)

α2=α+(1-α)·P(k，l)

(13)

式(13)中，心音存在條件概率P(k，l)， 可通過貝葉斯公式[22]計算得到，如式(14)和式(15)所示

(14)

(15)

第一次迭代過程中，根據(jù)判決準則對是否存在基礎心音進行粗略估計[19]，如式(16)所示

(16)

其中，γ0和ζ0為常數(shù)閾值參數(shù)，基礎心音存在時I(k，l) 為1，基礎心音不存在時I(k，l) 為0。γmin(k，l) 和ζ2(k，l) 的定義如式(17)和式(18)所示

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

其中，γ3為常數(shù)閾值。

2.3.2 基于HSAD-IMCRA的噪聲功率估計

傳統(tǒng)的基于MMSE-LSA的降噪算法采用IMCRA方法對信號噪聲功率譜進行估計，但IMCRA噪聲估計對基礎心音開始邊緣部位的先/后驗信噪比參數(shù)估計往往偏低，易造成降噪后的心音信號出現(xiàn)基礎心音成分失真。此外，若心動周期內(nèi)存在人體其它器官引起的額外噪聲，也易造成先/后驗信噪比偏高，使得噪聲殘留。而心音信號主要是由第一心音(S1)、收縮期(Sys)、第二心音(S2)、舒張期(Dia)這4個狀態(tài)組成[9]，且具備一定的周期性[7，9]，可對應圖2所示。其中，Sys和Dia一般不包含基礎心音信號成分，表現(xiàn)為噪聲特性，故可根據(jù)遞歸平滑法來估計并更新噪聲功率譜[20]，獲得更為準確的先/后驗信噪比估計結果。而S1和S2持續(xù)時間段是由噪聲和基礎心音信號混疊構成，故采用IMCRA對噪聲進行有效追蹤，動態(tài)估計噪聲功率譜并更新[22]。需要額外指出的是，IMCRA噪聲估計過程中，采用2.2節(jié)所述非因果信噪比計算方式替換因果信噪比的計算方式，減小心音失真情況的出現(xiàn)。

圖2 心音信號時域波形

故本文算法中采用HSAD-IMCRA的算法來估計當前含噪心音幀的噪聲功率譜，能進一步降低估計的干凈心音的失真和噪聲殘留。所提算法噪聲功率譜估計流程如圖3所示。

圖3 HSAD-IMCRA算法噪聲功率譜估計流程

活動性檢測是聲學信號處理領域里較為常用的處理手段，且類別繁多[13，20]。為了在有效檢測心音活動性的前提下進一步考慮低的算法復雜度，本文采用對數(shù)頻譜距離(logarithmic amplitude spectrum distance，LASD)來設計適合心音信號的HSAD方法，用于判別基礎心音幀和非基礎心音幀[21]?；谇皫讕瑸榉腔A心音幀的假設，可通過式(28)將LASD計算如下

(28)

(29)

T(l)=τD(l-1)+(τ-1)D(l)，τ>1&T(l)>D(1)

(30)

(31)

其中，ρ(l-1)=1表示前一幀為基礎心音幀，ρ(l-1)=0表示前一幀為非基礎心音幀，NIMCRA(k，l) 為IMCRA算法估計的噪聲幅度譜，τ為自適應閾值平滑參數(shù)，αn為平滑因子[20]。

此外，根據(jù)式(31)的原理，噪聲功率譜的平滑估計λD(k，l) 可表示如下

(32)

其中，λIMCRA(k，l) 為IMCRA估計的噪聲功率譜，在獲得估計的噪聲功率譜后，即可計算先/后驗信噪比，并按式(5)獲得增益函數(shù)GLSA，從而重構出干凈心音信號。實驗中所涉及的參數(shù)設置見表1。

表1 參數(shù)的取值

3 實驗仿真與結果分析

3.1 實驗數(shù)據(jù)

實驗仿真所用數(shù)據(jù)集由公開數(shù)據(jù)集(https：//github.com/yaseen21khan/Classification-of-Heart-Sound-Si-gnal-Using-Multiple-Features-)[23]和自采集心音數(shù)據(jù)集[24]兩部分組成，共計1030條。參考文獻[10]指出，高斯白噪聲是心音信號中的主要噪聲成分，故本文主要以去除高斯白噪聲為主。在實驗仿真中，從構建的數(shù)據(jù)集中抽取出80條相對干凈的正常心音、320條相對干凈的異常心音，并將高斯白噪聲分別以0 dB、5 dB、10 dB、15 dB、20 dB的SNR疊加到抽取的心音信號中[10]，構建含噪心音數(shù)據(jù)集，且所有心音信號均降采樣至2000 Hz。

3.2 HSAD-IMCRA噪聲估計性能分析

為了驗證HSAD-IMCRA算法的噪聲估計性能，將其與IMCRA噪聲估計性能進行了對比。選用的評測指標為對數(shù)似然比(log likelihood ratio，LLR)和加權譜斜率(weighted spectral slope，WSS)[13，20]。其中，LLR是對頻譜距離的度量，取值范圍在0～2之間，估計的干凈心音與原始干凈心音間的LLR越小，即噪聲估計越準確[13]；WSS是對頻帶譜斜率之間加權差距的度量，WSS越小，則估計的干凈心音與原始干凈心音越接近，即噪聲估計越準確[20]。噪聲估計算法相關指標對比結果如圖4所示。

圖4 噪聲估計算法相關指標結果對比

圖4(a)和圖4(b)分別是基于IMCRA算法和基于HSAD-IMCRA算法估計的干凈心音與真實的干凈心音間的LLR與WSS的對比結果。顯然，基于HSAD-IMCRA算法的WSS和LLR評測結果均低于IMCRA算法，這表明基于HSAD-IMCRA算法的噪聲估計性能要優(yōu)于IMCRA算法。

3.3 各類算法降噪結果的時頻圖對比

提出算法與基于小波閾值的降噪算法[8]、基于OMLSA-小波的降噪算法[11]進行實驗對比。其中，基于小波閾值降噪算法中采用的小波基為coif-5，分解層數(shù)為5層，閾值函數(shù)為自適應非線性閾值；基于OMLSA-小波降噪算法中，先通過OMLSA對含噪心音降噪，再利用小波降噪算法對OMLSA降噪后的心音繼續(xù)降噪處理，其中小波降噪算法中的相關參數(shù)參考文獻[11]。算法降噪性能的評測指標為輸出SNR和均方根誤差(root mean square error，RMSE)[2]。

采用3種降噪算法對公開數(shù)據(jù)集和自采集的心音信號進行降噪，來驗證各類算法的降噪效果。其中，公開數(shù)據(jù)集中的數(shù)據(jù)相對比較干凈，故對該類數(shù)據(jù)集中的心音進行降噪之前，將高斯白噪聲以10 dB的SNR分別疊加到所選的心音中，來構建含噪心音。而部分自采集數(shù)據(jù)集中本身含有大量的噪聲，故不對其進行加噪處理。各類算法對不同類型心音的降噪結果如圖5所示。

圖5 3種算法對不同類型心音的降噪結果

圖5(a)、圖5(b)和圖5(c)分別是上述3種算法對正常/異常心音降噪后的時域波形圖和頻譜圖的對比結果。其中，圖5(a)和圖5(b)分別使用的是公開數(shù)據(jù)集中的正常心音和異常心音，而圖5(c)使用的是自采集心音數(shù)據(jù)集的正常心音。為了便于對比分析，圖中用橢圓進行標記。從圖5(a)和圖5(c)的時域波形圖中可以看出，基于OMLSA-小波算法因?qū)⒉糠诌吘売杏眯盘柡图捌湮⑷醯碾s音當成噪聲信號進行估計造成后期估計的基礎心音成分有所丟失，而基于小波閾值算法和本文算法降噪后的波形與原始心音波形差別不大；然而，從圖5(a)和圖5(c)的頻譜圖中可以看出，基于小波閾值算法能通過頻帶去除高頻和低頻噪聲使非基礎心音部分殘留的噪聲較少，但基礎心音部分卻因為丟失的高頻和低頻分量中包含部分有用心音成分而造成了FHS的失真。另外從圖5(b)的時域波形圖和頻譜圖中也可以看出，基于小波閾值算法和基于OMLSA-小波算法在對噪聲去除的同時將噪聲和基礎心音信號混疊構成的異常心音信號成分一同丟失，導致出現(xiàn)較為明顯的心音失真，且基于OMLSA-小波算法降噪后的心音在非基礎心音部分仍存在部分中頻噪聲。相比上述兩種算法，本文算法可使心音信號在減小FHS失真的前提下，更為有效地抑制噪聲，重構出一個更為平滑自然的頻譜，更加符合后續(xù)分類處理的需求。綜上時頻圖的對比分析來看，提出算法較基于小波閾值的算法和基于OMLSA-小波的算法更優(yōu)。

3.4 各類算法降噪結果的指標評測

此外，本文對3種算法的降噪性能進行了客觀的對比，評測所用心音包括正常心音(N)與異常心音。其中，異常心音為二尖瓣反流(mitral regurgitation，MR)心音、二尖瓣脫垂(mitral valve prolapse，MVP)心音、二尖瓣狹窄(mitral stenosis，MS)心音和主動脈狹窄(aortic stenosis，AS)心音4類[23]。3種算法對不同信噪比的正常與異常心音降噪后輸出SNR的統(tǒng)計平均值見表2。輸出SNR值越大，表示算法的降噪性能越好[2，6，7]。

表2 對不同信噪比的含噪心音降噪后的SNR

從表2中可以看出，基于OMLSA-小波算法和基于小波閾值算法要明顯低于提出算法的輸出SNR；原因是前兩種算法降噪后的心音信號殘留了較多的中頻噪聲，且存在FHS失真，尤其是異常心音更為明顯。而提出算法對SNR為0 dB的MVP所得SNR值要低于基于OMLSA-小波算法，這是源于提出算法中的IMCRA為了防止噪聲過估計，使部分噪聲殘留，而MVP相對其它類型的心音殘留的噪聲較多，但整體上對多種類型的心音有著較好的降噪效果，且能在保證心音信號不失真的前提下，對噪聲最大程度的抑制。

鑒于SNR的部分不良表現(xiàn)，故進一步通過RMSE來檢驗3種算法在降噪中對正常心音與異常心音的失真與降噪平衡度。RMSE越小，表示降噪后的心音越接近于真實干凈心音，即降噪性能越好[2，6，7]，反之降噪性能越差。從數(shù)據(jù)中抽選N、MVP、MS和MR這4種類型的心音，利用3種算法分別對不同信噪比下的4類心音進行降噪處理，其RMSE對比結果[2]如圖6所示。

圖6 不同信噪比的4類心音降噪后的RMSE對比曲線

從圖6中可以看出，隨著噪聲水平的降低，提出算法的RMSE也逐漸減小。圖6(a)顯示，SNR為20 dB時，提出算法降噪后的RMSE為0.0124，說明降噪后的信號與原始信號之間相似度較高，即降噪效果較好；SNR為0 dB時，提出算法降噪后的RMSE為0.0399，降噪后的信號與原始信號之間相似度降低，降噪性能下降；這是因為低SNR下，為了確保FHS無失真引起了部分噪聲殘留。從圖6(a)、圖6(b)、圖6(c)和圖6(d)可以看出，除圖6(c)中因防止噪聲過估計造成MVP信號內(nèi)殘留部分噪聲，從而使提出算法在0 dB～5 dB間的RMSE值略高于基于OMLSA-小波的算法之外，其余情況下，提出算法對N、MR、MVP和MS這4種類型心音進行降噪所得RMSE值均低于對比算法，說明提出算法降噪后的信號更接近于原始干凈心音信號，能較好地保留FHS成分，更為適用與不同噪聲水平下的多種類型心音降噪。

4 結束語

本文根據(jù)心音信號的基本特征，設計HSAD方法判決當前心音幀是否為基礎心音幀，并分別采用遞歸平滑和IMCRA動態(tài)對非基礎心音和基礎心音的噪聲功率進行估計與更新。該方法一方面結合HSAD，來減小收縮期和舒張期雜音對噪聲估計的影響，彌補IMCRA的不足，減小雜音和噪聲的殘留；另一方面采用非因果信噪比替代因果信噪比，避免了由非基礎心音向基礎心音過渡時先/后驗信噪比偏小的問題，減小了基礎心音成分的失真。