腦電信號(hào)分析方法及其應(yīng)用

周晶晶，葉繼倫,，張旭,，李晨洋，檀雪

1 深圳大學(xué) 醫(yī)學(xué)部 生物醫(yī)學(xué)工程系，深圳市，518060

2 深圳市生物醫(yī)學(xué)工程重點(diǎn)實(shí)驗(yàn)室，深圳市，518060

3 廣東省生物醫(yī)學(xué)信號(hào)檢測(cè)與超聲成像重點(diǎn)實(shí)驗(yàn)室、深圳市生物醫(yī)學(xué)重點(diǎn)實(shí)驗(yàn)室，深圳市，518060

0 引言

EEG信號(hào)是研究大腦科學(xué)活動(dòng)的重要手段之一，可使用電生理指標(biāo)來(lái)記錄反映大腦頭皮神經(jīng)元興奮性和抑制性突觸后電位的潛在活動(dòng)[1]。用ECG來(lái)分析大腦的活動(dòng)，具有無(wú)創(chuàng)、實(shí)時(shí)、連續(xù)、時(shí)間分辨率高以及認(rèn)知意識(shí)相關(guān)性強(qiáng)等特性，是目前臨床醫(yī)學(xué)研究中不可或缺的手段，也是實(shí)驗(yàn)室重點(diǎn)研究的生理信號(hào)。利用腦電信號(hào)反映麻醉深度能夠指導(dǎo)臨床手術(shù)，研究腦電信號(hào)在疲勞狀態(tài)下的性質(zhì)，能夠判斷司機(jī)疲勞駕駛程度等，因此研究腦電信號(hào)具有很重要的作用。

EEG信號(hào)的頻率分布主要集中在 0.5～30 Hz，EEG信號(hào)幅值大多在5～200 μV之間[2-3]，是典型的極其微弱的非平穩(wěn)隨機(jī)信號(hào)。EEG信號(hào)的臨床應(yīng)用通常根據(jù)頻率段將EEG分成幾個(gè)不同的信號(hào)段，依次稱為δ波（1～3 Hz）、θ波（4～7 Hz）、α波（8～13 Hz）、β波（14～30 Hz）[4]。δ波在睡眠、缺氧、疲勞、麻醉等狀態(tài)下會(huì)出現(xiàn)，θ波會(huì)出現(xiàn)在疲憊時(shí)，α波出現(xiàn)在清醒和閉上眼睛時(shí)，β波出現(xiàn)在精神波動(dòng)和情緒緊張或興奮時(shí)。不同大腦狀態(tài)，各頻率段所占能量不同，為研究腦電活動(dòng)提供了基礎(chǔ)。

本研究主要從三個(gè)方面對(duì)腦電信號(hào)進(jìn)行處理，分別從時(shí)域、頻域和功率譜估計(jì)等方面對(duì)腦電信號(hào)做特征分析，研究腦電信號(hào)特征變化，為處理腦電信號(hào)提供參考價(jià)值。

1 功率譜分析

功率譜分析方法是一種非常重要的頻譜分析方法，功率譜分析主要探討信號(hào)能量隨頻率的變化情況，即在頻域的信號(hào)能量分布狀況。功率譜估計(jì)主要分為兩類：經(jīng)典功率譜（如直接法、自相關(guān)法）和參數(shù)模型（AR模型）功率譜[5]。經(jīng)典功率譜估計(jì)方法的性能較差，分辨率較低。性能差的主要原因是不能實(shí)現(xiàn)功率譜密度原始定義中平均值和極限的運(yùn)算，分辨率低的原因是周期圖法假設(shè)數(shù)據(jù)窗外的數(shù)據(jù)全是零，對(duì)于自相關(guān)方法是假設(shè)在延遲窗以外的自相關(guān)函數(shù)全是零。參數(shù)模型法是現(xiàn)代譜估計(jì)的主要內(nèi)容，具有較高的頻率分辨率和性能。

1.1 經(jīng)典功率譜

（1）直接法

直接法又叫做周期圖法，它將隨機(jī)信號(hào)x(n)的N點(diǎn)觀察數(shù)據(jù)xN(n)視為一個(gè)能量有限信號(hào)，直接對(duì)xN(n)進(jìn)行傅里葉變換得xN(ejw)。然后，取其幅值大小的平方，并除以N，作為對(duì)x(n)真實(shí)的功率譜P(ejw)的估計(jì)。由周期圖法估計(jì)出來(lái)的功率譜PPER(ejw)表示:

直接法截短無(wú)限長(zhǎng)平穩(wěn)信號(hào)，相當(dāng)于對(duì)信號(hào)添加矩形窗，使窗函數(shù)和功率譜卷積。該卷積導(dǎo)致頻譜泄露，并且當(dāng)數(shù)據(jù)長(zhǎng)度N太大時(shí)，譜曲線波動(dòng)加劇，并且如果N太小，譜的分辨率會(huì)較差。

（2）Welch法

Welch法[6]是對(duì)周期圖方法的改進(jìn)，利用平均法改進(jìn)其方差特性，得到加權(quán)交疊平均法，即Welch平均周期圖法。它的指導(dǎo)思想是將xN(n)數(shù)據(jù)長(zhǎng)度N分成L段，并分別計(jì)算每一段的功率譜，然后對(duì)其平均，在分割時(shí)允許每個(gè)數(shù)據(jù)存在部分重疊。每段的窗口可以是矩形窗口，可以是漢明窗也可以是漢寧窗。根據(jù)Welch法，每段的功率譜記為即

式中M為每段數(shù)據(jù)，也是窗長(zhǎng)度，是歸一化因子，它用于保證得到的頻譜是漸進(jìn)式無(wú)偏估計(jì)。因此，Welch平均法得到的平均功率譜為：

1.2 參數(shù)模型功率譜

參數(shù)譜估計(jì)法假定所研究的信號(hào)是由輸入激勵(lì)線性系統(tǒng)的輸出模型產(chǎn)生，然后對(duì)模型的參數(shù)進(jìn)行估計(jì)，再?gòu)闹械玫阶V特性。當(dāng)假設(shè)模型非常接近實(shí)際情況時(shí)，參數(shù)模型很容易響應(yīng)譜中的峰值，能得到更準(zhǔn)確的譜估計(jì)。

無(wú)論x(n)是確定性信號(hào)還是隨機(jī)信號(hào)，參數(shù)模型的線性系統(tǒng)都有如下輸入輸出關(guān)系：

u(n)是方差為σ2的白噪聲序列，對(duì)兩邊各取Z變換，并假定b0=1，可以得到傳遞函數(shù)：

輸出功率譜密度為：

設(shè)b1,b2,...,bq全為0，那么：

我們把（8）（9）（10）三式稱作自回歸模型，即AR模型[7]，它是一個(gè)全極點(diǎn)模型?！白曰貧w”的含義是：模型的輸出是當(dāng)前輸入和過(guò)去p個(gè)輸出的加權(quán)和。由于AR模型是有理分式，因此估計(jì)出的譜比經(jīng)典法的譜光滑無(wú)毛刺并且具有高分辨率的特征。

1.3 功率譜估計(jì)分析

該研究采用直接法和Welch法兩種經(jīng)典模型的功率譜估計(jì)和AR參數(shù)模型的功率譜估計(jì)對(duì)腦電信號(hào)進(jìn)行功率譜分析，圖1是原始時(shí)域EEG信號(hào)及其幅值譜，圖2是三種功率譜分析圖，從圖2中可以明顯看出，直接法效果較差，Welch法和AR模型估計(jì)的功率譜效果較好，從功率譜分布圖可以得出該信號(hào)能量主要集中在δ段，可能是處于麻醉深度或疲勞狀態(tài)，對(duì)臨床上判斷麻醉深度具有參考作用。

圖1 時(shí)域和頻域圖Fig.1 Time domain and frequency domain diagram

圖2 功率譜估計(jì)圖Fig.2 Power spectrum estimation diagram

2 時(shí)頻分析

2.1 短時(shí)傅里葉變換

由于EEG是時(shí)變非平穩(wěn)信號(hào)，因此信號(hào)的頻率隨時(shí)間改變而改變，傳統(tǒng)的傅里葉變換缺少時(shí)域定位的功能，為了解EEG信號(hào)時(shí)頻信息，我們采用短時(shí)傅里葉變換對(duì)腦電信號(hào)進(jìn)行處理。

給定一信號(hào)x(t)∈L2(R)：

式中：

短時(shí)傅里葉變換的思想是使用窗函數(shù)g(τ)來(lái)截取x(τ)，通過(guò)對(duì)截?cái)嗟木植啃盘?hào)進(jìn)行傅里葉變換，可以獲得不同時(shí)間t信號(hào)的傅里葉變換。然后不斷移動(dòng)時(shí)間t，即連續(xù)移動(dòng)窗函數(shù)g(τ)的中心位置，即可連續(xù)得到不同t時(shí)刻的傅里葉變換。這些移動(dòng)時(shí)間t獲得的傅里葉變換的組合就稱為短時(shí)傅里葉變換[8]。

如圖3所示，是Matlab下處理腦電信號(hào)的短時(shí)傅里葉變換圖。從圖中可看出信號(hào)頻率隨時(shí)間變換的關(guān)系，然而短時(shí)傅里葉變換也具有局限性，不能同時(shí)優(yōu)化頻率與時(shí)間分辨率，且時(shí)頻窗的面積不能小于1/2。

圖3 短時(shí)傅里葉變換Fig.3 Short-time Fourier transform

2.2 Wigner分布

Wigner分布又稱Wigner-Ville分布，簡(jiǎn)稱為WVD。令信號(hào)x(t),y(t)的傅里葉變換分別是X(jΩ),Y(jΩ),那么x(t),y(t)的聯(lián)合Wigner分布定義為：

信號(hào)x(t)的自Wigner分布定義為：

Wigner分布反映信號(hào)能量隨時(shí)間和頻率分布的關(guān)系，具有一系列良好的時(shí)頻性質(zhì)[9]，如對(duì)稱、時(shí)移、組合和復(fù)共扼等，不會(huì)損失信號(hào)的幅值與相位信息，對(duì)瞬時(shí)頻率和群延時(shí)有清晰的概念，所以在時(shí)頻分析中具有很重要的地位。但由于Wigner分布是雙線性變換的形式，因此兩個(gè)信號(hào)和的分布將存在交叉項(xiàng)干擾[10]，它會(huì)影響信號(hào)時(shí)頻行為的識(shí)別，這也是二次時(shí)頻分布的必然結(jié)果。

圖4是Matlab下EEG的Wigner三維分布圖，可以看出它的時(shí)頻分辨率優(yōu)于短時(shí)傅里葉變換。

圖4 腦電信號(hào)wigner分布三維圖Fig.4 Wigner distribution of EEG

2.3 小波變換[11-12]

小波變換可以分析局部信號(hào)的時(shí)間和頻率[13]，通過(guò)平移和伸縮對(duì)信號(hào)（函數(shù)）逐步進(jìn)行多尺度細(xì)化，最終可以使分析的信號(hào)在高頻處具有高時(shí)間分辨率，低頻處具有高頻率分辨率，因此小波變換能專注于信號(hào)的任意細(xì)節(jié)，解決了傅里葉變換的難題，并且可以適應(yīng)不同頻率范圍內(nèi)不同分辨率的信號(hào)分析的基本要求。

小波變換的定義是給定基本函數(shù)ψ(t)：

式中：a,b均為常數(shù)，且a＞0。ψa,b(t)是對(duì)基本函數(shù)做ψ(t)平移和伸縮之后得到的函數(shù)。如果a,b連續(xù)變化，可得ψa,b(t)函數(shù)的集合，給定平方可積的信號(hào)x(t)，即x(t)∈L2(R)，則x(t)的小波變換定義為：

式中a,b和t均是連續(xù)變量。所以該式也稱為連續(xù)小波變換（CWT），信號(hào)x(t)的小波變換是a和b的函數(shù)，b是時(shí)移，a是尺度因子。ψ(t)又叫基本小波，或者母小波，ψa, b(t)是母小波經(jīng)平移和伸縮后產(chǎn)生的函數(shù)集合，稱為小波基函數(shù)。因此WT又可解釋為信號(hào)x(t)和一組小波基的內(nèi)積。

小波變換能夠在時(shí)域和頻域中描述信號(hào)局部的特征。小波變換可以將信號(hào)的頻譜分解成不同的頻率段，然后把信號(hào)的能量集中到某些頻帶的幾個(gè)系數(shù)上，通過(guò)把其他頻帶上分解的系數(shù)置零或是給小一點(diǎn)的權(quán)重，能夠較好抑制噪聲[14]。該研究使用db4小波對(duì)EEG信號(hào)去噪，可以取得良好的效果。如圖5所示，小波去噪后的信號(hào)比原始信號(hào)更平滑，邊界更清晰，而且可以去掉很多毛刺。

圖5 腦電信號(hào)小波去噪Fig.5 Wavelet denoising of EEG signal

3 雙譜分析[15]

因?yàn)槟X電信號(hào)中存在大量非線性耦合現(xiàn)象，所以研究腦電信號(hào)的非高斯特征、非線性特征也很重要[16]。雖然功率譜分析可以反映信號(hào)的二階信息，但是會(huì)丟失包括相位信息在內(nèi)的高階信息[17]，因此研究腦電信號(hào)的非線性特征最常用的方法是對(duì)信號(hào)進(jìn)行高階累計(jì)量的分析。

雙譜密度函數(shù)定義為：

雙譜分析具有很多性質(zhì)，如雙周期、對(duì)稱性、復(fù)函數(shù)，可通過(guò)振幅與相位來(lái)表示等。

根據(jù)觀測(cè)數(shù)據(jù)長(zhǎng)度的限制，實(shí)際信號(hào)的雙譜估計(jì)通常分為兩大類[18]，第一類是參數(shù)模型估計(jì)法，第二類為非參數(shù)估計(jì)法，其中非參數(shù)估計(jì)法包括直接估計(jì)法和間接估計(jì)法。直接估計(jì)法使用的是離散傅里葉變化的三階周期法，間接估計(jì)法通過(guò)對(duì)三階累積量加窗并進(jìn)行二維傅里葉變換得到。該研究采用雙譜直接估計(jì)法。

圖6為在Matlab環(huán)境下對(duì)腦電信號(hào)做雙譜估計(jì)的三維圖，雙譜圖中包含信號(hào)的能量、幅度信息、信號(hào)的相位，可以作為檢測(cè)EEG的非高斯和非線性特性的重要工具[19]。雙譜估計(jì)三維圖不能準(zhǔn)確描述固定頻率與其他頻率的耦合，而二維切片可以清楚地得到該信號(hào)和其他頻率的非線性耦合。如圖7所示，是該EEG雙譜估計(jì)三維圖的平面切圖，能夠反映雙譜分析的一部分信息，提取雙譜特征值將有益于分類器基于特征值的分類[20]。

圖6 腦電信號(hào)雙譜估計(jì)三維圖Fig.6 Three-dimensional estimation of EEG signal bispectrum

圖7 雙譜估計(jì)切片圖Fig.7 Bispectral estimation section diagram

4 應(yīng)用

該研究建立了初步的腦電測(cè)量系統(tǒng)，根據(jù)以上的信號(hào)處理介紹，選擇了相關(guān)的方法進(jìn)行噪聲及特征參數(shù)的分析，其中功率譜計(jì)算主要用于分析腦電信號(hào)各頻段的能量所占的比例，對(duì)于麻醉和疲勞等狀態(tài)下EEG不同頻段分期具有良好的參考價(jià)值，可用于麻醉深度的檢測(cè)和疲勞程度的劃分，對(duì)于指導(dǎo)腦電信號(hào)的應(yīng)用具有重要的意義。時(shí)頻分析多用于光滑原始信號(hào)，獲得較高質(zhì)量的信號(hào)，從而進(jìn)行特征提取等。雙譜分析可顯示傳統(tǒng)腦電圖無(wú)法顯示的信息，以及通過(guò)腦電信號(hào)三階能量在雙頻域中各頻段的分布，研究不同腦功能狀態(tài)下的腦電信號(hào)，為我們了解大腦功能提供了一種新的分析方法。