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

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

2 深圳市生物醫(yī)學工程重點實驗室，深圳市，518060

3 廣東省生物醫(yī)學信號檢測與超聲成像重點實驗室、深圳市生物醫(yī)學重點實驗室，深圳市，518060

0 引言

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

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

本研究主要從三個方面對腦電信號進行處理，分別從時域、頻域和功率譜估計等方面對腦電信號做特征分析，研究腦電信號特征變化，為處理腦電信號提供參考價值。

1 功率譜分析

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

1.1 經(jīng)典功率譜

（1）直接法

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

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

（2）Welch法

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

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

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

參數(shù)譜估計法假定所研究的信號是由輸入激勵線性系統(tǒng)的輸出模型產(chǎn)生，然后對模型的參數(shù)進行估計，再從中得到譜特性。當假設模型非常接近實際情況時，參數(shù)模型很容易響應譜中的峰值，能得到更準確的譜估計。

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

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

輸出功率譜密度為：

設b1,b2,...,bq全為0，那么：

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

1.3 功率譜估計分析

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

圖1 時域和頻域圖Fig.1 Time domain and frequency domain diagram

圖2 功率譜估計圖Fig.2 Power spectrum estimation diagram

2 時頻分析

2.1 短時傅里葉變換

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

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

式中：

短時傅里葉變換的思想是使用窗函數(shù)g(τ)來截取x(τ)，通過對截斷的局部信號進行傅里葉變換，可以獲得不同時間t信號的傅里葉變換。然后不斷移動時間t，即連續(xù)移動窗函數(shù)g(τ)的中心位置，即可連續(xù)得到不同t時刻的傅里葉變換。這些移動時間t獲得的傅里葉變換的組合就稱為短時傅里葉變換[8]。

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

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

2.2 Wigner分布

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

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

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

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

圖4 腦電信號wigner分布三維圖Fig.4 Wigner distribution of EEG

2.3 小波變換[11-12]

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

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

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

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

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

圖5 腦電信號小波去噪Fig.5 Wavelet denoising of EEG signal

3 雙譜分析[15]

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

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

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

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

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

圖6 腦電信號雙譜估計三維圖Fig.6 Three-dimensional estimation of EEG signal bispectrum

圖7 雙譜估計切片圖Fig.7 Bispectral estimation section diagram

4 應用

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