咳嗽的自動分類在臨床上具有重要的輔助診斷作用。傳統的Mel頻率倒譜系數(MFCC)采用Mel均勻濾波器組, 高頻段的濾波器分布較稀疏, 未能最大程度反映兩類咳嗽的特征差別。針對這個問題, 本文在分析干性咳嗽和濕性咳嗽頻譜能量分布特點的基礎上, 提出了一種改進的反向MFCC提取方法, 采用反向Mel刻度上的均勻濾波器組, 并放置在兩類咳嗽都具有高頻譜能量的頻段, 使得特征提取集中在兩類咳嗽特征信息豐富且差別顯著的頻段進行。基于隱馬爾可夫模型的咳嗽干濕性自動分類實驗結果表明, 該方法獲得了優于傳統MFCC的分類性能, 總體分類準確率從89.76%提高到了93.66%。
引用本文： 朱春媚, 劉保軍, 黎萍, 莫鴻強, 鄭則廣. 基于改進反向Mel頻率倒譜系數的咳嗽干濕性自動分類. 生物醫學工程學雜志, 2016, 33(2): 239-243, 254. doi: 10.7507/1001-5515.20160042 復制
引言
咳嗽是呼吸系統疾病常見的癥狀[1]，咳嗽的監測和自動識別在臨床上具有重要的診斷指導意義。目前，咳嗽聲音的分析和識別主要參考語音識別系統，采用神經網絡[2-3]、隱馬爾可夫模型[4-5]、隱馬爾可夫模型與神經網絡的混合模型[6]或分類樹[7]等方法來實現。通過借鑒語音識別方面的技術，咳嗽識別的準確率已達到較高的水平，目前已有多種類型的便攜式咳嗽監測儀，如Hull自動化咳嗽計數器[8]、Leicester咳嗽監測系統[9]、Lifeshirt系統[10-11]和VitaloJAK系統[12]等，均在一定程度上實現了咳嗽的自動或半自動識別。
咳嗽按性質可分為干性咳嗽和濕性咳嗽兩種，干性咳嗽無痰或痰量少，常見于慢性喉炎、氣管炎、氣管受壓或管內異物；濕性咳嗽有痰，常見于慢性支氣管炎、肺炎、支氣管擴張和肺膿腫等[13-14]。咳嗽的特征取決于潛在疾病的種類，因此咳嗽的類型在一定程度上能夠反映呼吸道疾病的類型及其嚴重程度。咳嗽類型的自動分類有助于醫生快速診斷病因和提高療效，是計算機輔助診斷應用于呼吸系統疾病的一個重要組成內容。
相對于咳嗽的自動識別，國內外關于咳嗽分類的研究報道不多。文獻[15]采用Mel頻率倒譜系數(Mel frequency cepstrum coefficients, MFCC)和動態時間規整的方法來對咳嗽信號進行自動干濕性分類，得到干性咳嗽和濕性咳嗽的分類準確率分別為94.55%和73.91%；文獻[16]采用MFCC、共振峰頻率和過零率等特征，通過Logistic回歸模型識別兒科患者的濕性咳嗽，準確率為84%；文獻[17]提出了直接利用咳嗽的時域特征和頻域特征進行干濕性分類的方法，并用16個典型樣本測試了其有效性。可見，高準確率并通過大量測試樣本驗證的咳嗽分類方法還有待進一步的研究。
在以上咳嗽識別和咳嗽分類的研究中，一般以MFCC作為特征量。傳統MFCC的提取采用Mel刻度上均勻分布的濾波器組，由于實際頻率與Mel頻率的非線性對應關系，濾波器組在低頻區域分布較密集，在高頻區域分布較稀疏，即特征提取集中在低頻區域進行。但事實上，根據咳嗽的發聲機制[18-19]，當發生濕性咳嗽時，由于呼吸道中有痰導致氣流震蕩加劇，使得咳嗽聲的中頻和高頻成分增多，因此干性咳嗽和濕性咳嗽的頻譜特征在中頻和高頻區域有較顯著的差別，如圖 1所示。傳統MFCC由于濾波器過多集中在低頻段，中頻和高頻段提取的信息不足，未能最大程度反映兩類咳嗽的差別，進而影響咳嗽分類效果。針對這個問題，本文提出一種改進的反向MFCC特征提取方法，將特征提取集中在兩類咳嗽頻譜能量高且差別顯著的頻段進行，并采用隱馬爾可夫模型實現咳嗽的干濕性自動分類。
圖1
咳嗽的時域波形和語譜圖
Figure1.
Time domain waveform and spectrogram of cough
1 改進特征提取方法
1.1 傳統MFCC提取
傳統MFCC的提取流程如圖 2所示。
圖2
傳統MFCC提取流程
Figure2.
Process of traditional MFCC extraction
具體包含以下步驟：
(1)?對輸入的咳嗽信號進行分幀、加窗，然后作離散傅里葉變換，獲得頻譜分布信息。
設第i幀咳嗽信號的離散傅里葉變換為：
| $ {X_i}\left(k \right)=\sum\limits_{n=1}^{N-1} {x\left(n \right){{\rm{e}}^{-j2\pi k/N}}, 0 \le k \le N} $ |
式中x(n)為輸入的咳嗽信號，N表示傅立葉變換的點數。
(2)?求頻譜幅度的平方，得到能量譜：
| $ {E_i}\left(k \right)={\left| {{X_i}\left(k \right)} \right|^2} $ |
(3)?將能量譜通過均勻Mel濾波器組。
定義一個在Mel刻度上均勻分布的三角濾波器組，包含M個濾波器，中心頻率為f(m), m=1, 2, …, M。當M=20時，均勻濾波器組的結構如圖 3所示。
圖3
Mel濾波器組結構
Figure3.
Structure of Mel filter bank
濾波器的頻率響應定義為:
| $ {H_m}\left(k \right)=\left\{ {\begin{array}{*{20}{c}} {\frac{{2\left({k-f\left[ {m-1} \right]} \right)}}{{\left({f\left[ {m + 1} \right]-f\left[ {m-1} \right]} \right)\left({f\left[ m \right]-f\left[ {m + 1} \right]} \right)}}, f\left[ {m-1} \right] < k \le f\left[ m \right]}\\ {\frac{{2\left({f\left[ {m-1} \right]-k} \right)}}{{\left({f\left[ {m + 1} \right]-f\left[ {m-1} \right]} \right)\left({f\left[ {m + 1} \right]-f\left[ m \right]} \right)}}, f\left[ m \right] < k < f\left[ {m + 1} \right]}\\ {0, k \le f\left[ {m-1} \right]或k \ge f\left[ {m + 1} \right]} \end{array}} \right. $ |
(4)?計算每個濾波器組輸出的對數能量為:
| $ S\left(m \right)=\ln \left({\sum\limits_{k=0}^{N-1} {{E_i}\left(k \right){H_m}\left(k \right)} } \right), 0 \le m < M $ |
(5)?經離散余弦變換得到靜態MFCC系數：
| $ C\left(n \right)=\sum\limits_{m=0}^{N-1} {S\left(m \right)\cos \left({\frac{{\pi l\left({m-0.5} \right)}}{M}} \right)}, l=1, 2, \cdots, L $ |
式中L為MFCC系數的階數，通常取L=12。
(6)?差分參數的計算可以采用下面的公式[20]：
| $ {d_t}=\left\{ {\begin{array}{*{20}{c}} {{C_{t + 1}}-{C_t}, t < K}\\ {\frac{{\sum\limits_{k=1}^K {k\left({{C_{t + k}}-{C_{t-k}}} \right)} }}{{\sqrt {2\sum\limits_{k-1}^K {{k^2}} } }}, 其它}\\ {{C_t}-{C_{t-1}}, t \ge Q-K} \end{array}} \right. $ |
式中，dt表示第t個一階差分；Ct表示第t個倒譜系數；Q表示倒譜系數的階數；K表示一階導數的時間差，可取1或2。
1.2 反向MFCC提取
反向MFCC的基本思想是在反向Mel頻率刻度上放置均勻的濾波器組，從而使得高頻區域的濾波器分布比低頻區域密集。反向Mel頻率刻度fI-mel與實際頻率之間的轉換如式(7)所示[21]：
| $ {f_{I-mel}}=2195.29-2595{\log _{10}}\left({1 + \frac{{4031.25-g}}{{700}}} \right) $ |
其中，f為實際頻率，單位為Hz。
記反向Mel濾波器響應為，與均勻Mel濾波器相應Hm(k)對應關系如下：
| $ \overline {{H_m}} \left(k \right)={H_{M-m + 1}}\left({\frac{N}{2}-k + 1} \right) $ |
用取代傳統MFCC提取流程中的Hm(k)，即可得到反向MFCC。
反向Mel濾波器組的結構如圖 4所示，可見，與圖 3中的Mel濾波器組相反，反向Mel濾波器組在高頻區域放置較多的濾波器，分布較密集；在低頻區域分布較稀疏。
圖4
反向Mel濾波器組結構
Figure4.
Structure of reverse Mel filter bank
1.3 改進反向MFCC提取
當采樣頻率為8 000 Hz時，在咳嗽信號的整個[0, 4 000]Hz頻域范圍內，0 Hz和4 000 Hz附近頻段的頻譜能量很小[18-19]，低能量頻段包含特征信息少且容易受噪聲影響，咳嗽的特征主要集中在頻譜能量較高的頻段。因此，在進行咳嗽的干濕性分類時，為了最大程度地反映干性咳嗽和濕性咳嗽的區別，有效的特征提取應集中在兩類咳嗽頻譜能量高且差別較顯著的頻段進行。
雖然反向MFCC濾波器組有效地獲取了高頻區域的頻譜信息，但濾波器過多地集中在了咳嗽頻譜能量較小的4 000 Hz附近的頻段上，仍不能滿足咳嗽分類的要求。本文根據咳嗽頻譜能量分布對反向MFCC濾波器組進行改進，在咳嗽頻譜能量較大的[fl, fh]Hz的頻率范圍放置反向Mel濾波器組，其中fl和fh分別對應于高頻譜能量頻段的下限頻率和上限頻率，且fl > 0 Hz，fh < 4 000 Hz。
相對于反向Mel濾波器組的響應，計算改進反向Mel濾波器組響應的區別主要在于中心頻率f′[m]的計算。改進反向Mel濾波器組響應的計算步驟如下：
(1)?根據實際頻率與Mel頻率的對應關系fmel=2 595log10(1+f/700)，得到[fl, fh]Hz對應的Mel頻率范圍為[fmel-l, fmel-h]；
(2)?將Mel頻率范圍[fmel-l, fmel-h]平均劃分為(M+1)份，M為濾波器個數，得到Mel刻度下第m個濾波器的中心頻率fMel[m]=m(fmel-h-fmel-l)/M；
(3)?將Mel刻度的中心頻率fMed[m]轉為以Hz為單位，得到改進反向Mel濾波器的中心頻率f′[m]。
(4)?將f′[m]代入公式(3)和(8)計算得到的，即為改進反向Mel濾波器組的響應。
(5)?用改進反向Mel濾波器組的響應取代傳統MFCC提取流程中的Hm(k)，得到改進反向MFCC。
文中取M=20, fl=55 Hz，fh=3 040 Hz，得到改進的反向Mel濾波器組結構如圖 5所示。可見，改進后的反向濾波器組集中分布在[55, 3 040] Hz頻率范圍的中、高頻區域，即特征提取集中在兩類咳嗽頻譜能量高且差別較顯著的頻段進行。
圖5
改進的反向Mel濾波器組結構
Figure5.
Structure of improved reverse Mel filter bank
2 實驗方法和結果
2.1 數據來源
實驗樣本來自廣州醫學院第一附屬醫院病房，使用DELL原裝機的聲卡信號采集設備進行單聲道采樣，采樣頻率為8 000 Hz，每個錄音保存為16位數字量化的WAV文件。選取不同性別、年齡和呼吸感染情況的患者13例(其中以干性咳嗽為主的患者7例，男4例，女3例，平均年齡53歲；以濕性咳嗽為主的患者例數6例，男3例，女3例，平均年齡59歲), 在患者知情同意的情況下，進行連續、長時間的日常錄音。將錄音人工分割成單個咳嗽信號并在醫生指導下分類，作為訓練樣本和測試樣本。
2.2 隱馬爾可夫模型咳嗽分類
隱馬爾可夫模型由于其優異的動態建模能力在語言識別中獲得了廣泛的應用[22]，本文采用隱馬爾可夫模型來實現干性咳嗽和濕性咳嗽嗽的分類。分別取60個干性咳嗽(其中男32個，女28個；無噪聲環境下46個，有噪聲環境下14個)和60個濕性咳嗽(其中男36個，女24個；無噪聲環境下41個，噪聲環境下19個)作為隱馬爾可夫模型訓練樣本，再分別取120個干性咳嗽(其中男63個，女57個；無噪聲環境下96個，有噪聲環境下24個)和85個濕性咳嗽(其中男49個，女36個；無噪聲環境下61個，有噪聲環境下24個)作為測試樣本，測試樣本不包含訓練樣本。
采用隱馬爾可夫模型進行咳嗽干濕性分類的步驟如下：
(1)?為干性咳嗽信號和濕性咳嗽信號分別定義一個隱馬爾可夫模型，記為H={h1, h2}。采用含2個隱含狀態的自左向右隱馬爾可夫模型來實現，每個狀態為包含3個高斯元的混合模型，采用C均值的方法初始化高斯混合模型參數。
(2)?取60個干性咳嗽和60個濕性咳嗽作為訓練樣本，分別訓練hi(i=1, 2)，獲得最佳模型μi。采用Baum-Welch算法[22]進行訓練，迭代次數設置為50次。
(3)?采用訓練好的模型μi(i=1, 2)進行咳嗽分類。采用Viterbi算法[22]對未知的輸入信號O，計算Pr(O｜μi)(i=1, 2)，若Pr(O｜μi)=max(Pr(O, μi), i=1, 2))，則將輸入信號O分類為hj類。
2.3 實驗結果
分別采用包含一階和二級差分系數的39維傳統MFCC、反向MFCC和改進反向MFCC作為特征量，采用隱馬爾可夫模型進行咳嗽的干濕性分類，結果如表 1所示，其中干性咳嗽的分類準確率=正確分類的干性咳嗽個數/干性咳嗽總數；濕性咳嗽的分類準確率=正確分類的濕性咳嗽個數/濕性咳嗽總數；總體分類準確率=正確分類的咳嗽個數/咳嗽總數。
由表 1可見，三種特征提取方法中，反向MFCC對應的干性咳嗽、濕性咳嗽和總體分類準確率都最低，而改進反向MFCC相應的三種分類準確率都最高。另外采用三種不同的特征時，濕性咳嗽分類準確率皆低于干性咳嗽，反映了濕性咳嗽的頻譜構成更復雜和多樣化，這與文獻[19]報道的咳嗽特點一致。
表1
三種特征提取方法的分類結果比較
Table1.
Comparison of cough classification results respectively by the three feature extraction methods
3 討論
在咳嗽識別和咳嗽分類研究中，以往的研究大多參考語音識別系統，采用傳統的MFCC作為特征量，傳統MFCC在Mel刻度放置均勻濾波器，特征提取集中在信號的低頻和中頻區域進行。然而通過分析咳嗽的機制和聲學特征發現，在咳嗽信號的整個0～4 000 Hz頻域范圍內，0 Hz和4 000 Hz附近頻段的頻譜能量很小，咳嗽的頻譜能量主要分布在50～3 050 Hz范圍內，且濕性咳嗽由于呼吸道中有痰導致氣流震蕩加劇，其頻譜的中頻和高頻成分比干性咳嗽更加豐富。咳嗽分類的目標是區分兩類咳嗽，因此有效的特征提取應集中在兩類咳嗽頻譜能量高且差別較顯著的頻段進行。本文提出的改進反向MFCC，將濾波器組集中放置在兩類咳嗽頻譜能量高且差別較顯著的頻段，實現了有效的特征提取，從而獲得了優于傳統MFCC的分類結果，使得總體分類準確率從89.76%提高到了93.66%；該方法主要提取了咳嗽中頻和高頻段的特征，使得該頻段特征明顯的濕性咳嗽更容易被正確識別，因而采用該方法時濕性咳嗽分類準確率的提升比干性咳嗽顯著。另一方面，與傳統MFCC一樣，改進反向MFCC仍是反映信號頻譜能量的特征表示方法，因而容易受到噪聲的影響，當干性咳嗽疊加了噪聲，尤其是中頻和高頻的噪聲時，由于中頻和高頻成分增多，容易被錯分類為濕性咳嗽。如何進一步提高特征提取方法的魯棒性是我們將來的研究重點。
引言
咳嗽是呼吸系統疾病常見的癥狀[1]，咳嗽的監測和自動識別在臨床上具有重要的診斷指導意義。目前，咳嗽聲音的分析和識別主要參考語音識別系統，采用神經網絡[2-3]、隱馬爾可夫模型[4-5]、隱馬爾可夫模型與神經網絡的混合模型[6]或分類樹[7]等方法來實現。通過借鑒語音識別方面的技術，咳嗽識別的準確率已達到較高的水平，目前已有多種類型的便攜式咳嗽監測儀，如Hull自動化咳嗽計數器[8]、Leicester咳嗽監測系統[9]、Lifeshirt系統[10-11]和VitaloJAK系統[12]等，均在一定程度上實現了咳嗽的自動或半自動識別。
咳嗽按性質可分為干性咳嗽和濕性咳嗽兩種，干性咳嗽無痰或痰量少，常見于慢性喉炎、氣管炎、氣管受壓或管內異物；濕性咳嗽有痰，常見于慢性支氣管炎、肺炎、支氣管擴張和肺膿腫等[13-14]。咳嗽的特征取決于潛在疾病的種類，因此咳嗽的類型在一定程度上能夠反映呼吸道疾病的類型及其嚴重程度。咳嗽類型的自動分類有助于醫生快速診斷病因和提高療效，是計算機輔助診斷應用于呼吸系統疾病的一個重要組成內容。
相對于咳嗽的自動識別，國內外關于咳嗽分類的研究報道不多。文獻[15]采用Mel頻率倒譜系數(Mel frequency cepstrum coefficients, MFCC)和動態時間規整的方法來對咳嗽信號進行自動干濕性分類，得到干性咳嗽和濕性咳嗽的分類準確率分別為94.55%和73.91%；文獻[16]采用MFCC、共振峰頻率和過零率等特征，通過Logistic回歸模型識別兒科患者的濕性咳嗽，準確率為84%；文獻[17]提出了直接利用咳嗽的時域特征和頻域特征進行干濕性分類的方法，并用16個典型樣本測試了其有效性。可見，高準確率并通過大量測試樣本驗證的咳嗽分類方法還有待進一步的研究。
在以上咳嗽識別和咳嗽分類的研究中，一般以MFCC作為特征量。傳統MFCC的提取采用Mel刻度上均勻分布的濾波器組，由于實際頻率與Mel頻率的非線性對應關系，濾波器組在低頻區域分布較密集，在高頻區域分布較稀疏，即特征提取集中在低頻區域進行。但事實上，根據咳嗽的發聲機制[18-19]，當發生濕性咳嗽時，由于呼吸道中有痰導致氣流震蕩加劇，使得咳嗽聲的中頻和高頻成分增多，因此干性咳嗽和濕性咳嗽的頻譜特征在中頻和高頻區域有較顯著的差別，如圖 1所示。傳統MFCC由于濾波器過多集中在低頻段，中頻和高頻段提取的信息不足，未能最大程度反映兩類咳嗽的差別，進而影響咳嗽分類效果。針對這個問題，本文提出一種改進的反向MFCC特征提取方法，將特征提取集中在兩類咳嗽頻譜能量高且差別顯著的頻段進行，并采用隱馬爾可夫模型實現咳嗽的干濕性自動分類。
圖1
咳嗽的時域波形和語譜圖
Figure1.
Time domain waveform and spectrogram of cough
1 改進特征提取方法
1.1 傳統MFCC提取
傳統MFCC的提取流程如圖 2所示。
圖2
傳統MFCC提取流程
Figure2.
Process of traditional MFCC extraction
具體包含以下步驟：
(1)?對輸入的咳嗽信號進行分幀、加窗，然后作離散傅里葉變換，獲得頻譜分布信息。
設第i幀咳嗽信號的離散傅里葉變換為：
| $ {X_i}\left(k \right)=\sum\limits_{n=1}^{N-1} {x\left(n \right){{\rm{e}}^{-j2\pi k/N}}, 0 \le k \le N} $ |
式中x(n)為輸入的咳嗽信號，N表示傅立葉變換的點數。
(2)?求頻譜幅度的平方，得到能量譜：
| $ {E_i}\left(k \right)={\left| {{X_i}\left(k \right)} \right|^2} $ |
(3)?將能量譜通過均勻Mel濾波器組。
定義一個在Mel刻度上均勻分布的三角濾波器組，包含M個濾波器，中心頻率為f(m), m=1, 2, …, M。當M=20時，均勻濾波器組的結構如圖 3所示。
圖3
Mel濾波器組結構
Figure3.
Structure of Mel filter bank
濾波器的頻率響應定義為:
| $ {H_m}\left(k \right)=\left\{ {\begin{array}{*{20}{c}} {\frac{{2\left({k-f\left[ {m-1} \right]} \right)}}{{\left({f\left[ {m + 1} \right]-f\left[ {m-1} \right]} \right)\left({f\left[ m \right]-f\left[ {m + 1} \right]} \right)}}, f\left[ {m-1} \right] < k \le f\left[ m \right]}\\ {\frac{{2\left({f\left[ {m-1} \right]-k} \right)}}{{\left({f\left[ {m + 1} \right]-f\left[ {m-1} \right]} \right)\left({f\left[ {m + 1} \right]-f\left[ m \right]} \right)}}, f\left[ m \right] < k < f\left[ {m + 1} \right]}\\ {0, k \le f\left[ {m-1} \right]或k \ge f\left[ {m + 1} \right]} \end{array}} \right. $ |
(4)?計算每個濾波器組輸出的對數能量為:
| $ S\left(m \right)=\ln \left({\sum\limits_{k=0}^{N-1} {{E_i}\left(k \right){H_m}\left(k \right)} } \right), 0 \le m < M $ |
(5)?經離散余弦變換得到靜態MFCC系數：
| $ C\left(n \right)=\sum\limits_{m=0}^{N-1} {S\left(m \right)\cos \left({\frac{{\pi l\left({m-0.5} \right)}}{M}} \right)}, l=1, 2, \cdots, L $ |
式中L為MFCC系數的階數，通常取L=12。
(6)?差分參數的計算可以采用下面的公式[20]：
| $ {d_t}=\left\{ {\begin{array}{*{20}{c}} {{C_{t + 1}}-{C_t}, t < K}\\ {\frac{{\sum\limits_{k=1}^K {k\left({{C_{t + k}}-{C_{t-k}}} \right)} }}{{\sqrt {2\sum\limits_{k-1}^K {{k^2}} } }}, 其它}\\ {{C_t}-{C_{t-1}}, t \ge Q-K} \end{array}} \right. $ |
式中，dt表示第t個一階差分；Ct表示第t個倒譜系數；Q表示倒譜系數的階數；K表示一階導數的時間差，可取1或2。
1.2 反向MFCC提取
反向MFCC的基本思想是在反向Mel頻率刻度上放置均勻的濾波器組，從而使得高頻區域的濾波器分布比低頻區域密集。反向Mel頻率刻度fI-mel與實際頻率之間的轉換如式(7)所示[21]：
| $ {f_{I-mel}}=2195.29-2595{\log _{10}}\left({1 + \frac{{4031.25-g}}{{700}}} \right) $ |
其中，f為實際頻率，單位為Hz。
記反向Mel濾波器響應為，與均勻Mel濾波器相應Hm(k)對應關系如下：
| $ \overline {{H_m}} \left(k \right)={H_{M-m + 1}}\left({\frac{N}{2}-k + 1} \right) $ |
用取代傳統MFCC提取流程中的Hm(k)，即可得到反向MFCC。
反向Mel濾波器組的結構如圖 4所示，可見，與圖 3中的Mel濾波器組相反，反向Mel濾波器組在高頻區域放置較多的濾波器，分布較密集；在低頻區域分布較稀疏。
圖4
反向Mel濾波器組結構
Figure4.
Structure of reverse Mel filter bank
1.3 改進反向MFCC提取
當采樣頻率為8 000 Hz時，在咳嗽信號的整個[0, 4 000]Hz頻域范圍內，0 Hz和4 000 Hz附近頻段的頻譜能量很小[18-19]，低能量頻段包含特征信息少且容易受噪聲影響，咳嗽的特征主要集中在頻譜能量較高的頻段。因此，在進行咳嗽的干濕性分類時，為了最大程度地反映干性咳嗽和濕性咳嗽的區別，有效的特征提取應集中在兩類咳嗽頻譜能量高且差別較顯著的頻段進行。
雖然反向MFCC濾波器組有效地獲取了高頻區域的頻譜信息，但濾波器過多地集中在了咳嗽頻譜能量較小的4 000 Hz附近的頻段上，仍不能滿足咳嗽分類的要求。本文根據咳嗽頻譜能量分布對反向MFCC濾波器組進行改進，在咳嗽頻譜能量較大的[fl, fh]Hz的頻率范圍放置反向Mel濾波器組，其中fl和fh分別對應于高頻譜能量頻段的下限頻率和上限頻率，且fl > 0 Hz，fh < 4 000 Hz。
相對于反向Mel濾波器組的響應，計算改進反向Mel濾波器組響應的區別主要在于中心頻率f′[m]的計算。改進反向Mel濾波器組響應的計算步驟如下：
(1)?根據實際頻率與Mel頻率的對應關系fmel=2 595log10(1+f/700)，得到[fl, fh]Hz對應的Mel頻率范圍為[fmel-l, fmel-h]；
(2)?將Mel頻率范圍[fmel-l, fmel-h]平均劃分為(M+1)份，M為濾波器個數，得到Mel刻度下第m個濾波器的中心頻率fMel[m]=m(fmel-h-fmel-l)/M；
(3)?將Mel刻度的中心頻率fMed[m]轉為以Hz為單位，得到改進反向Mel濾波器的中心頻率f′[m]。
(4)?將f′[m]代入公式(3)和(8)計算得到的，即為改進反向Mel濾波器組的響應。
(5)?用改進反向Mel濾波器組的響應取代傳統MFCC提取流程中的Hm(k)，得到改進反向MFCC。
文中取M=20, fl=55 Hz，fh=3 040 Hz，得到改進的反向Mel濾波器組結構如圖 5所示。可見，改進后的反向濾波器組集中分布在[55, 3 040] Hz頻率范圍的中、高頻區域，即特征提取集中在兩類咳嗽頻譜能量高且差別較顯著的頻段進行。
圖5
改進的反向Mel濾波器組結構
Figure5.
Structure of improved reverse Mel filter bank
2 實驗方法和結果
2.1 數據來源
實驗樣本來自廣州醫學院第一附屬醫院病房，使用DELL原裝機的聲卡信號采集設備進行單聲道采樣，采樣頻率為8 000 Hz，每個錄音保存為16位數字量化的WAV文件。選取不同性別、年齡和呼吸感染情況的患者13例(其中以干性咳嗽為主的患者7例，男4例，女3例，平均年齡53歲；以濕性咳嗽為主的患者例數6例，男3例，女3例，平均年齡59歲), 在患者知情同意的情況下，進行連續、長時間的日常錄音。將錄音人工分割成單個咳嗽信號并在醫生指導下分類，作為訓練樣本和測試樣本。
2.2 隱馬爾可夫模型咳嗽分類
隱馬爾可夫模型由于其優異的動態建模能力在語言識別中獲得了廣泛的應用[22]，本文采用隱馬爾可夫模型來實現干性咳嗽和濕性咳嗽嗽的分類。分別取60個干性咳嗽(其中男32個，女28個；無噪聲環境下46個，有噪聲環境下14個)和60個濕性咳嗽(其中男36個，女24個；無噪聲環境下41個，噪聲環境下19個)作為隱馬爾可夫模型訓練樣本，再分別取120個干性咳嗽(其中男63個，女57個；無噪聲環境下96個，有噪聲環境下24個)和85個濕性咳嗽(其中男49個，女36個；無噪聲環境下61個，有噪聲環境下24個)作為測試樣本，測試樣本不包含訓練樣本。
采用隱馬爾可夫模型進行咳嗽干濕性分類的步驟如下：
(1)?為干性咳嗽信號和濕性咳嗽信號分別定義一個隱馬爾可夫模型，記為H={h1, h2}。采用含2個隱含狀態的自左向右隱馬爾可夫模型來實現，每個狀態為包含3個高斯元的混合模型，采用C均值的方法初始化高斯混合模型參數。
(2)?取60個干性咳嗽和60個濕性咳嗽作為訓練樣本，分別訓練hi(i=1, 2)，獲得最佳模型μi。采用Baum-Welch算法[22]進行訓練，迭代次數設置為50次。
(3)?采用訓練好的模型μi(i=1, 2)進行咳嗽分類。采用Viterbi算法[22]對未知的輸入信號O，計算Pr(O｜μi)(i=1, 2)，若Pr(O｜μi)=max(Pr(O, μi), i=1, 2))，則將輸入信號O分類為hj類。
2.3 實驗結果
分別采用包含一階和二級差分系數的39維傳統MFCC、反向MFCC和改進反向MFCC作為特征量，采用隱馬爾可夫模型進行咳嗽的干濕性分類，結果如表 1所示，其中干性咳嗽的分類準確率=正確分類的干性咳嗽個數/干性咳嗽總數；濕性咳嗽的分類準確率=正確分類的濕性咳嗽個數/濕性咳嗽總數；總體分類準確率=正確分類的咳嗽個數/咳嗽總數。
由表 1可見，三種特征提取方法中，反向MFCC對應的干性咳嗽、濕性咳嗽和總體分類準確率都最低，而改進反向MFCC相應的三種分類準確率都最高。另外采用三種不同的特征時，濕性咳嗽分類準確率皆低于干性咳嗽，反映了濕性咳嗽的頻譜構成更復雜和多樣化，這與文獻[19]報道的咳嗽特點一致。
表1
三種特征提取方法的分類結果比較
Table1.
Comparison of cough classification results respectively by the three feature extraction methods
3 討論
在咳嗽識別和咳嗽分類研究中，以往的研究大多參考語音識別系統，采用傳統的MFCC作為特征量，傳統MFCC在Mel刻度放置均勻濾波器，特征提取集中在信號的低頻和中頻區域進行。然而通過分析咳嗽的機制和聲學特征發現，在咳嗽信號的整個0～4 000 Hz頻域范圍內，0 Hz和4 000 Hz附近頻段的頻譜能量很小，咳嗽的頻譜能量主要分布在50～3 050 Hz范圍內，且濕性咳嗽由于呼吸道中有痰導致氣流震蕩加劇，其頻譜的中頻和高頻成分比干性咳嗽更加豐富。咳嗽分類的目標是區分兩類咳嗽，因此有效的特征提取應集中在兩類咳嗽頻譜能量高且差別較顯著的頻段進行。本文提出的改進反向MFCC，將濾波器組集中放置在兩類咳嗽頻譜能量高且差別較顯著的頻段，實現了有效的特征提取，從而獲得了優于傳統MFCC的分類結果，使得總體分類準確率從89.76%提高到了93.66%；該方法主要提取了咳嗽中頻和高頻段的特征，使得該頻段特征明顯的濕性咳嗽更容易被正確識別，因而采用該方法時濕性咳嗽分類準確率的提升比干性咳嗽顯著。另一方面，與傳統MFCC一樣，改進反向MFCC仍是反映信號頻譜能量的特征表示方法，因而容易受到噪聲的影響，當干性咳嗽疊加了噪聲，尤其是中頻和高頻的噪聲時，由于中頻和高頻成分增多，容易被錯分類為濕性咳嗽。如何進一步提高特征提取方法的魯棒性是我們將來的研究重點。
