基于小波分解的鼾音信號功率譜估計方法研究*

劉 勇 袁立飛 袁麗峰

（河北省眼科醫(yī)院 邢臺 054000）

1 引言

OSAHS 病人的阻塞位置需要經(jīng)睡眠呼吸監(jiān)測結(jié)合纖維喉鏡檢查、影像學(xué)檢查才能確定，成本高且需要在醫(yī)院才能進(jìn)行，打鼾是OSAHS 病人的典型癥狀，研究不同阻塞位置鼾音的功率譜估計［1］特點，為診斷OSAHS 阻塞位置提供一種簡單方便的檢測方法。功率譜估計方法分為傳統(tǒng)譜估計法［2］和現(xiàn)代譜估計法［3］。傳統(tǒng)譜估計法原理簡單，便于實現(xiàn)，在實際中應(yīng)用較廣，常用的方法有Welch法［4~6］、多窗譜法［7~8］等?，F(xiàn)代譜估計是先將實驗信號建立參數(shù)模型，再依據(jù)模型估計信號的功率譜。線性模型在信號處理中應(yīng)用最多，具體分為AR 模型［9］、MA模型、ARMA模型。常用的方法有基于AR模型的自相關(guān)法［10］和Burg 法［11］。分析鼾音和阻塞位置的不同功率譜估計特點，為診斷OSAHS 阻塞位置提供幫助。

針對上述問題，提出了基于小波分解的鼾音功率譜估計方法，首先將鼾音信號進(jìn)行小波分解［12~14］，子帶信號代表鼾音信號的不同時頻域信息，對各子帶信號分別采用Welch法、多窗譜法、基于AR 模型的自相關(guān)法和Burg 法。結(jié)果表明，基于小波分解的Welch 法得到的鼾音功率譜估計結(jié)果分辨率和方差可以滿足需要，同時初步分析了上部阻塞和下部阻塞功率譜估計特征。

2 小波分解原理

設(shè)ψ(t)∈L2(R) ，傅里葉變換表示ψ?(ω) 。當(dāng)時，ψ(t)表示母小波。小波序列是指ψ(t)伸縮、平移所得的序列。

離散小波序列為

對于二進(jìn)離散小波變換，信號{X(t)} 通過小波函數(shù)ψj，k(t)，尺度函數(shù)φj，k(t)進(jìn)行J層分解：

dj，k為尺度j上的細(xì)節(jié)信息，稱為小波系數(shù)；aJ，k為尺度J上的逼近信息，為近似系數(shù)。

選擇不同的尺度，其時頻分辨率不同，小波分解存在逆變換，將小波系數(shù)和近似系數(shù)重構(gòu)［15］，可得到信號中的不同頻率成分。

3 鼾音功率譜估計方法

3.1 Welch法

Welch 法是將整個信號分段處理，相鄰的每段進(jìn)行重疊，然后對數(shù)據(jù)加窗進(jìn)行譜估計，是周期圖法的改進(jìn)。

主要步驟如下：

1）將數(shù)據(jù)xN(n) 分成L段，每段數(shù)據(jù)長度為M，即N=LM。對xN(n)分段時，假設(shè)50%重疊，則第i段數(shù)據(jù)表示：

2）對每段數(shù)據(jù)加窗w(n)傅里葉變換，表示為

3）求分段功率譜：

4）求平均功率譜：

通過分段重疊和加窗的方法，譜估計結(jié)果是逐步一致估計，并且可以減少方差。窗函數(shù)是Welch法的關(guān)鍵因素，為提高分辨率和減小譜泄漏，理想是選擇旁瓣幅值小，主瓣寬度小的窗函數(shù)，但是實際中兩者不能同時小，需要根據(jù)信號特點和實際情況，選擇合適的窗函數(shù)。

3.2 多窗譜法

多窗譜法是對信號采用多個互相正交的窗函數(shù)，然后進(jìn)行傅里葉變換，最后得到信號的譜估計。

設(shè)時間序列為x(1)，…，x(N)，S(f)為真實譜。存在K個互相正交的數(shù)據(jù)窗，{ }ht，k:t1，…，N代表第k個窗，K個直接譜的均值為該序列的多窗譜Smt(f)。表示為

其中Smt(f)是對數(shù)據(jù)加第k個窗的直接譜。表示為

多正弦窗，即正交窗和最小偏差窗的集合。對于正弦窗：

多正弦窗主瓣窄，分析式簡單，就有較低的偏差。

為防止譜泄漏，采用一組相互正交的數(shù)據(jù)窗代替單個數(shù)據(jù)窗，這樣可以利用信號兩端的數(shù)據(jù)，進(jìn)而減小方差，提升譜估計性能。

3.3 AR模型功率譜估計

原理是現(xiàn)在的輸入和過去p個輸出的加權(quán)之和表示該模型現(xiàn)在的輸出。

對于p階AR模型，有

其中w(n)為白噪聲；a1，a2，…ap為p階AR 模型的參數(shù)。

該AR模型系統(tǒng)的傳遞參數(shù)為

x(n)的功率譜為

AR模型是全極點的模型。平穩(wěn)信號可以通過AR模型表示，AR模型計算簡單，是最常用的模型，在語音信號處理等領(lǐng)域應(yīng)用廣泛。

3.3.1 自相關(guān)法

自相關(guān)法先估計自相關(guān)序列，然后通過萊文森遞歸公式算出AR系數(shù)，計算出功率譜。

式中，r(1)，r(2)，…，r(n+1) 為相 關(guān)系數(shù)，a(2)，…，a(n+1)為自回歸系數(shù)。

自相關(guān)算法功率譜公式為

式中，e(f)為負(fù)數(shù)正弦曲線。

萊文森遞推是從低到高直到p階，p階之前的參數(shù)都會求出，可以尋找合適的階數(shù)，自相關(guān)法相當(dāng)于加窗處理，短數(shù)據(jù)分辨率低。

3.3.2 Burg法

Burg 法是選擇反射系數(shù)使得前向和后向總的均方誤差的和最小，然后通過萊文森遞推算法算出模型系數(shù)。

數(shù)據(jù)x(1)，…，x(N)，前向預(yù)測誤差表示：

后向預(yù)測誤差表示：

各階前向和后向預(yù)測誤差為

前向和后向預(yù)測誤差的平均功率為

為使Pk最小，令?Pk?λk=0 ，得到反射系數(shù)為

利用萊文森遞歸公式求解AR模型的參數(shù)為

Burg 法中階次的選擇影響譜估計的質(zhì)量，階次越高，頻率分辨率越高，但是容易產(chǎn)生譜線開裂和譜峰偏差的問題，同時對于低噪正弦信號敏感。

4 基于小波分解的鼾音譜估計

小波分解可以把采集到鼾音信號投影到小波基函數(shù)的空間中，將指定尺度上的小波系數(shù)和近似系數(shù)重構(gòu)鼾音信號，可獲得對應(yīng)的頻率子帶信號。將鼾音信號小波分解，提取相應(yīng)的子帶實現(xiàn)鼾音的時頻域精細(xì)化分析，對不同子帶信號進(jìn)行功率譜估計。

主要步驟如下：

1）將鼾音信號分幀處理。類似于語言信號，將10ms~30ms鼾音信號看作是相對平穩(wěn)信號。

2）對每幀的鼾音信號進(jìn)行小波4 層分解，選取3、4 層的小波系數(shù)重構(gòu)，分離出相應(yīng)的頻率成分，子帶信號表示為S1、S2 、S3，對應(yīng)的頻率為0~1000Hz，1000Hz~2000Hz，2000Hz~4000Hz。

3）對子帶信號采用Welch 法、多窗譜法、自相關(guān)法、Burg法進(jìn)行功率譜估計得到各子帶鼾音信號的功率譜估計。

5 實驗與分析

采集河北省眼科醫(yī)院150 名OSAHS 病人的鼾音信息，電容麥克風(fēng)位于病人口部上方15cm 處，采集的音頻為WAV 格式，采樣頻率為32000Hz，精度為16Bit。

抽取一次鼾音信號作為數(shù)據(jù)樣本進(jìn)行仿真實驗。對鼾音信號分幀，幀長為1024，根據(jù)鼾音信號的特點，采用dB4 小波基，對鼾音信號進(jìn)行小波4層分解，對第1 個子帶鼾音信號S1 分別采用Welch法、多窗譜法、自相關(guān)法、Burg 法進(jìn)行功率譜估計。Welch法中采用Blackman窗和Hamming窗，窗長為1024，段重疊為512，F(xiàn)FT長度為1024。多窗譜法中采用多正弦窗，時間帶寬積為2 和4。自相關(guān)法中階數(shù)為80 和300。Burg 法中階數(shù)為80 和300。子帶信號S1 的不同功率譜估計結(jié)果分別圖1、圖2、圖3、圖4所示。

圖2 多窗譜法功率譜估計圖

圖3 自相關(guān)法功率譜估計圖

圖4 Burg法功率譜估計圖

對于子帶信號S1，不同功率譜估計方法的方差如表1所示。

從圖1~4、表1可知：

表1 不同鼾音功率譜估計方法的方差

圖1 Welch法功率譜估計圖

1）Welch 法 中Hamming 窗 的 頻 率 分 辨 率 比Blackman窗高，而且方差更小。窗函數(shù)的選擇至關(guān)重要，由于Hamming 窗的主瓣寬度更小，因此具有更好的分辨率性能。

2）多窗譜法中時間帶寬積等于2比等于4時分辨率高，方差增大。隨著時間帶寬積增大，窗的數(shù)目會增多，導(dǎo)致頻率分辨率降低，時間分辨率增加，同時譜泄露會增加。

3）自相關(guān)法和Burg 法性能非常接近，需要提高階次才能提高分辨率，階數(shù)越高有更多的譜細(xì)節(jié)，譜峰信息比較尖銳，其余地方比較平滑，但是高階次容易導(dǎo)致譜峰偏移。

4）從整體來看，分辨率方面，采用Hamming 窗Welch 法分辨率最高，其次是采用Blackman 窗Welch法，然后是階數(shù)為300的自相關(guān)法和Burg法，最后是時間帶寬積為2 的多窗譜法；方差方面，在分辨率相對高的情況下，方差由大到小依次為采用Blackman 窗Welch 法、采用Hamming 窗Welch 法、階數(shù)為300 的自相關(guān)法、階數(shù)為300 的Burg 法、時間帶寬積為2的多窗譜法。

5）綜上所述，采用hamming 窗Welch 法對鼾音的功率譜估計性能最好，分辨率和方差可以滿足需要。

鼾音是睡眠時呼吸道產(chǎn)生的聲音，沒有標(biāo)準(zhǔn)的數(shù)字模型，只能借助語音信號［16］的功率譜估計方法。在語音處理領(lǐng)域Welch 法是最常用的功率譜估計方法，從結(jié)果上看依然適用于鼾音信號功率譜估計；多窗譜法適合短數(shù)據(jù)、頻譜波動大的信號，如地震監(jiān)測等，在鼾音領(lǐng)域不如Welch 法；AR 模型是全極點模型，鼾音信號類似于語音信號，有極點和零點。 基于AR 模型的自相關(guān)法和Burg 法結(jié)果是傳統(tǒng)譜估計的包絡(luò)線，階數(shù)越高有更多得譜細(xì)節(jié)，更重視譜峰信息，其余地方比較平滑。

臨床上不同的阻塞位置采用不同的手術(shù)方法，分析上部阻塞和下部阻塞的功率譜特性來判斷阻塞位置。根據(jù)上述結(jié)論小波分解的基礎(chǔ)上選擇Welch法，對3個子帶信號進(jìn)行功率譜估計，上部阻塞和下部阻塞功率譜估計結(jié)果如圖5所示。

圖5 不同阻塞位置子帶功率譜估計圖

150例OSAHS病人的鼾音信號采用本文方法，上部阻塞和下部阻塞子帶信號功率譜估計的平均方差如表2所示。

表2 不同阻塞位置子帶信號功率譜估計的平均方差

綜合圖5、表2 可以看出上部阻塞和下部阻塞鼾音信號功率譜估計在方差方面，0~1000Hz 子帶和1000Hz~2000Hz 子帶差別不明顯，2000Hz~4000Hz 子帶的方差差別顯著，下部阻塞信號譜估計方差遠(yuǎn)大于上部阻塞的鼾音信號，可以作為判斷阻塞位置的特征值。

6 結(jié)語

提出了基于小波分解的鼾音譜估計方法，將鼾音信號進(jìn)行小波分解，再對提取的三個子帶信號分別采用Welch 法、多窗譜法、自相關(guān)法、Burg 法進(jìn)行功率譜估計。研究表明，基于小波分解Welch 法適合于鼾音領(lǐng)域。利用本文方法分析150 例OSAHS阻塞位置，實驗發(fā)現(xiàn)上部阻塞信號在2000Hz~4000Hz 子帶的平均方差遠(yuǎn)小于下部阻塞信號，可以將該特征信息作為判斷阻塞位置的標(biāo)準(zhǔn)，提供了一種簡單方便的檢測方法。