劉 勇 袁立飛 袁麗峰

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

1 引言

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

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

2 小波分解原理

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

離散小波序列為

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

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

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

3 鼾音功率譜估計(jì)方法

3.1 Welch法

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

主要步驟如下：

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

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

3）求分段功率譜：

4）求平均功率譜：

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

3.2 多窗譜法

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

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

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

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

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

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

3.3 AR模型功率譜估計(jì)

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

對(duì)于p階AR模型，有

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

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

x(n)的功率譜為

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

3.3.1 自相關(guān)法

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

式中，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ù)都會(huì)求出，可以尋找合適的階數(shù)，自相關(guān)法相當(dāng)于加窗處理，短數(shù)據(jù)分辨率低。

3.3.2 Burg法

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

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

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

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

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

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

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

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

4 基于小波分解的鼾音譜估計(jì)

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

主要步驟如下：

1）將鼾音信號(hào)分幀處理。類(lèi)似于語(yǔ)言信號(hào)，將10ms~30ms鼾音信號(hào)看作是相對(duì)平穩(wěn)信號(hào)。

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

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

5 實(shí)驗(yàn)與分析

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

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

圖2 多窗譜法功率譜估計(jì)圖

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

圖4 Burg法功率譜估計(jì)圖

對(duì)于子帶信號(hào)S1，不同功率譜估計(jì)方法的方差如表1所示。

從圖1~4、表1可知：

表1 不同鼾音功率譜估計(jì)方法的方差

圖1 Welch法功率譜估計(jì)圖

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

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

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

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

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

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

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

圖5 不同阻塞位置子帶功率譜估計(jì)圖

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

表2 不同阻塞位置子帶信號(hào)功率譜估計(jì)的平均方差

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

6 結(jié)語(yǔ)

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