基于HHT-MFCC和短時(shí)能量的慢性阻塞性肺病患者呼吸聲識(shí)別

常 崢，羅 萍，楊 波，張曉曉

（重慶郵電大學(xué)自動(dòng)化學(xué)院，重慶 400065）

（*通信作者電子郵箱1213884879@qq.com）

0 引言

呼吸系統(tǒng)是身體與外界進(jìn)行氣體交換的器官總稱，肺部是系統(tǒng)中重要的組成部分，負(fù)責(zé)氧氣輸送和二氧化碳排出［1］。肺作為身體的重要器官之一，人們從未停止對(duì)肺部相關(guān)疾病的研究。隨著年齡增長(zhǎng)和環(huán)境惡化，60 歲以上老年人的肺部情況較為惡劣，尤其是慢性阻塞性肺?。–hronic Obstructive Pulmonary Disease，COPD），這種疾病具有發(fā)病周期長(zhǎng)且難以治療的特點(diǎn)，是老年人主要醫(yī)療經(jīng)濟(jì)負(fù)擔(dān)和死亡原因。

聲音是信息的主要載體，人們可以通過聲音來判斷事物的某種狀態(tài)，通過呼吸聲信號(hào)診斷COPD 能夠解決患者到醫(yī)院就診所帶來的診斷周期長(zhǎng)、時(shí)間成本高等問題，因此在國內(nèi)外出現(xiàn)很多關(guān)于COPD 呼吸聲信號(hào)方面的研究。Jaber 等［2］利用改進(jìn)的隨機(jī)森林分類器對(duì)肺部聲音信號(hào)進(jìn)行肺部疾病有效診斷，其中包括COPD、哮喘等肺部疾病；Azam 等［3］采用CEEMD（Complete Ensemble Empirical Mode Decomposition）并結(jié)合支持向量機(jī)（Support Vector Machine，SVM）的方法區(qū)分COPD、哮喘和支氣管炎等疾?。籔ramono 等［4］通過計(jì)算復(fù)雜度較低的線性分類器區(qū)分健康、哮喘和COPD 呼吸聲；Islam 等［5］通過提取功率譜密度子帶中特征向量并輸入到人工神經(jīng)網(wǎng)絡(luò)（Artificial Neural Network，ANN）分類器實(shí)現(xiàn)識(shí)別健康人、哮喘和COPD 患者；Altan 等［6］采用集合經(jīng)驗(yàn)?zāi)B(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）和深度信念網(wǎng)絡(luò)的方法區(qū)別COPD 患者和健康人呼吸聲，EEMD 解決了EMD 中模態(tài)混疊問題；Haider 等［7］提取中值頻率和線性預(yù)測(cè)編碼特征并結(jié)合SVM，實(shí)現(xiàn)結(jié)合呼吸聲測(cè)量數(shù)據(jù)和呼吸聲信號(hào)特征來區(qū)分COPD 患者呼吸聲。這些研究用不同方法通過分析呼吸聲信號(hào)實(shí)現(xiàn)分類識(shí)別，但特征信號(hào)比較單一，基本采用從頻域中提取特征向量作為識(shí)別診斷的依據(jù)，沒有考慮到時(shí)域中信號(hào)的部分信息，達(dá)不到理想的識(shí)別效果。

本文對(duì)采集的呼吸聲信號(hào)進(jìn)行預(yù)處理，提取信號(hào)的短時(shí)能量，同時(shí)將希爾伯特黃變換（Hilbert-Huang Transform，HHT）代替梅爾頻率倒譜系數(shù)（Mel-Frequency Cepstral Coefficients，MFCC）中的快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT），對(duì)信號(hào)進(jìn)行特征提取，再將兩者進(jìn)行融合得到新特征向量。在Matlab R2018a 環(huán)境下，本文融合算法與MFCC 和HHT-MFCC 利用SVM 進(jìn)行實(shí)驗(yàn)，通過實(shí)驗(yàn)比較證明了該融合算法的合理性和優(yōu)越性，并完成了COPD 患者和健康人呼吸聲的分類識(shí)別。

1 特征提取

特征提取的過程就是用較少維數(shù)表征音頻信號(hào)特征的過程，因此COPD 識(shí)別率的高低很大程度上取決于對(duì)呼吸聲信號(hào)提取的特征參數(shù)是否準(zhǔn)確。特征參數(shù)既要把呼吸聲信號(hào)中的特性盡可能提取出來，還需要將COPD 患者和健康人呼吸聲的特征參數(shù)有效地區(qū)分開，從而提高識(shí)別率。在特征提取之前，一般都需要進(jìn)行預(yù)處理操作，減少無用信息和噪聲的干擾，提高特征提取時(shí)的準(zhǔn)確性。

1.1 預(yù)處理

音頻信號(hào)是時(shí)變信號(hào)，但具有短時(shí)平穩(wěn)性的特點(diǎn)，即在短時(shí)范圍內(nèi)可以作為穩(wěn)態(tài)信號(hào)處理。通常將一段音頻信號(hào)分為若干等長(zhǎng)幀，幀的長(zhǎng)度一般為10～30 ms。在音頻信號(hào)識(shí)別系統(tǒng)中的研究與運(yùn)算都是基于信號(hào)短時(shí)分析基礎(chǔ)上實(shí)現(xiàn)的。為更好地提取特征參數(shù)，對(duì)音頻信號(hào)進(jìn)行模數(shù)轉(zhuǎn)換（Analogueto-Digital Conversion，ADC）、預(yù)加重和分幀加窗等預(yù)處理操作。

1）模數(shù)轉(zhuǎn)換，即將獲取到的模擬音頻信號(hào)轉(zhuǎn)換為易處理的數(shù)字信號(hào)。

2）預(yù)加重。由于呼吸聲的頻率為60～1 600 Hz［8］，而頻率高于800 Hz 時(shí)會(huì)有6 dB 的衰減，因此通過預(yù)加重，提高呼吸聲的高頻部分，其傳遞函數(shù)為：

其中：μ∈[0.9，1]，一般情況下取μ=0.937 5。

3）分幀加窗。由于音頻信號(hào)在短時(shí)是平穩(wěn)的［9］，便于特征分析，需要對(duì)信號(hào)分幀，同時(shí)為保證每相鄰兩幀可以平滑過渡，需要對(duì)幀信號(hào)進(jìn)行疊加，再對(duì)信號(hào)加入窗函數(shù)。幀移公式如式（2）所示，加窗信號(hào)yn(m)如式（3）所示。

式中：o為重疊長(zhǎng)度；wlen為幀長(zhǎng)度；c為幀位移量。

式中：w(k)為一個(gè)窗函數(shù)，本文采用漢明窗；x(m)為音頻信號(hào)。

圖1（a）、（b）分別為健康人和COPD 患者呼吸聲的原始信號(hào)（以下呼吸聲信號(hào)均采用一個(gè)呼吸周期內(nèi)的呼吸聲）；圖1（c）、（d）分別是上述兩種呼吸聲經(jīng)過預(yù)處理后的信號(hào)。

圖1 健康人和COPD患者呼吸聲的原始信號(hào)與預(yù)處理后信號(hào)Fig.1 Original and preprocessed respiratory sound signals of healthy people and COPD patients

1.2 短時(shí)能量

短時(shí)能量是對(duì)音頻信號(hào)的時(shí)域分析，能量大小能夠反映出一段信號(hào)的特征。n時(shí)刻音頻信號(hào)的短時(shí)能量En為：

式中：w(k)為窗函數(shù)，N為窗長(zhǎng)，窗函數(shù)采用漢明窗，窗長(zhǎng)取200個(gè)采樣點(diǎn)。

圖2（a）、（b）分別是經(jīng)過預(yù)處理后的健康人和COPD 患者呼吸聲的短時(shí)能量。

1.3 HHT-MFCC

希爾伯特黃變換是一種對(duì)非線性、非平穩(wěn)信號(hào)的自適應(yīng)性時(shí)頻分析方法［10］，主要由經(jīng)驗(yàn)?zāi)B(tài)分析（Empirical Mode Decomposition，EMD）和希爾伯特譜分析（Hilbert Spectrum Analysis，HSA）兩部分構(gòu)成。EMD 是將預(yù)處理后的信號(hào)分解成若干固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF）和殘余分量；Hilbert 譜分析是對(duì)每個(gè)IMF 進(jìn)行分析得到相應(yīng)的Hilbert譜，反映了信號(hào)的瞬時(shí)頻率，是傅里葉展開的一般化［11］，所有IMF 的Hilbert 譜之和即為原始信號(hào)的Hilbert 譜。HHT 具有自適應(yīng)和能體現(xiàn)信號(hào)非平穩(wěn)特性的優(yōu)點(diǎn)，克服了FFT 必須選取特定基函數(shù)以及只能全局分析的缺陷。

圖2 經(jīng)過預(yù)處理后的健康人和COPD患者呼吸聲的短時(shí)能量Fig.2 Short-term energies of respiratory sound signals of healthy people and COPD patients after pretreatment

無論是線性的、非線性的、平穩(wěn)的還是非平穩(wěn)的復(fù)雜信號(hào)，都可以通過EMD得到若干IMF，但每個(gè)IMF需要滿足兩個(gè)條件［12］：

1）在整個(gè)數(shù)據(jù)序列中，極值點(diǎn)和過零點(diǎn)的數(shù)目必須相等或至多相差一個(gè)；

2）在任意時(shí)刻，局部極大值包絡(luò)與局部極小值包絡(luò)定義的局部平均包絡(luò)線必須為零。

對(duì)于一個(gè)經(jīng)過預(yù)處理后的呼吸聲信號(hào)x(t)，其EMD 步驟如下：

1）求出信號(hào)的局部極大值和局部極小值，通過三次樣條插值得到由局部極小值構(gòu)成的下包絡(luò)線xlow(t)和由局部極大值構(gòu)成的上包絡(luò)線xhigh(t)，上下包絡(luò)線通過公式m1=計(jì)算局部平均包絡(luò)線m1。

2）求出原信號(hào)x(t)與m1的差值，即h10=x(t) -m10。h10如果滿足IMF 條件，則是一個(gè)IMF 分量；如果不滿足，將其作為新的原信號(hào)，重復(fù)上述兩個(gè)步驟，得到h11=h10-m11。每得到新差值，都要判斷是否滿足條件，若仍不滿足則繼續(xù)重復(fù)步驟，直到重復(fù)k次 時(shí)h1k=k1(k-1)-m1k滿足IMF篩選停止原則：

式中：SD的典型值在0.2 到0.3 之間。第一個(gè)IMF 分量h1k記為c1。

3）由信號(hào)x(t)與c1可得到r1=x(t) -c1，將r1作為新信號(hào)重復(fù)上述三個(gè)步驟，即可得到第二個(gè)IMF，記作c2。重復(fù)步驟3），直到殘余信號(hào)rn為單調(diào)曲線或第n個(gè)IMF 分量小于設(shè)定值時(shí)，則停止EMD 操作，這樣就可以得到n個(gè)IMF 分量c1，c1，…，cn。因此得到原信號(hào)x(t)與IMF 分量和殘余信號(hào)的關(guān)系：

圖3是預(yù)處理后健康人呼吸聲信號(hào)經(jīng)EMD分解得到的從高頻到低頻的前五個(gè)IMF 分量，橫坐標(biāo)為時(shí)間，縱坐標(biāo)為幅值。

圖3 預(yù)處理后健康人呼吸聲信號(hào)經(jīng)EMD分解得到的從高頻到低頻的前五個(gè)IMF分量Fig.3 First five IMF components from high to low frequency obtained by EMD of healthy people respiratory sound signal after pretreatment

通過EMD 得到IMF 后，需要對(duì)分量進(jìn)行Hilbert 譜分析，但主要信息多集中于高頻分量，為減少計(jì)算量，選取前五個(gè)IMF 進(jìn)行分析。同樣地，由于殘余信號(hào)屬于低頻信號(hào)，而高頻信號(hào)的信息較為重要，一般情況會(huì)在譜分析中省略殘余信號(hào)。因此經(jīng)Hilbert譜分析的信號(hào)為：

式（7）的傅里葉變換為：

式（7）表示前五個(gè)IMF在時(shí)域中的頻率和幅值，若在時(shí)間和頻率平面中表示IMF 的幅值稱為Hilbert 譜，用H(w，t)表示，表達(dá)式為：

得出Hilbert譜后，便可定義Hilbert邊際譜H(w)：

Hilbert譜是對(duì)每個(gè)IMF 在時(shí)頻平面中求局部幅值并將每一時(shí)刻所有幅值加權(quán)，而Hilbert 邊際譜是對(duì)整個(gè)信號(hào)的幅值的累積，因此它是對(duì)信號(hào)總幅值或總能量的度量。

在Hilbert 邊際譜中，某一頻率的幅值是每個(gè)IMF 幅值的累加，某一頻率存在能量意味著在Hilbert 譜中時(shí)頻平面上可能存在該頻率的振動(dòng)波，而在傅里葉頻譜中某一頻率處的幅值意味著信號(hào)中存在一個(gè)該頻率的三角函數(shù)，某一頻率存在能量意味著時(shí)間軸中存在一個(gè)正弦波或余弦波。這意味著FFT 只能反映整個(gè)時(shí)間軸中固定的基函數(shù)，而Hilbert 邊際譜能更好地體現(xiàn)信號(hào)的非平穩(wěn)性和局部特性。因此，用HHT 代替依據(jù)人耳聽覺特性分析的MFCC 中FFT 可以提取更好的特征參數(shù)。

圖4為HHT-MFCC的提取流程詳細(xì)步驟為：

1）對(duì)采集的信號(hào)x(t)進(jìn)行預(yù)處理，包括預(yù)加重、分幀加窗等操作，使信號(hào)x(t)變?yōu)檩^為平滑的信號(hào)x(m)并保持其短時(shí)穩(wěn)定性。

2）將處理后的信號(hào)x(m)通過EMD 得到若干IMF，本文取前五個(gè)IMF進(jìn)行HHT變換。

3）將變換后的信號(hào)幅值累加，得到Hilbert 邊際譜H(w)，再取平方計(jì)算出每幀Hilbert邊際譜的能量。

4）通過Mel濾波器，將每幀譜能量轉(zhuǎn)換成在Mel頻域中的能量，進(jìn)而求出n維特征向量。

5）對(duì)特征向量取對(duì)數(shù)和離散余弦變換（Discrete Cosine Transform，DCT），便可得到n維HHT-MFCC特征向量。

HHT-MFCC特征提取公式為：

其中：

式中：n是DCT 后的譜線；m是指第m個(gè)Mel 濾波器（共有M個(gè)）；k是頻率中的第k條譜線；N是信號(hào)x(t)的長(zhǎng)度；Hm(k)是Mel濾波器的頻域響應(yīng)；xi(m)是經(jīng)預(yù)處理后第i幀音頻信號(hào)。

圖4 HHT-MFCC提取流程Fig.4 HHT-MFCC extraction flowchart

1.4 特征融合

呼吸聲信號(hào)是一類較為復(fù)雜的信號(hào)，它包含了呼吸聲音、心跳生理音和環(huán)境噪聲。在經(jīng)過預(yù)處理后，雖然噪聲干擾降低了，但僅從頻域分析提取信號(hào)特征不足以多角度提取信號(hào)特征向量，導(dǎo)致部分重要特征丟失，識(shí)別率不高。因此，本文采用HHT-MFCC 與短時(shí)能量融合算法得到呼吸聲信號(hào)的新特征矢量，從兩個(gè)不同角度分析呼吸聲信號(hào)能更好地反映其特征。同時(shí)，由于計(jì)算短時(shí)能量的運(yùn)算復(fù)雜度低，并不會(huì)給融合后的特征帶來明顯特征矢量維數(shù)。

特征融合就是將兩種或兩種以上音頻信號(hào)的特征向量通過首尾相連的方式進(jìn)行組合。通過HHT-MFCC 特征提取算法，可以得到呼吸聲信號(hào)的特征向量為：

式中：FNM為第N幀、第M維HHT-MFCC；T1為N行M列的矩陣。該段呼吸聲信號(hào)的短時(shí)能量特征向量為：

式中：N為幀數(shù)，T2為N行1 列矩陣。混合之前分別把每個(gè)特征矩陣進(jìn)行重組，得到矩陣為：

2 SVM

SVM 是以統(tǒng)計(jì)學(xué)習(xí)理論為基礎(chǔ)解決機(jī)器學(xué)習(xí)問題的方法，是從線性可分問題下的最優(yōu)分類面發(fā)展而來的，解決了小樣本、非線性和“過學(xué)習(xí)”等難題［13］。同時(shí)，它能根據(jù)有限的樣本在模型的復(fù)雜性和學(xué)習(xí)能力之間尋求最佳折中。本文通過短時(shí)能量和HHT-MFCC 的融合算法提取特征向量，采用SVM進(jìn)行分類識(shí)別，并與MFCC和HHT-MFCC算法進(jìn)行比較。

假設(shè)訓(xùn)練樣本集為Z={(x1，y1)，(x2，y2)，…，(xi，yi)}，其中xi∈Rn，yi∈{1，-1}，i=1，2，…，n，xi為第i個(gè)訓(xùn)練樣本特征向量，yi是第i個(gè)訓(xùn)練樣本分類標(biāo)簽［14］。具體訓(xùn)練步驟如下：

1）給定訓(xùn)練樣本集滿足式（17）且令Q(λ)為最大值：

式中：λi為拉格朗日乘子，當(dāng)Q(λ)為最大值，λ取值為λ*。

3）根據(jù)w*和b*可得到超平面f(x)，表達(dá)式為：

式中，K(xi，yi)為核函數(shù)，是SVM 的核心算法。核函數(shù)能將二維線性不可分的問題轉(zhuǎn)化為高維線性可分問題，本文采用的核函數(shù)是徑向基核函數(shù)（Radial Basis Function，RBF），表達(dá)式為：

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

本文采用的數(shù)據(jù)集是由葡萄牙和希臘的兩個(gè)研究團(tuán)隊(duì)為生物醫(yī)學(xué)健康信息學(xué)國際聯(lián)盟（International Confederation on Biomedical Health Informatics，ICBHI）組織的挑戰(zhàn)賽而制作的［15］，其中包含了健康人和COPD 患者呼吸聲。本文選取兩類呼吸聲來自60 名60～80 周歲的測(cè)試人員，共計(jì)429 條呼吸聲樣本，包括187 條健康人呼吸聲樣本和242 條COPD 患者呼吸聲樣本。為保證兩類數(shù)據(jù)保持一致，呼吸聲信號(hào)均選取8 000 Hz 采樣率、幀長(zhǎng)取20 ms、幀移取10 ms，窗函數(shù)采用漢明窗，窗長(zhǎng)取200 個(gè)采樣點(diǎn)。由于每段呼吸聲均超過20 s，為方便計(jì)算，截取2～17 s 作為呼吸聲信號(hào)。兩類疾病各隨機(jī)挑選出20、40、60、80 和100 個(gè)訓(xùn)練樣本，其余作為測(cè)試樣本使用。在Matlab R2018a 環(huán)境下，分別對(duì)MFCC、HHT-MFCC 和HHT-MFCC+Energy 三種算法進(jìn)行比較，并對(duì)不同特征矢量的影響進(jìn)行判斷，表1和表2分別表示針對(duì)不同維數(shù)的三種特征提取算法的識(shí)別性能，維數(shù)表示HHT-MFCC 的特征矢量也是Mel 濾波器的個(gè)數(shù)，若為融合特征提取方法，則計(jì)算出倒譜系數(shù)后還需要加上一維短時(shí)能量特征。

表1 不同特征提取方法識(shí)別率（24維）Tab.1 Recognition rates of different feature extraction methods（24 dimensions）

表2 不同特征提取方法識(shí)別率（12維）Tab.2 Recognition rates of different feature extraction methods（12 dimensions）

從實(shí)驗(yàn)結(jié)果來看，當(dāng)特征提取的維數(shù)相同時(shí)，HHTMFCC+Energy 的融合算法識(shí)別率最高，這是因?yàn)樵撍惴饶軓亩喾矫娅@取信號(hào)特征，同時(shí)包含時(shí)域和頻域中的特點(diǎn)，使得特征描述更完整，也能改善MFCC 中只能全局分析的缺陷，因而識(shí)別率更高；對(duì)于融合算法，維數(shù)增加意味著獲取到更多信號(hào)細(xì)節(jié)，兩類聲音能更好區(qū)分，因此24 維的特征提取方法要比12維的識(shí)別率高。同樣地，相較于文獻(xiàn)［6］中93.67%的識(shí)別率和文獻(xiàn)［7］中83.6%的識(shí)別率，本文提出的方法分別提升了4.13 個(gè)百分點(diǎn)和14.2 個(gè)百分點(diǎn)。實(shí)驗(yàn)結(jié)果表明，HHTMFCC+Energy 特征提取算法結(jié)合SVM 的識(shí)別率最高，對(duì)COPD 患者呼吸聲的識(shí)別更準(zhǔn)確，特別是當(dāng)特征矢量為24，且樣本數(shù)量為100時(shí)，對(duì)COPD患者和健康人呼吸聲的識(shí)別率分別能達(dá)到97.5%和98.1%，基本實(shí)現(xiàn)了利用呼吸聲信號(hào)診斷被測(cè)人員是否患有COPD的目標(biāo)。

4 結(jié)語

為了提高呼吸聲信號(hào)識(shí)別率，增加COPD 疾病的識(shí)別率，本文提出了一種融合短時(shí)能量的HHT-MFCC 特征提取算法HHT-MFCC+Energy。該算法能準(zhǔn)確地反映呼吸聲信號(hào)非平穩(wěn)特性和局部特征，也解決了僅從頻域中提取特征的方法不能較為全面捕獲音頻信號(hào)特征的問題，基本實(shí)現(xiàn)了COPD 患者與健康人呼吸聲的識(shí)別。本文預(yù)處理部分未能徹底分離呼吸聲和心跳生理音，導(dǎo)致特征提取的過程中仍存在部分噪聲，若能去除心跳聲的干擾，則可以進(jìn)一步提高特征向量表述信號(hào)的準(zhǔn)確性，實(shí)現(xiàn)更好的識(shí)別效果。