本文通過分析嗓音的發音機制,提取正常與病態嗓音的傳統聲學參數[基頻、Mel倒譜系數(MFCC)、線性預測系數(LPCC)、頻率微擾、振幅微擾]與非線性動力學特征參數[熵(樣本熵、模糊熵、多尺度熵)、計盒維數、計維截距和Hurst參數],作為病態嗓音識別的特征矢量集。應用支持向量機(SVM)對/a/音的78例正常嗓音與73例病態嗓音和/i/音的78例正常嗓音與80例病態嗓音進行建模與識別。結果表明,相對于傳統的聲學特征參數,非線性特征參數能更好地區分正常與病態嗓音;實驗提取的所有參數中,除了多尺度熵,/a/音的正常與病態嗓音的識別率均高于/i/音,因此為了達到識別病態嗓音的目的,國內外相關研究大多采用/a/音數據;多尺度熵特征對/i/音的正常與病態嗓音的識別率較/a/音高,它或能為評價聲帶代償功能狀態的研究提供有益的啟發。
引用本文： 甘德英, 胡維平, 趙冰心. 傳統聲學特征和非線性特征用于病態嗓音的比較研究. 生物醫學工程學雜志, 2014, 31(5): 1149-1154. doi: 10.7507/1001-5515.20140217 復制
引言
嗓音的產生主要靠聲帶振動，因此聲帶的質量、韌性、黏膜病變等因素均會影響其振動的頻率、振幅和模式[1]。嗓音的聲學檢測是喉功能檢查的客觀手段。發聲是喉的重要功能之一,喉部疾病往往出現發聲障礙，故嗓音檢查一直為喉科醫生、言語病理學家所重視[2]。通過分析嗓音的發音機制，提取嗓音信號的Mel頻率倒譜系數(Mel frequency cepstral coefficient,MFCC)、線性預測倒譜系數(linear prediction cepstrum coefficient,LPCC)、基音頻率(fundamental frequency,F0)、頻率微擾及振幅微擾等傳統聲學特征,可對嗓音客觀地進行評價[3]。但是語音信號是一種非線性、非平穩的信號，而傳統聲學特征參數的提取都是將語音信號近似為線性信號，因此在檢測病態嗓音時忽略了語音信號的非線性特性，存在一定的局限性[4]。近年來，非線性動力學技術在病態嗓音檢測中得到廣泛使用，研究者致力于探尋能有效描述嗓音病理特征的非線性聲學特征參數，以彌補傳統聲學特征參數的不足，同時提高評估病態嗓音的能力，為臨床聲學檢測提供可靠的理論基礎[6]。
聲帶正常振動時，通過調節聲帶本體層張力改變其基頻，發出不同音調的正常聲音；聲帶受到病變影響而造成嗓音嘶啞和發聲困難時，通過拉緊聲帶的本體層提供額外的張力以達到恢復控制聲帶的目的[7]。發/a/音時聲帶本體層自然放松，而發/i/音時聲帶本體層張緊，聲帶拉緊變薄，張力增加，比發/a/音更易受到代償能力的影響，掩蓋嗓音的病理特征，導致其病態嗓音的檢出率相對/a/音較低[8]。本文擬將單特征應用于病態嗓音/a/音、/i/音的識別比較，從而進行聲帶代償功能的探索性研究。
1 數據來源
本文試驗數據來源于臨床病例，采集數據時的環境要求在安靜的室內進行；采樣頻率為16 kHz，時間1.5~3 s；受試聲樣為漢語元音/a/及/i/，分別對正常人和各類喉科疾病的患者進行語音采樣。正常對照組/a/音78例、/i/音78例，被采集者年齡18~40周歲，平均年齡25周歲，經詢問近期無喉部疾病；病態嗓音組/a/音73例、/i/音80例，被采集者年齡15~50周歲，平均年齡27周歲，為前來醫院就診的患者。試驗時，分別計算數據庫中/a/和/i/的正常嗓音與病態嗓音的各種聲學特征參數值，隨后由支持向量機(support vector machine,SVM)自動分配正常嗓音和病態嗓音的訓練集與測試集，分別對/a/音、/i/音進行識別[9]。
2 特征提取
2.1 傳統聲學特征參數的提取
用于嗓音分析與識別的傳統聲學特征參數有很多，并在嗓音的病理檢測與診斷方面取得了一定的成果。本試驗選取了在病態嗓音臨床醫學檢測中廣泛使用的F0、頻率微擾、振幅微擾、MFCC、LPCC[10]。當聲帶發生病變時，這些特征參數值將會偏離正常值，從而有效識別正常嗓音與病態嗓音。
2.1.1 Mel頻域倒譜系數與線性預測倒譜系數
MFCC的分析著眼于人耳的聽覺特性，因為人耳聽到聲音的高低與聲音的頻率并不成線性正比關系，而用Mel頻率尺度更符合人耳的聽覺特性。Mel頻率與實際頻率的具體關系式如下[11]：
Mel(f)=2595lg(1+f/700)
LPCC由于利用了線性預測中聲道系統函數的最小相位特性，避免了相位卷繞問題，并且運算量小，含有語音信號頻譜的包絡信息。MFCC與LPCC的計算方法見文獻[11]。
2.1.2 基音頻率、頻率微擾與振幅微擾
F0是指聲帶振動的頻率，在發音系統疾病診斷中應用廣泛。本試驗采用自相關法提取嗓音信號的F0，根據計算得到的F0及其對應的振幅可計算頻率微擾與振幅微擾參數，其中頻率微擾參數包括絕對頻率抖動(Jita)、頻率抖動百分比(Jitt)、頻率微擾商(pitch period perturbation quotient，PPQ)，振幅微擾包括振幅抖動百分比(Shim)、振幅抖動(ShdB)、振幅微擾商(amplitude perturbation quotient，APQ)[12]。文獻[12]驗證了這些韻律特征能有效地描述病態嗓音的病理特征。
2.2 熵的提取
2.2.1 樣本熵
樣本熵計算步驟[13]：
(1)給定長度為N的時間序列{u(i),1≤i≤N},按式Xim={u(i),u(i+1),…,u(i+m-1)}重構m維向量Xim，1≤i≤N-m+1,Xim代表從第i個矢量開始的連續m個u值，m稱為嵌入維數。
(2)定義Xim與Xjm(1≤j≤N-m且j≠i)間的距離dijm為兩向量對應元素中差值最大的一個，即
dijm=d[Xim,Xjm]=maxk∈(0,m-1)|u(i+k)-u(j+k)|, 式中1≤i≤N-m,1≤j≤N-m且j≠i。
(3)給定閾值r，對每個1≤i≤N-m的值，計算dijm小于r的數目(模版匹配數)與距離總數(矢量總個數)的比值，記做：Brm(i)=(N-m-1)-1·(dijm-r)，Θ是Heaviside函數，。
(4)對所有的Brm(i)求平均，記為Bm(r)=(N-m)-1·。
(5)增加一維，對m+1維矢量繼續重復步驟(1)~(4)，得到Bm+1(r)。
(6)理論上，此序列的樣本熵為：SampEn(m,r,N)=}。
當N為有限值時，樣本熵的估計值為： SampEn(m,r,N)=。
根據文獻[13]和[14]，一般情況下，m=1或2，r=0.1SD～0.25SD(SD為原始序列的標準差)，計算得到的樣本熵具有較合理的統計特性。
2.2.2 多尺度熵
2002年，Costa等[15]提出多尺度熵(multiscale entropy,MSE)的概念，構造粗粒化時間序列{yτ}，yj(τ)=。式中τ是尺度因子，N是原始序列長度,可計算不同尺度上的樣本熵。多尺度熵粗粒化過程實質是頻率篩分的過程，隨尺度增大，頻率逐漸降低，考慮了信號在不同尺度上的低頻信息。
多尺度熵定義為樣本熵值在多個尺度下的集合[16]：
| $\begin{align} & MSE=\{\tau Sampen(\tau ,m,r)= \\ & -ln[C\tau ,m+1\left( r \right)/C\tau ,m\left( r \right)] \\ \end{align}$ |
2.2.3 模糊熵
根據樣本熵的定義，矢量間的相似性是基于Heaviside函數。Heaviside函數的主要特征是只用是0或1的函數判斷就確定輸入信號是否在閾值內。這個函數是不連續的，因為當值接近0或1它會發生突變。在模糊熵[17-18]中，引用了模糊隸屬函數u(dijm,r)來模糊測量序列Xim與Xjm的相似性。
通常高斯函數、Sigmoid函數、鐘形函數或其它模糊隸屬函數都可做為兩個向量間的相似性度量，模糊隸屬函數要滿足兩條性質：① 連續函數；② 凸函數[18]。
2.2.4 香農熵與二階Renyi熵
香農熵：假設將一個系統分成p個單元格，pi表示某一時刻系統軌道在第i格的概率，則香農熵的定義為[19]
| $H=-\sum\limits_{i}{{{p}_{i}}ln{{p}_{i}}}$ |
Renyi熵：基于轉移概率的Renyi熵能夠較好地度量信息的損失量。Pi1,i2,…,im為任意時刻n觀察值落入單元格i1及時刻n+1觀察值落入間隔i2等的聯合概率，m維r-Renyi熵的計算表達式為[20]
| ${{H}_{r}}\left( m \right)=\frac{1}{1-r}ln\sum\limits_{{{i}_{1}},{{i}_{2}},\ldots ,{{i}_{m}}}{{{p}_{{{i}_{1}}^{r},{{i}_{2}},\ldots {{i}_{m}}}}}~r0,r\ne 1$ |
Renyi熵是香農熵的擴展，當r=1時，r-Renyi熵為香農熵。
2.3 非線性動力學特征的提取
2.3.1 計盒維數與計維截距
語音信號是一個復雜的非線性過程，其中存在著產生混沌的機制[21]。分形理論是描述混沌信號的一種有效手段，是一種分析非線性動態系統的新方法[22]。分形的度量有多種描述方法，常用的有計盒維數、容量維數、關聯維數、Hausdorff維等[23]。本文采用計盒維數的描述方法。
計盒維數：N維歐氏空間子集F的計盒維數Db定義為[24]
| ${{D}_{b}}=\underset{\sigma \to 0}{\mathop{lim}}\,\frac{logN\sigma \left( F \right)}{log(1/\sigma )}$ |
假定上述極限存在，其中Nσ(F)在本文中定義為：覆蓋F的邊長為σ的最小立方體個數。
本實驗中采用計算量較小的多點擬合來計算Db。具體過程如下：
(1)取M個σ值，記σi(i=1，…，M);
(2)計算Nσi(F)。先將語音信號歸一化到單位正方形區域[21]，即
F={s(t)|min(s(t))=0,max(s(t))=1,0≤t≤1}。s(t)表示歸一化后的語音信號，取最小正方形的邊長為σ，易知在區間In=[(n-1)σ,nσ]內，覆蓋F的最小正方形數為
Nn=[(max s(t)-min s(t))/σ],(t∈In) ,
式中[X]表示不小于X的最小正整數，因此在[0, 1]區間內覆蓋F的總正方形個數為：
| ${{N}_{\sigma }}\left( F \right)=\sum\limits_{n=1}^{1/\sigma }{{{N}_{n}}}$ |
(3)令xi=log(1/σ)，yi=logNσ(F)，分別計算點(xi,yi)，按最小均方誤差擬合一條直線y=kx+b,則此直線的斜率k為計盒維數Db,b為截距，即：
| ${{D}_{b}}=k=\frac{(\sum\limits_{i=1}^{M}{{{y}_{i}}})(\sum\limits_{i=1}^{M}{{{x}_{i}}})-M(\sum\limits_{i=1}^{M}{{{y}_{i}}{{x}_{i}}})}{{{(\sum\limits_{i=1}^{M}{{{x}_{i}}})}^{2}}-M\sum\limits_{i=1}^{M}{{{x}_{i}}^{2}}}$ |
2.3.2 Hurst參數
聲帶的病變會導致聲音產生過程中混沌現象的出現，而嗓音信號在廣泛的時間尺度上顯示出自相似特性，自相似過程的特征之一即為呈現出長相關性[25]。Hurst參數可以描述嗓音信號自相似的長相關特性的程度，因此本實驗采用Hurst參數檢測喉部病理狀況。
設樣本嗓音信號時間序列為X(i)，將時間序列分成Y個大小為m的非重疊區塊，對于每一個非重疊區塊的R/S統計量定義為
| $\begin{align} & \frac{R}{S}\left( m \right)=\frac{1}{{{S}^{2}}\left( m \right)}\left[ \underset{0\le p\le m}{\mathop{max}}\,\left( Y\left( p \right)-\frac{p}{m}Y\left( m \right) \right) \right. \\ & \left. -\underset{0\le p\le m}{\mathop{min}}\,\left( Y\left( p \right)-\frac{p}{m}Y\left( m \right) \right) \right], \\ \end{align}$ |
式中，S2(m)為樣本方差，則Hurst參數值為log(R/m)/S(m))與log(m)的最小二乘擬合值[26]。
3 實驗及結果
本試驗使用的數據含/a/音正常嗓音78例，病態嗓音73例，/i/音正常嗓音78例，病態嗓音80例。每例數據取前5 000個采樣點進行特征提取和模式識別。使用MATLAB中crossvalind函數，得到訓練集。/a/音中正常嗓音39例和病態嗓音37例作為訓練集，其余的作為測試集；/i/音中正常嗓音39例和病態嗓音40例作為訓練集，其余的作為測試集。隨機選擇100次訓練集和對應的測試集，可得到這100次SVM識別結果的最大識別率和平均識別率。SVM的核函數選擇徑向基核函數，參數為MATLAB默認值[27]。
傳統聲學特征參數的提取及識別： 首先將語音作分幀處理，幀長512，幀移256，提取每幀語音的MFCC、LPCC、F0、頻率微擾(Jitt、Jita、PPQ)、振幅微擾(Shim、ShdB、APQ)。采用SVM分別對/a/音和/i/音的上述九種特征參數進行訓練和識別，得到各自正常與病態嗓音的識別率如表 1所示。
表1
/a/和/i/傳統聲學特征參數的SVM識別率
Table1.
Recognition rates of SVM on traditional acoustic parameters for /a/ and /i/
表 1中，對應相同的傳統聲學特征，/a/音的嗓音識別率全部高于/i/音，采用基音頻率對/a/音進行訓練與識別時取得最高識別率90.67%，/i/音的最高識別率為82.28%；采用MFCC聲學特征參數/a/音最高識別率為88.00%，/i/音為84.81%。實驗結果表明傳統聲學特征參數對病態和正常嗓音具有較好的可區分度。
計算熵特征時，設置參數m=2,r=0.25SD，τ=2。計算/a/音78例正常嗓音和73例病態嗓音與/i/音78例正常嗓音和80例病態嗓音的的樣本熵、多尺度熵、模糊熵、香農熵及第二階Renyi熵。計算模糊熵時，采用隸屬函數u(dijm,r)=exp(-dij2/r)。提取的非線性動力學參數有計盒維數、計維截距及Hurst參數。
圖 1為 /a/音和/i/音正常與病態嗓音的熵與非線性動力學特征值分布，矩形框下邊界為下四分位數，上邊界為上四分位數，觸須的上、下邊界分別為特征參數的最大值和最小值。矩形框內的橫線為特征參數的均值，矩形框外“+”表示離群點。a_nor: /a/音正常嗓音；a_ab: /a/音病態嗓音；i_nor: /i/音正常嗓音；i_ab: /i/音病態嗓音。
圖1
/a/音和/i/音正常與病態嗓音的特征參數值分布
Figure1.
Data distributions of normal and pathological voice for /a/ and /i/
模糊熵特征分布圖中，/a/音與/i/音正常嗓音參數分布的上邊界與病態嗓音參數分布的下邊界無重合部分并相離較遠，正常與病態特征數據區分度較明顯，表明模糊熵對/a/音與/i/音均有較好的區分效果，但模糊熵對/a/音正常與病態嗓音的區分度較/i/音更為明顯。
Hurst參數描述混沌時間序列在時間尺度上的自相似，系統越混沌，Hurst值越大，因此Hurst參數特征值表明正常嗓音的相關程度較高，病態嗓音信號具有較強的混沌特性，信號較復雜；/i/音較/a/音更復雜，并具較強的混沌特性。Hurst參數值能有效區分正常與病態嗓音，但對/i/音的可區分度較/a/音小。此外，計盒維數、計維截距、樣本熵及第二階Renyi熵對正常與病態嗓音也具有較好的區分度，但不如模糊熵與Hurst參數。
/a/音和/i/音非線性聲學特征參數的SVM識別率如表 2所示，與表 1中/a/音和/i/音傳統聲學特征參數的SVM識別率比較，傳統聲學單個特征的平均識別率分布在50.87%~80.32%和48.37%~73.81%,而非線性聲學單特征的平均識別率分布在65.16%~87.58%和57.86%~76.39%，這說明采用非線性聲學特征參數得到的正常與病態嗓音的識別率明顯高于傳統聲學特征參數。采用模糊熵、Hurst參數兩種特征值識別/a/音和/i/音的正常與病態嗓音均取得較高的識別率，表明這兩種非線性聲學特征對病態與正常嗓音具有較好的區分度，而香農熵、第二階Renyi熵、樣本熵、計盒維數及計維截距亦能較好區分正常與病態嗓音，但不如模糊熵和Hurst參數，這與前面特征值分布分析的結果一致。此外，除了多尺度熵，其余各參數對/a/音的正常、病態嗓音的識別率均高于/i/音，這與采用傳統聲學特征參數得到的識別結果一致，因此為了達到識別病態嗓音的目的，目前國內外相關研究幾乎都采用/a/音數據，而不采用/i/音數據。
表2
/a/和/i/非線性聲學特征參數的SVM識別率
Table2.
Recognition rate of SVM on nonlinear characteristic parameters for /a/ and /i/
發/a/音時聲帶本體層自然放松其F0相對較低，而發/i/音時聲帶本體層張緊F0相對較高。因此，發/i/音時，環甲肌位置決定了其相對張緊的聲帶本體層能夠提供甲杓肌更大的代償空間，正是這種相對較高程度的代償功能的介入掩蓋了病態嗓音的特征，導致其病態嗓音的檢出率相對/a/音較低。當采用多尺度熵特征值時，可提取/i/音不同尺度的低頻信息，把隱藏在拉緊聲帶中的病理嗓音特征有效地提取出來，得到的識別率較/a/音高。多尺度熵以及/a/音、/i/音非線性特征的相對關系或能為研究聲帶代償功能的狀態提供有益的啟發,而對病態嗓音聲帶代償功能的狀態進行研究正是我們下一步的研究目標。
4 結束語
近年來，非線性特征成為描述語音信號的一個有力工具。本試驗通過比較傳統聲學特征參數和非線性聲學特征參數對/a/音及/i/音的正常與病態嗓音的SVM識別率，充分證明非線性聲學特征能更好地描述嗓音信號的病理特征，并能有效區分正常與病態嗓音。尤其是模糊熵對嗓音的識別，在/a/音取得94.67%的最高識別率，/i/音取得87.58%的最高識別率。除多尺度熵外，所有聲學特征參數對/a/音的識別率均高于/i/音，因此為了達到識別病態嗓音的目的，國內外的相關研究大多采用/a/音數據。多尺度熵對/i/音的嗓音識別率較/a/音高，表明多尺度熵能夠把隱藏在拉緊聲帶中的病理嗓音特征有效地提取出來，而發/i/音比發/a/音更易受到聲帶代償能力的影響，因此多尺度熵以及/a/音、/i/音非線性特征的相對關系或能為評價聲帶代償功能狀態的研究提供有益的啟發。
引言
嗓音的產生主要靠聲帶振動，因此聲帶的質量、韌性、黏膜病變等因素均會影響其振動的頻率、振幅和模式[1]。嗓音的聲學檢測是喉功能檢查的客觀手段。發聲是喉的重要功能之一,喉部疾病往往出現發聲障礙，故嗓音檢查一直為喉科醫生、言語病理學家所重視[2]。通過分析嗓音的發音機制，提取嗓音信號的Mel頻率倒譜系數(Mel frequency cepstral coefficient,MFCC)、線性預測倒譜系數(linear prediction cepstrum coefficient,LPCC)、基音頻率(fundamental frequency,F0)、頻率微擾及振幅微擾等傳統聲學特征,可對嗓音客觀地進行評價[3]。但是語音信號是一種非線性、非平穩的信號，而傳統聲學特征參數的提取都是將語音信號近似為線性信號，因此在檢測病態嗓音時忽略了語音信號的非線性特性，存在一定的局限性[4]。近年來，非線性動力學技術在病態嗓音檢測中得到廣泛使用，研究者致力于探尋能有效描述嗓音病理特征的非線性聲學特征參數，以彌補傳統聲學特征參數的不足，同時提高評估病態嗓音的能力，為臨床聲學檢測提供可靠的理論基礎[6]。
聲帶正常振動時，通過調節聲帶本體層張力改變其基頻，發出不同音調的正常聲音；聲帶受到病變影響而造成嗓音嘶啞和發聲困難時，通過拉緊聲帶的本體層提供額外的張力以達到恢復控制聲帶的目的[7]。發/a/音時聲帶本體層自然放松，而發/i/音時聲帶本體層張緊，聲帶拉緊變薄，張力增加，比發/a/音更易受到代償能力的影響，掩蓋嗓音的病理特征，導致其病態嗓音的檢出率相對/a/音較低[8]。本文擬將單特征應用于病態嗓音/a/音、/i/音的識別比較，從而進行聲帶代償功能的探索性研究。
1 數據來源
本文試驗數據來源于臨床病例，采集數據時的環境要求在安靜的室內進行；采樣頻率為16 kHz，時間1.5~3 s；受試聲樣為漢語元音/a/及/i/，分別對正常人和各類喉科疾病的患者進行語音采樣。正常對照組/a/音78例、/i/音78例，被采集者年齡18~40周歲，平均年齡25周歲，經詢問近期無喉部疾病；病態嗓音組/a/音73例、/i/音80例，被采集者年齡15~50周歲，平均年齡27周歲，為前來醫院就診的患者。試驗時，分別計算數據庫中/a/和/i/的正常嗓音與病態嗓音的各種聲學特征參數值，隨后由支持向量機(support vector machine,SVM)自動分配正常嗓音和病態嗓音的訓練集與測試集，分別對/a/音、/i/音進行識別[9]。
2 特征提取
2.1 傳統聲學特征參數的提取
用于嗓音分析與識別的傳統聲學特征參數有很多，并在嗓音的病理檢測與診斷方面取得了一定的成果。本試驗選取了在病態嗓音臨床醫學檢測中廣泛使用的F0、頻率微擾、振幅微擾、MFCC、LPCC[10]。當聲帶發生病變時，這些特征參數值將會偏離正常值，從而有效識別正常嗓音與病態嗓音。
2.1.1 Mel頻域倒譜系數與線性預測倒譜系數
MFCC的分析著眼于人耳的聽覺特性，因為人耳聽到聲音的高低與聲音的頻率并不成線性正比關系，而用Mel頻率尺度更符合人耳的聽覺特性。Mel頻率與實際頻率的具體關系式如下[11]：
Mel(f)=2595lg(1+f/700)
LPCC由于利用了線性預測中聲道系統函數的最小相位特性，避免了相位卷繞問題，并且運算量小，含有語音信號頻譜的包絡信息。MFCC與LPCC的計算方法見文獻[11]。
2.1.2 基音頻率、頻率微擾與振幅微擾
F0是指聲帶振動的頻率，在發音系統疾病診斷中應用廣泛。本試驗采用自相關法提取嗓音信號的F0，根據計算得到的F0及其對應的振幅可計算頻率微擾與振幅微擾參數，其中頻率微擾參數包括絕對頻率抖動(Jita)、頻率抖動百分比(Jitt)、頻率微擾商(pitch period perturbation quotient，PPQ)，振幅微擾包括振幅抖動百分比(Shim)、振幅抖動(ShdB)、振幅微擾商(amplitude perturbation quotient，APQ)[12]。文獻[12]驗證了這些韻律特征能有效地描述病態嗓音的病理特征。
2.2 熵的提取
2.2.1 樣本熵
樣本熵計算步驟[13]：
(1)給定長度為N的時間序列{u(i),1≤i≤N},按式Xim={u(i),u(i+1),…,u(i+m-1)}重構m維向量Xim，1≤i≤N-m+1,Xim代表從第i個矢量開始的連續m個u值，m稱為嵌入維數。
(2)定義Xim與Xjm(1≤j≤N-m且j≠i)間的距離dijm為兩向量對應元素中差值最大的一個，即
dijm=d[Xim,Xjm]=maxk∈(0,m-1)|u(i+k)-u(j+k)|, 式中1≤i≤N-m,1≤j≤N-m且j≠i。
(3)給定閾值r，對每個1≤i≤N-m的值，計算dijm小于r的數目(模版匹配數)與距離總數(矢量總個數)的比值，記做：Brm(i)=(N-m-1)-1·(dijm-r)，Θ是Heaviside函數，。
(4)對所有的Brm(i)求平均，記為Bm(r)=(N-m)-1·。
(5)增加一維，對m+1維矢量繼續重復步驟(1)~(4)，得到Bm+1(r)。
(6)理論上，此序列的樣本熵為：SampEn(m,r,N)=}。
當N為有限值時，樣本熵的估計值為： SampEn(m,r,N)=。
根據文獻[13]和[14]，一般情況下，m=1或2，r=0.1SD～0.25SD(SD為原始序列的標準差)，計算得到的樣本熵具有較合理的統計特性。
2.2.2 多尺度熵
2002年，Costa等[15]提出多尺度熵(multiscale entropy,MSE)的概念，構造粗粒化時間序列{yτ}，yj(τ)=。式中τ是尺度因子，N是原始序列長度,可計算不同尺度上的樣本熵。多尺度熵粗粒化過程實質是頻率篩分的過程，隨尺度增大，頻率逐漸降低，考慮了信號在不同尺度上的低頻信息。
多尺度熵定義為樣本熵值在多個尺度下的集合[16]：
| $\begin{align} & MSE=\{\tau Sampen(\tau ,m,r)= \\ & -ln[C\tau ,m+1\left( r \right)/C\tau ,m\left( r \right)] \\ \end{align}$ |
2.2.3 模糊熵
根據樣本熵的定義，矢量間的相似性是基于Heaviside函數。Heaviside函數的主要特征是只用是0或1的函數判斷就確定輸入信號是否在閾值內。這個函數是不連續的，因為當值接近0或1它會發生突變。在模糊熵[17-18]中，引用了模糊隸屬函數u(dijm,r)來模糊測量序列Xim與Xjm的相似性。
通常高斯函數、Sigmoid函數、鐘形函數或其它模糊隸屬函數都可做為兩個向量間的相似性度量，模糊隸屬函數要滿足兩條性質：① 連續函數；② 凸函數[18]。
2.2.4 香農熵與二階Renyi熵
香農熵：假設將一個系統分成p個單元格，pi表示某一時刻系統軌道在第i格的概率，則香農熵的定義為[19]
| $H=-\sum\limits_{i}{{{p}_{i}}ln{{p}_{i}}}$ |
Renyi熵：基于轉移概率的Renyi熵能夠較好地度量信息的損失量。Pi1,i2,…,im為任意時刻n觀察值落入單元格i1及時刻n+1觀察值落入間隔i2等的聯合概率，m維r-Renyi熵的計算表達式為[20]
| ${{H}_{r}}\left( m \right)=\frac{1}{1-r}ln\sum\limits_{{{i}_{1}},{{i}_{2}},\ldots ,{{i}_{m}}}{{{p}_{{{i}_{1}}^{r},{{i}_{2}},\ldots {{i}_{m}}}}}~r0,r\ne 1$ |
Renyi熵是香農熵的擴展，當r=1時，r-Renyi熵為香農熵。
2.3 非線性動力學特征的提取
2.3.1 計盒維數與計維截距
語音信號是一個復雜的非線性過程，其中存在著產生混沌的機制[21]。分形理論是描述混沌信號的一種有效手段，是一種分析非線性動態系統的新方法[22]。分形的度量有多種描述方法，常用的有計盒維數、容量維數、關聯維數、Hausdorff維等[23]。本文采用計盒維數的描述方法。
計盒維數：N維歐氏空間子集F的計盒維數Db定義為[24]
| ${{D}_{b}}=\underset{\sigma \to 0}{\mathop{lim}}\,\frac{logN\sigma \left( F \right)}{log(1/\sigma )}$ |
假定上述極限存在，其中Nσ(F)在本文中定義為：覆蓋F的邊長為σ的最小立方體個數。
本實驗中采用計算量較小的多點擬合來計算Db。具體過程如下：
(1)取M個σ值，記σi(i=1，…，M);
(2)計算Nσi(F)。先將語音信號歸一化到單位正方形區域[21]，即
F={s(t)|min(s(t))=0,max(s(t))=1,0≤t≤1}。s(t)表示歸一化后的語音信號，取最小正方形的邊長為σ，易知在區間In=[(n-1)σ,nσ]內，覆蓋F的最小正方形數為
Nn=[(max s(t)-min s(t))/σ],(t∈In) ,
式中[X]表示不小于X的最小正整數，因此在[0, 1]區間內覆蓋F的總正方形個數為：
| ${{N}_{\sigma }}\left( F \right)=\sum\limits_{n=1}^{1/\sigma }{{{N}_{n}}}$ |
(3)令xi=log(1/σ)，yi=logNσ(F)，分別計算點(xi,yi)，按最小均方誤差擬合一條直線y=kx+b,則此直線的斜率k為計盒維數Db,b為截距，即：
| ${{D}_{b}}=k=\frac{(\sum\limits_{i=1}^{M}{{{y}_{i}}})(\sum\limits_{i=1}^{M}{{{x}_{i}}})-M(\sum\limits_{i=1}^{M}{{{y}_{i}}{{x}_{i}}})}{{{(\sum\limits_{i=1}^{M}{{{x}_{i}}})}^{2}}-M\sum\limits_{i=1}^{M}{{{x}_{i}}^{2}}}$ |
2.3.2 Hurst參數
聲帶的病變會導致聲音產生過程中混沌現象的出現，而嗓音信號在廣泛的時間尺度上顯示出自相似特性，自相似過程的特征之一即為呈現出長相關性[25]。Hurst參數可以描述嗓音信號自相似的長相關特性的程度，因此本實驗采用Hurst參數檢測喉部病理狀況。
設樣本嗓音信號時間序列為X(i)，將時間序列分成Y個大小為m的非重疊區塊，對于每一個非重疊區塊的R/S統計量定義為
| $\begin{align} & \frac{R}{S}\left( m \right)=\frac{1}{{{S}^{2}}\left( m \right)}\left[ \underset{0\le p\le m}{\mathop{max}}\,\left( Y\left( p \right)-\frac{p}{m}Y\left( m \right) \right) \right. \\ & \left. -\underset{0\le p\le m}{\mathop{min}}\,\left( Y\left( p \right)-\frac{p}{m}Y\left( m \right) \right) \right], \\ \end{align}$ |
式中，S2(m)為樣本方差，則Hurst參數值為log(R/m)/S(m))與log(m)的最小二乘擬合值[26]。
3 實驗及結果
本試驗使用的數據含/a/音正常嗓音78例，病態嗓音73例，/i/音正常嗓音78例，病態嗓音80例。每例數據取前5 000個采樣點進行特征提取和模式識別。使用MATLAB中crossvalind函數，得到訓練集。/a/音中正常嗓音39例和病態嗓音37例作為訓練集，其余的作為測試集；/i/音中正常嗓音39例和病態嗓音40例作為訓練集，其余的作為測試集。隨機選擇100次訓練集和對應的測試集，可得到這100次SVM識別結果的最大識別率和平均識別率。SVM的核函數選擇徑向基核函數，參數為MATLAB默認值[27]。
傳統聲學特征參數的提取及識別： 首先將語音作分幀處理，幀長512，幀移256，提取每幀語音的MFCC、LPCC、F0、頻率微擾(Jitt、Jita、PPQ)、振幅微擾(Shim、ShdB、APQ)。采用SVM分別對/a/音和/i/音的上述九種特征參數進行訓練和識別，得到各自正常與病態嗓音的識別率如表 1所示。
表1
/a/和/i/傳統聲學特征參數的SVM識別率
Table1.
Recognition rates of SVM on traditional acoustic parameters for /a/ and /i/
表 1中，對應相同的傳統聲學特征，/a/音的嗓音識別率全部高于/i/音，采用基音頻率對/a/音進行訓練與識別時取得最高識別率90.67%，/i/音的最高識別率為82.28%；采用MFCC聲學特征參數/a/音最高識別率為88.00%，/i/音為84.81%。實驗結果表明傳統聲學特征參數對病態和正常嗓音具有較好的可區分度。
計算熵特征時，設置參數m=2,r=0.25SD，τ=2。計算/a/音78例正常嗓音和73例病態嗓音與/i/音78例正常嗓音和80例病態嗓音的的樣本熵、多尺度熵、模糊熵、香農熵及第二階Renyi熵。計算模糊熵時，采用隸屬函數u(dijm,r)=exp(-dij2/r)。提取的非線性動力學參數有計盒維數、計維截距及Hurst參數。
圖 1為 /a/音和/i/音正常與病態嗓音的熵與非線性動力學特征值分布，矩形框下邊界為下四分位數，上邊界為上四分位數，觸須的上、下邊界分別為特征參數的最大值和最小值。矩形框內的橫線為特征參數的均值，矩形框外“+”表示離群點。a_nor: /a/音正常嗓音；a_ab: /a/音病態嗓音；i_nor: /i/音正常嗓音；i_ab: /i/音病態嗓音。
圖1
/a/音和/i/音正常與病態嗓音的特征參數值分布
Figure1.
Data distributions of normal and pathological voice for /a/ and /i/
模糊熵特征分布圖中，/a/音與/i/音正常嗓音參數分布的上邊界與病態嗓音參數分布的下邊界無重合部分并相離較遠，正常與病態特征數據區分度較明顯，表明模糊熵對/a/音與/i/音均有較好的區分效果，但模糊熵對/a/音正常與病態嗓音的區分度較/i/音更為明顯。
Hurst參數描述混沌時間序列在時間尺度上的自相似，系統越混沌，Hurst值越大，因此Hurst參數特征值表明正常嗓音的相關程度較高，病態嗓音信號具有較強的混沌特性，信號較復雜；/i/音較/a/音更復雜，并具較強的混沌特性。Hurst參數值能有效區分正常與病態嗓音，但對/i/音的可區分度較/a/音小。此外，計盒維數、計維截距、樣本熵及第二階Renyi熵對正常與病態嗓音也具有較好的區分度，但不如模糊熵與Hurst參數。
/a/音和/i/音非線性聲學特征參數的SVM識別率如表 2所示，與表 1中/a/音和/i/音傳統聲學特征參數的SVM識別率比較，傳統聲學單個特征的平均識別率分布在50.87%~80.32%和48.37%~73.81%,而非線性聲學單特征的平均識別率分布在65.16%~87.58%和57.86%~76.39%，這說明采用非線性聲學特征參數得到的正常與病態嗓音的識別率明顯高于傳統聲學特征參數。采用模糊熵、Hurst參數兩種特征值識別/a/音和/i/音的正常與病態嗓音均取得較高的識別率，表明這兩種非線性聲學特征對病態與正常嗓音具有較好的區分度，而香農熵、第二階Renyi熵、樣本熵、計盒維數及計維截距亦能較好區分正常與病態嗓音，但不如模糊熵和Hurst參數，這與前面特征值分布分析的結果一致。此外，除了多尺度熵，其余各參數對/a/音的正常、病態嗓音的識別率均高于/i/音，這與采用傳統聲學特征參數得到的識別結果一致，因此為了達到識別病態嗓音的目的，目前國內外相關研究幾乎都采用/a/音數據，而不采用/i/音數據。
表2
/a/和/i/非線性聲學特征參數的SVM識別率
Table2.
Recognition rate of SVM on nonlinear characteristic parameters for /a/ and /i/
發/a/音時聲帶本體層自然放松其F0相對較低，而發/i/音時聲帶本體層張緊F0相對較高。因此，發/i/音時，環甲肌位置決定了其相對張緊的聲帶本體層能夠提供甲杓肌更大的代償空間，正是這種相對較高程度的代償功能的介入掩蓋了病態嗓音的特征，導致其病態嗓音的檢出率相對/a/音較低。當采用多尺度熵特征值時，可提取/i/音不同尺度的低頻信息，把隱藏在拉緊聲帶中的病理嗓音特征有效地提取出來，得到的識別率較/a/音高。多尺度熵以及/a/音、/i/音非線性特征的相對關系或能為研究聲帶代償功能的狀態提供有益的啟發,而對病態嗓音聲帶代償功能的狀態進行研究正是我們下一步的研究目標。
4 結束語
近年來，非線性特征成為描述語音信號的一個有力工具。本試驗通過比較傳統聲學特征參數和非線性聲學特征參數對/a/音及/i/音的正常與病態嗓音的SVM識別率，充分證明非線性聲學特征能更好地描述嗓音信號的病理特征，并能有效區分正常與病態嗓音。尤其是模糊熵對嗓音的識別，在/a/音取得94.67%的最高識別率，/i/音取得87.58%的最高識別率。除多尺度熵外，所有聲學特征參數對/a/音的識別率均高于/i/音，因此為了達到識別病態嗓音的目的，國內外的相關研究大多采用/a/音數據。多尺度熵對/i/音的嗓音識別率較/a/音高，表明多尺度熵能夠把隱藏在拉緊聲帶中的病理嗓音特征有效地提取出來，而發/i/音比發/a/音更易受到聲帶代償能力的影響，因此多尺度熵以及/a/音、/i/音非線性特征的相對關系或能為評價聲帶代償功能狀態的研究提供有益的啟發。
