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

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

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

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

0 引言

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

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

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

1 特征提取

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

1.1 預(yù)處理

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

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

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

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

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

式中：o為重疊長度；wlen為幀長度；c為幀位移量。

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

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

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

1.2 短時(shí)能量

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

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

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

1.3 HHT-MFCC

希爾伯特黃變換是一種對非線性、非平穩(wěn)信號的自適應(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ù)處理后的信號分解成若干固有模態(tài)函數(shù)（Intrinsic Mode Function，IMF）和殘余分量；Hilbert 譜分析是對每個IMF 進(jìn)行分析得到相應(yīng)的Hilbert譜，反映了信號的瞬時(shí)頻率，是傅里葉展開的一般化［11］，所有IMF 的Hilbert 譜之和即為原始信號的Hilbert 譜。HHT 具有自適應(yīng)和能體現(xiàn)信號非平穩(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ù)雜信號，都可以通過EMD得到若干IMF，但每個IMF需要滿足兩個條件［12］：

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

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

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

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

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

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

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

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

圖3 預(yù)處理后健康人呼吸聲信號經(jīng)EMD分解得到的從高頻到低頻的前五個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 后，需要對分量進(jìn)行Hilbert 譜分析，但主要信息多集中于高頻分量，為減少計(jì)算量，選取前五個IMF 進(jìn)行分析。同樣地，由于殘余信號屬于低頻信號，而高頻信號的信息較為重要，一般情況會在譜分析中省略殘余信號。因此經(jīng)Hilbert譜分析的信號為：

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

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

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

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

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

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

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

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

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

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

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

HHT-MFCC特征提取公式為：

其中：

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

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

1.4 特征融合

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

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

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

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

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

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

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

4 結(jié)語

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