基于活動(dòng)性檢測(cè)動(dòng)態(tài)估計(jì)噪聲的心音降噪算法

許春冬，辛鵬麗，閔 源+，應(yīng)冬文，2，周 靜，李海兵

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

0 引 言

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

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

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

1 算法框架

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

圖1 心音降噪算法流程

2 算法原理

2.1 MMSE-LSA算法原理

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

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

(1)

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

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

(2)

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

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

(3)

其中，l和k分別代表幀數(shù)和頻點(diǎn)數(shù)的索引，Y(k，l)、S(k，l)和D(k，l) 分別是y(n)、s(n) 和d(n) 進(jìn)行STFT后的結(jié)果，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) 為先驗(yàn)信噪比，γ1(k，l)=|Y(k，l)|2/λD(k，l) 為后驗(yàn)信噪比，λS(k，l) 為心音功率譜，λD(k，l) 為噪聲功率譜。

2.2 非因果先驗(yàn)信噪比估計(jì)

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

(7)

(8)

(9)

2.3 噪聲估計(jì)

2.3.1 IMCRA噪聲估計(jì)算法

(10)

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

(11)

(12)

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

(13)

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

(14)

(15)

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

(16)

其中，γ0和ζ0為常數(shù)閾值參數(shù)，基礎(chǔ)心音存在時(shí)I(k，l) 為1，基礎(chǔ)心音不存在時(shí)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的噪聲功率估計(jì)

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

圖2 心音信號(hào)時(shí)域波形

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

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

活動(dòng)性檢測(cè)是聲學(xué)信號(hào)處理領(lǐng)域里較為常用的處理手段，且類別繁多[13，20]。為了在有效檢測(cè)心音活動(dòng)性的前提下進(jìn)一步考慮低的算法復(fù)雜度，本文采用對(duì)數(shù)頻譜距離(logarithmic amplitude spectrum distance，LASD)來(lái)設(shè)計(jì)適合心音信號(hào)的HSAD方法，用于判別基礎(chǔ)心音幀和非基礎(chǔ)心音幀[21]。基于前幾幀為非基礎(chǔ)心音幀的假設(shè)，可通過(guò)式(28)將LASD計(jì)算如下

(28)

(29)

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

(30)

(31)

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

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

(32)

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

表1 參數(shù)的取值

3 實(shí)驗(yàn)仿真與結(jié)果分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

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

3.2 HSAD-IMCRA噪聲估計(jì)性能分析

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

圖4 噪聲估計(jì)算法相關(guān)指標(biāo)結(jié)果對(duì)比

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

3.3 各類算法降噪結(jié)果的時(shí)頻圖對(duì)比

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

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

圖5 3種算法對(duì)不同類型心音的降噪結(jié)果

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

3.4 各類算法降噪結(jié)果的指標(biāo)評(píng)測(cè)

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

表2 對(duì)不同信噪比的含噪心音降噪后的SNR

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

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

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

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

4 結(jié)束語(yǔ)

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