咳嗽識別在臨床上具有重要的診斷指導意義。針對咳嗽頻譜能量的分布特點,本文提出了一種新的梅爾(Mel)頻率倒譜系數(MFCC)提取方法。將咳嗽頻譜劃分為若干個頻段,采用主元分析方法計算各頻段的能量強度系數,根據強度系數的插值曲線分配濾波器個數,設計Mel刻度上非均勻分布的濾波器組進行MFCC特征提取。基于隱馬爾可夫模型(HMM)的咳嗽識別實驗表明,該方法可以有效改善咳嗽識別的效果。
引用本文： 朱春媚, 莫鴻強, 田聯房, 鄭則廣. 基于主元分析法和非均勻濾波器組的咳嗽信號特征提取. 生物醫學工程學雜志, 2015, 32(4): 746-750. doi: 10.7507/1001-5515.20150136 復制
引言
咳嗽是呼吸系統疾病常見的癥狀，咳嗽信號的記錄和自動識別在臨床上具有重要的診斷指導意義。咳嗽記錄主要通過硬件實現，目前已有多種類型的便攜式咳嗽監測儀，能有效實現對咳嗽的長時間記錄[1-3]。咳嗽記錄結果的分析目前仍主要靠人工分析，而咳嗽信號的自動識別是咳嗽監測儀推廣應用的關鍵，有效的特征提取方法則是實現高精度識別的基礎。
過去幾年里，研究者主要參照語音識別的方法，以傳統梅爾(Mel)頻率倒譜系數(Mel frequency cepstrum coefficients,MFCC)為特征參數，采用動態時間規整、隱馬爾可夫模型(hidden Markov model,HMM)、高斯混合模型或神經網絡為分類器來實現咳嗽識別[4-8]。咳嗽信號與語音信號相比，在時域和頻域特性都有其自身的特點，在特征提取或識別的時候充分考慮咳嗽信號的特性將會使算法更加有效。
傳統MFCC的計算首先用快速傅立葉變換將時域信號轉化成頻域，之后對其對數能量譜依照Mel刻度上均勻分布的三角濾波器組進行卷積，最后對各個濾波器的輸出構成的向量進行離散余弦變換，取前N個系數[9]。雖然Mel刻度反映了人耳對不同信號的感知特點，但事實上，咳嗽頻譜能量在Mel刻度上并非均勻分布，即并不隨頻率的增加而遞減，而是在某些頻段出現峰值[10]，例如慢性支氣管炎患者咳嗽的頻率集中在500、700和1 200 Hz附近[11]，又如圖 1所示是一個典型干性咳嗽的時域波形和語譜圖,其能量峰值出現在[300,1 000]Hz和[2 000，3 000] Hz范圍。這種情況下采用Mel刻度上均勻分布的濾波器組，會出現頻譜能量低的頻段反而放置了較多濾波器的問題，不能最大程度地獲得咳嗽特征。基于這一點，本文將整個頻域范圍劃分為若干個頻段，采用主元分析法分析咳嗽信號的頻譜能量分布，按各頻段的能量強度系數分配濾波器，設計Mel刻度上的非均勻濾波器組，在高能量頻段放置較多濾波器，在低能量頻段放置較少濾波器，改進傳統MFCC參數的提取過程。
圖1
咳嗽的時域波形和語譜圖
Figure1.
Time domain waveform and spectrogram of a typical cough
1 數據來源
實驗樣本來自廣州醫學院第一附屬醫院病房，使用DELL原裝機的聲卡信號采集設備進行單聲道采樣,16位數字量化。盡管文獻[11]和[12]都報道觀察到了5 000～6 000 Hz的咳嗽分量，但3 500 Hz以上頻段的咳嗽分量很少[10]，忽略這部分高頻成分對咳嗽特征提取影響不大，因此為了減少計算量，本文采樣頻率取8 000 Hz，對應信號的處理范圍為0～4 000 Hz。選取不同性別、年齡和呼吸道感染情況的患者進行連續、長時間的日常錄音，人工分割成單個咳嗽信號作為咳嗽樣本，并隨機截取時長0.2～3 s的非咳嗽信號作為非咳嗽樣本。
2 非均勻濾波器組設計
2.1 總體流程
非均勻濾波器組設計的基本思想是根據咳嗽的頻譜能量分布來分配濾波器個數，在高能量頻段放置較多濾波器，在低能量頻段放置較少濾波器。首先將整個頻域范圍劃分為若干個頻段，采用主元分析法對大量的咳嗽樣本進行分析，得到咳嗽信號的頻譜能量分布系數，然后插值得到頻譜能量曲線，對能量曲線的積分進行均勻劃分，分界點即為濾波器的中心點。非均勻濾波器組的設計總體流程如圖 2所示。
圖2
非均勻濾波器組的總體設計流程
Figure2.
Overall design process of non-uniform filter-bank
2.2 主元分析法求頻譜能量分布
按幀長240(30 ms)、幀移80(10 ms)對咳嗽信號進行分幀，將每幀[0,4 000] Hz的頻率范圍劃分為n個頻段，每個頻段與Mel刻度上n個濾波器的均勻濾波器組的通帶范圍對應，則每幀信號由這n個頻段的頻譜按一定方式組合而成，記第i幀語音信號為fi，則fi可表示為n個頻段能量的線性組合:
| ${f_i} = \sum\limits_{k = 1}^n {{\omega _{ik}}{x_{ik}},} $ |
其中xik為第i幀信號第k個頻段的頻譜能量和；ωik為第i幀信號第k個頻段的能量系數，它反映了第k個頻段的能量在該幀信號頻譜中所占的比例。對多幀的咳嗽信號可得到能量矩陣A：
| $A = \left( {\begin{array}{*{20}{c}} {{x_{11}}}&{{x_{12}}}& \cdots &{{x_{1n}}}\\ {{x_{21}}}&{{x_{22}}}& \cdots &{{x_{2n}}}\\ \cdots & \cdots & \cdots & \cdots \\ {{x_{m1}}}&{{x_{m2}}}& \cdots &{{x_{mn}}} \end{array}} \right)$ |
其中m是總的幀數，n為頻段數。主元分析法能夠找出數據中最重要的元素[13]，對A進行主元分析，找到主元對應的頻段系數即為統計意義上的頻段能量系數。具體步驟如下。
步驟1：計算能量矩陣A。對m幀咳嗽信號分別求快速傅里葉變換，將頻段在Mel刻度上平均劃分為n個頻段，計算每幀每個頻段的頻譜能量，得到m*n的能量矩陣A。
步驟2：計算協方差矩陣ATA，得到n*n的方陣。
步驟3：計算協方差矩陣ATA的特征值λ和特征向量P，滿足
| ${P^T}{A^T}AP = \Lambda ,$ |
| $A = \left( {\begin{array}{*{20}{c}} {\sqrt {{\lambda _1}} }&0& \cdots &0\\ 0&{\sqrt {{\lambda _2}} }& \cdots &0\\ \cdots & \cdots & \cdots & \cdots \\ 0&0& \cdots &{\sqrt {{\lambda _n}} } \end{array}} \right),$ |
其中λ按從大到小排序，按85%的累計貢獻率，得到主元個數為L。
步驟4：求各頻段的能量比重，第k個頻段的能量比重為
| ${\alpha _k} = \sum\limits_{j = 1}^L {{\lambda _j}{P_{kj}}} ,k = 1,2, \cdots ,n$ |
本文選取了92個典型的咳嗽信號，共4 015幀，按Mel刻度上均勻濾波器組的帶寬范圍將0～4 000 Hz劃分為24個頻段，每個頻段范圍如圖 3(a)所示，求得主元個數L=6，求得各頻段歸一化能量比例系數如圖 3(b)所示。
圖3
24個頻段的范圍和歸一化能量系數
(a)各頻段對應的頻率范圍;(b)各頻段的歸一化能量系數，橫軸對應各頻段的中心頻率
Figure3. Bandwidths and normalized energy coefficients of 24 sub-bands(a) bandwidths of 24 sub-bands; (b) normalized energy coefficients of 24 sub-bands,the horizontal axis is the central frequency of each band
2.3 非均勻濾波器的分配
對圖 3(b)中的歸一化能量系數曲線進行插值得到頻譜能量曲線，對能量曲線的積分進行均勻劃分，以分界點作為濾波器的中心點，每個濾波器的帶寬為其前一個濾波器中點至其后一個濾波器中點，第一個濾波器的下限頻率為0 Hz，最后一個濾波器的上限頻率為4 000 Hz。具體方法如下：
令n維向量y為能量比重系數，
| ${y_k} = {\alpha _k} = \sum\limits_{j = 1}^L {{\lambda _j}{P_{kj}}} ,k = 1,2, \cdots ,n$ |
n維向量x為各頻段的中心頻段對應的傅里葉變換點數，n=24。以256點傅里葉變換為例，對[x,y]采用三次樣條函數插值為128維向量。
若濾波器的個數為m，將該插值曲線積分平均分為(m+1)等份，分界點即為每個濾波器中心對應的快速傅里葉變換點，第k個濾波器中心對應的點數為Ck，滿足
| $\sum\limits_{i = {c_{k - 1}}}^{{c_k}} {{y_i}} = \frac{1}{m}\sum\limits_{i = 1}^{128} {{y_i}} $ |
以24個濾波器為例，得到的Mel刻度上非均勻濾波器組如圖 4所示。
圖4
非均勻濾波器組
Figure4.
Non-uniform filter-bank
3 特征提取
得到Mel刻度上非均勻分布的濾波器組后，用其取代傳統的均勻分布濾波器組進行改進的MFCC提取，流程如圖 5所示。
圖5
改進的MFCC特征提取流程
Figure5.
Improved process of MFCC extraction
本文采用12階MFCC參數和對數能量，并求其一階差分和二階差分，得到共39維參數作為一幀咳嗽信號的特征參數。
4 實驗方法與結果
采用4個隱含狀態的HMM作為分類器識別咳嗽和非咳嗽信號，每個狀態包含3個高斯混合模型，采用Kmeans方法初始化高斯混合模型參數，訓練過程采用Baum-Welch算法，識別采用Viterbi算法[14]。
HMM進行咳嗽識別的步驟如下：
(1)為咳嗽信號和非咳嗽信號分別定義一個HMM，記為V={v1,v2}；
(2)選取咳嗽和非咳嗽信號作為訓練樣本，訓練樣本分別取60個咳嗽信號和60個非咳嗽信號，測試樣本為92個咳嗽信號和100個非咳嗽信號；
(3)分別訓練vi(i=1,2)，獲得最佳模型λi；
(4)采用訓練好的模型λi(i=1,2)進行識別。對未知的輸入信號O計算Pr(O｜λi)(i=1,2)，若Pr(O｜λj)=max(Pr(O,λi),i=1,2)，則將輸入信號O識別為vj類。
改進方法和傳統MFCC的識別結果對比如表 1所示，其中，靈敏度=正確識別的咳嗽個數/咳嗽總數×100%；特異度=正確識別的非咳嗽個數/非咳嗽總數×100%；漏識別率＝誤識別為非咳嗽的咳嗽個數/咳嗽總數×100%；誤識別率＝誤識別為咳嗽的非咳嗽個數/非咳嗽總數×100%。
表1
咳嗽識別結果
Table1.
Results of cough recognition
由表 1可見，取不同濾波器個數時，本文方法得到的靈敏度、特異度、漏識別率和誤識別率相對傳統的MFCC在總體上都有改善，其中靈敏度改善最明顯，這是由于咳嗽頻譜的高能量頻段包含了更多的咳嗽特征信息，非均勻濾波器組在這些頻段放置了較多的濾波器，使得特征提取更加有效。同時，低能量頻段由于包含特征信息較少且容易受干擾信號的影響，放置較少的濾波器可以避免提取到較多的非特征信息，從而提高識別的準確率和抗干擾性。該方法中的頻率能量分布是通過對大量咳嗽樣本進行統計分析得到的，因此咳嗽樣本的選取對濾波器的分配有直接的影響，統計樣本的性質與待識別的咳嗽相近則改進后的識別效果較好，如何更加合理地選擇咳嗽樣本的類型、數量以提高方法的通用性還需要進一步探討。
5 小結
本文提出了根據咳嗽頻譜能量的分布特點設計Mel刻度上非均勻濾波器組進行MFCC特征提取的方法。該方法在提取特征時將濾波器主要放置在咳嗽頻譜能量高的頻段，突出了咳嗽具有高頻譜能量這一顯著特點，使得特征提取更加有效，從而改善了咳嗽識別的效果。該方法是基于大量咳嗽的統計特性提出的，對錄音環境和訓練樣本選取的要求較高，下一步工作是進一步減小系統對訓練樣本的依賴性和提高系統的抗干擾能力。
引言
咳嗽是呼吸系統疾病常見的癥狀，咳嗽信號的記錄和自動識別在臨床上具有重要的診斷指導意義。咳嗽記錄主要通過硬件實現，目前已有多種類型的便攜式咳嗽監測儀，能有效實現對咳嗽的長時間記錄[1-3]。咳嗽記錄結果的分析目前仍主要靠人工分析，而咳嗽信號的自動識別是咳嗽監測儀推廣應用的關鍵，有效的特征提取方法則是實現高精度識別的基礎。
過去幾年里，研究者主要參照語音識別的方法，以傳統梅爾(Mel)頻率倒譜系數(Mel frequency cepstrum coefficients,MFCC)為特征參數，采用動態時間規整、隱馬爾可夫模型(hidden Markov model,HMM)、高斯混合模型或神經網絡為分類器來實現咳嗽識別[4-8]。咳嗽信號與語音信號相比，在時域和頻域特性都有其自身的特點，在特征提取或識別的時候充分考慮咳嗽信號的特性將會使算法更加有效。
傳統MFCC的計算首先用快速傅立葉變換將時域信號轉化成頻域，之后對其對數能量譜依照Mel刻度上均勻分布的三角濾波器組進行卷積，最后對各個濾波器的輸出構成的向量進行離散余弦變換，取前N個系數[9]。雖然Mel刻度反映了人耳對不同信號的感知特點，但事實上，咳嗽頻譜能量在Mel刻度上并非均勻分布，即并不隨頻率的增加而遞減，而是在某些頻段出現峰值[10]，例如慢性支氣管炎患者咳嗽的頻率集中在500、700和1 200 Hz附近[11]，又如圖 1所示是一個典型干性咳嗽的時域波形和語譜圖,其能量峰值出現在[300,1 000]Hz和[2 000，3 000] Hz范圍。這種情況下采用Mel刻度上均勻分布的濾波器組，會出現頻譜能量低的頻段反而放置了較多濾波器的問題，不能最大程度地獲得咳嗽特征。基于這一點，本文將整個頻域范圍劃分為若干個頻段，采用主元分析法分析咳嗽信號的頻譜能量分布，按各頻段的能量強度系數分配濾波器，設計Mel刻度上的非均勻濾波器組，在高能量頻段放置較多濾波器，在低能量頻段放置較少濾波器，改進傳統MFCC參數的提取過程。
圖1
咳嗽的時域波形和語譜圖
Figure1.
Time domain waveform and spectrogram of a typical cough
1 數據來源
實驗樣本來自廣州醫學院第一附屬醫院病房，使用DELL原裝機的聲卡信號采集設備進行單聲道采樣,16位數字量化。盡管文獻[11]和[12]都報道觀察到了5 000～6 000 Hz的咳嗽分量，但3 500 Hz以上頻段的咳嗽分量很少[10]，忽略這部分高頻成分對咳嗽特征提取影響不大，因此為了減少計算量，本文采樣頻率取8 000 Hz，對應信號的處理范圍為0～4 000 Hz。選取不同性別、年齡和呼吸道感染情況的患者進行連續、長時間的日常錄音，人工分割成單個咳嗽信號作為咳嗽樣本，并隨機截取時長0.2～3 s的非咳嗽信號作為非咳嗽樣本。
2 非均勻濾波器組設計
2.1 總體流程
非均勻濾波器組設計的基本思想是根據咳嗽的頻譜能量分布來分配濾波器個數，在高能量頻段放置較多濾波器，在低能量頻段放置較少濾波器。首先將整個頻域范圍劃分為若干個頻段，采用主元分析法對大量的咳嗽樣本進行分析，得到咳嗽信號的頻譜能量分布系數，然后插值得到頻譜能量曲線，對能量曲線的積分進行均勻劃分，分界點即為濾波器的中心點。非均勻濾波器組的設計總體流程如圖 2所示。
圖2
非均勻濾波器組的總體設計流程
Figure2.
Overall design process of non-uniform filter-bank
2.2 主元分析法求頻譜能量分布
按幀長240(30 ms)、幀移80(10 ms)對咳嗽信號進行分幀，將每幀[0,4 000] Hz的頻率范圍劃分為n個頻段，每個頻段與Mel刻度上n個濾波器的均勻濾波器組的通帶范圍對應，則每幀信號由這n個頻段的頻譜按一定方式組合而成，記第i幀語音信號為fi，則fi可表示為n個頻段能量的線性組合:
| ${f_i} = \sum\limits_{k = 1}^n {{\omega _{ik}}{x_{ik}},} $ |
其中xik為第i幀信號第k個頻段的頻譜能量和；ωik為第i幀信號第k個頻段的能量系數，它反映了第k個頻段的能量在該幀信號頻譜中所占的比例。對多幀的咳嗽信號可得到能量矩陣A：
| $A = \left( {\begin{array}{*{20}{c}} {{x_{11}}}&{{x_{12}}}& \cdots &{{x_{1n}}}\\ {{x_{21}}}&{{x_{22}}}& \cdots &{{x_{2n}}}\\ \cdots & \cdots & \cdots & \cdots \\ {{x_{m1}}}&{{x_{m2}}}& \cdots &{{x_{mn}}} \end{array}} \right)$ |
其中m是總的幀數，n為頻段數。主元分析法能夠找出數據中最重要的元素[13]，對A進行主元分析，找到主元對應的頻段系數即為統計意義上的頻段能量系數。具體步驟如下。
步驟1：計算能量矩陣A。對m幀咳嗽信號分別求快速傅里葉變換，將頻段在Mel刻度上平均劃分為n個頻段，計算每幀每個頻段的頻譜能量，得到m*n的能量矩陣A。
步驟2：計算協方差矩陣ATA，得到n*n的方陣。
步驟3：計算協方差矩陣ATA的特征值λ和特征向量P，滿足
| ${P^T}{A^T}AP = \Lambda ,$ |
| $A = \left( {\begin{array}{*{20}{c}} {\sqrt {{\lambda _1}} }&0& \cdots &0\\ 0&{\sqrt {{\lambda _2}} }& \cdots &0\\ \cdots & \cdots & \cdots & \cdots \\ 0&0& \cdots &{\sqrt {{\lambda _n}} } \end{array}} \right),$ |
其中λ按從大到小排序，按85%的累計貢獻率，得到主元個數為L。
步驟4：求各頻段的能量比重，第k個頻段的能量比重為
| ${\alpha _k} = \sum\limits_{j = 1}^L {{\lambda _j}{P_{kj}}} ,k = 1,2, \cdots ,n$ |
本文選取了92個典型的咳嗽信號，共4 015幀，按Mel刻度上均勻濾波器組的帶寬范圍將0～4 000 Hz劃分為24個頻段，每個頻段范圍如圖 3(a)所示，求得主元個數L=6，求得各頻段歸一化能量比例系數如圖 3(b)所示。
圖3
24個頻段的范圍和歸一化能量系數
(a)各頻段對應的頻率范圍;(b)各頻段的歸一化能量系數，橫軸對應各頻段的中心頻率
Figure3. Bandwidths and normalized energy coefficients of 24 sub-bands(a) bandwidths of 24 sub-bands; (b) normalized energy coefficients of 24 sub-bands,the horizontal axis is the central frequency of each band
2.3 非均勻濾波器的分配
對圖 3(b)中的歸一化能量系數曲線進行插值得到頻譜能量曲線，對能量曲線的積分進行均勻劃分，以分界點作為濾波器的中心點，每個濾波器的帶寬為其前一個濾波器中點至其后一個濾波器中點，第一個濾波器的下限頻率為0 Hz，最后一個濾波器的上限頻率為4 000 Hz。具體方法如下：
令n維向量y為能量比重系數，
| ${y_k} = {\alpha _k} = \sum\limits_{j = 1}^L {{\lambda _j}{P_{kj}}} ,k = 1,2, \cdots ,n$ |
n維向量x為各頻段的中心頻段對應的傅里葉變換點數，n=24。以256點傅里葉變換為例，對[x,y]采用三次樣條函數插值為128維向量。
若濾波器的個數為m，將該插值曲線積分平均分為(m+1)等份，分界點即為每個濾波器中心對應的快速傅里葉變換點，第k個濾波器中心對應的點數為Ck，滿足
| $\sum\limits_{i = {c_{k - 1}}}^{{c_k}} {{y_i}} = \frac{1}{m}\sum\limits_{i = 1}^{128} {{y_i}} $ |
以24個濾波器為例，得到的Mel刻度上非均勻濾波器組如圖 4所示。
圖4
非均勻濾波器組
Figure4.
Non-uniform filter-bank
3 特征提取
得到Mel刻度上非均勻分布的濾波器組后，用其取代傳統的均勻分布濾波器組進行改進的MFCC提取，流程如圖 5所示。
圖5
改進的MFCC特征提取流程
Figure5.
Improved process of MFCC extraction
本文采用12階MFCC參數和對數能量，并求其一階差分和二階差分，得到共39維參數作為一幀咳嗽信號的特征參數。
4 實驗方法與結果
采用4個隱含狀態的HMM作為分類器識別咳嗽和非咳嗽信號，每個狀態包含3個高斯混合模型，采用Kmeans方法初始化高斯混合模型參數，訓練過程采用Baum-Welch算法，識別采用Viterbi算法[14]。
HMM進行咳嗽識別的步驟如下：
(1)為咳嗽信號和非咳嗽信號分別定義一個HMM，記為V={v1,v2}；
(2)選取咳嗽和非咳嗽信號作為訓練樣本，訓練樣本分別取60個咳嗽信號和60個非咳嗽信號，測試樣本為92個咳嗽信號和100個非咳嗽信號；
(3)分別訓練vi(i=1,2)，獲得最佳模型λi；
(4)采用訓練好的模型λi(i=1,2)進行識別。對未知的輸入信號O計算Pr(O｜λi)(i=1,2)，若Pr(O｜λj)=max(Pr(O,λi),i=1,2)，則將輸入信號O識別為vj類。
改進方法和傳統MFCC的識別結果對比如表 1所示，其中，靈敏度=正確識別的咳嗽個數/咳嗽總數×100%；特異度=正確識別的非咳嗽個數/非咳嗽總數×100%；漏識別率＝誤識別為非咳嗽的咳嗽個數/咳嗽總數×100%；誤識別率＝誤識別為咳嗽的非咳嗽個數/非咳嗽總數×100%。
表1
咳嗽識別結果
Table1.
Results of cough recognition
由表 1可見，取不同濾波器個數時，本文方法得到的靈敏度、特異度、漏識別率和誤識別率相對傳統的MFCC在總體上都有改善，其中靈敏度改善最明顯，這是由于咳嗽頻譜的高能量頻段包含了更多的咳嗽特征信息，非均勻濾波器組在這些頻段放置了較多的濾波器，使得特征提取更加有效。同時，低能量頻段由于包含特征信息較少且容易受干擾信號的影響，放置較少的濾波器可以避免提取到較多的非特征信息，從而提高識別的準確率和抗干擾性。該方法中的頻率能量分布是通過對大量咳嗽樣本進行統計分析得到的，因此咳嗽樣本的選取對濾波器的分配有直接的影響，統計樣本的性質與待識別的咳嗽相近則改進后的識別效果較好，如何更加合理地選擇咳嗽樣本的類型、數量以提高方法的通用性還需要進一步探討。
5 小結
本文提出了根據咳嗽頻譜能量的分布特點設計Mel刻度上非均勻濾波器組進行MFCC特征提取的方法。該方法在提取特征時將濾波器主要放置在咳嗽頻譜能量高的頻段，突出了咳嗽具有高頻譜能量這一顯著特點，使得特征提取更加有效，從而改善了咳嗽識別的效果。該方法是基于大量咳嗽的統計特性提出的，對錄音環境和訓練樣本選取的要求較高，下一步工作是進一步減小系統對訓練樣本的依賴性和提高系統的抗干擾能力。
