本發(fā)明涉及一種實用化在線腦電偽跡剔除方法，屬于生物醫(yī)學信息處理技術領域。

背景技術：

腦電(electroencephalography,eeg)蘊含了大量的心理、生理和病理信息，目前被廣泛運用于腦部疾病診斷和腦-機接口(brain-computerinterface,bci)等諸多研究領域中。正常腦電信號主要分布在δ頻帶(0-4hz)，θ頻帶(4-8hz)，α頻帶(8-13hz)，β頻帶(13-20hz)，γ頻帶(30-80hz)五個頻帶范圍內(nèi)，而eeg的幅度十分微弱，易受眼電(electrooculogram,eog)、肌電(electromyogram,emg)、50hz工頻以及線性漂移等各種偽跡的干擾，這些偽跡存在嚴重影響了bci的性能。

以往的一些偽跡剔除方法中，偽跡拒絕無可避免使腦電受到eog、emg等偽跡的污染；線性濾波無法適用于偽跡與正常腦電信號的頻帶發(fā)生重疊的情況；線性回歸的方法依賴于一個好的偽跡參考信號，通常我們很難找到一個適合于emg偽跡和非生理偽跡的參考信號。一些研究者直接使用帶通濾波器或設定閾值的方法來剔除原始腦電中高頻肌電、噪聲和眼電信號，這種方法雖然簡單，但是往往在剔除偽跡的同時也剔除了部分有用的腦電信號，導致預處理效果不理想。盲源分離(blindsignalseparation,bbs)技術是一種最有前景的腦電偽跡剔除方法，獨立成分分析(independentcomponentanalysis,ica)就是bbs方法中的一種。

ica是信號處理研究中常用的一種統(tǒng)計方法，它旨在尋找能最大化數(shù)據(jù)中各個分量獨立性的線性投影。它假設從電極陣列采集到的n維腦電信號x＝[x1,x2,…,xn]^t，該信號是由m個未知且統(tǒng)計獨立的信號源s＝[s1,s2,…,sm]^t的線性組合。ica算法的目標是尋找分離矩陣w，使得：s＝wx。一些改進的ica算法已經(jīng)被實現(xiàn)和免費使用，其中，fastica相比于其它改進的ica算法而言主要有以下優(yōu)點：a、fastica算法收斂速度快，能夠提高算法的實時性；b、該算法使用方便，數(shù)據(jù)選取時不需要設置步長參數(shù)；c、它具有類似于神經(jīng)網(wǎng)絡算法并行、分布的優(yōu)點，計算復雜度低，占用內(nèi)存空間小，適合于在線腦電信號的分析。通常，ica要求觀測信號的數(shù)量要不少于源信號的數(shù)量，當觀測信號的數(shù)量少于源信號數(shù)量時則會出現(xiàn)超完備ica現(xiàn)象。

目前，多通道腦電信號中偽跡的剔除算法已經(jīng)較為成熟。而在單/少通道腦電信號的偽跡剔除中，由于腦電采集通道數(shù)量較少，所包含的有效信息較少，且缺乏偽跡參考信號，目前尚無十分有效的偽跡剔除方法。此外，多導聯(lián)腦電信號設備由于價格昂貴，體型巨大，操作繁瑣，往往不能滿足便捷采集、處理腦電的要求。近年來，為了解決ica算法用于單/少通道腦電信號分析中存在的超完備問題，通道變換技術隨著便攜式腦電設備的研究而興起，bogdanmijovic針對單通道腦電信號提出經(jīng)驗模態(tài)分解(empiricalmodedecomposition,emd)和ica相結合的腦電偽跡剔除方法，即emd-ica方法，而李明愛針對少通道腦電信號提出將離散小波變換(discretewavelettransform,dwt)和獨立成分分析相結合的方法，即dwt-ica方法，該方法主要針對腦電信號中的眼電偽跡，而且沒有在線運用，本發(fā)明主要是在該方法的基礎上進行了改進，以提高偽跡成分自動識別和剔除能力，并進行在線運用。

小波變換(wavelettransform,wt)通過將信號分解成多尺度的小波函數(shù)來獲取信號在時域和頻率的局部特征，而dwt的計算速度快，廣泛的應用于腦電這種非線性、非平穩(wěn)信號的分析。傳統(tǒng)的ica算法直接對含高斯噪聲的觀測信號進行分解時容易產(chǎn)生較大的誤差，還增大了計算量，大幅降低了分離效果，為此，本發(fā)明在ica分離之前引入了dwt，dwt在將少通道腦電轉化為適合ica分析的多通道腦電的同時增強了ica的分離效果。雖然ica能夠成功的將多通道原始腦電信號中的真實腦電和偽跡成分分離，但是在這種方法中偽跡分量的識別往往依賴于專業(yè)人員的視覺檢查，費時且?guī)в泻艽蟮闹饔^性，單純的ica算法還遠遠達不到在線自動高效識別偽跡成分的要求。

在偽跡的在線自動識別方面，聚類技術是一種高效的方法，它通?；诿總€獨立分量的一些特定特征，用于自動從各個獨立成分中分離出偽跡相關的成分。一些研究者基于兩個獨立成分之間的相似度運用了k均值算法和模糊c均值聚類算法進行聚類，自動識別并剔除偽跡。而這兩種聚類算法都是一種迭代的方法，都需要事先給定聚類的目標數(shù)目來確定迭代的次數(shù)，實際聚類過程中聚類的目標數(shù)目往往是未知的。層次聚類(hierarchicalclustering,hc)算法相比于上述兩種聚類算法而言主要有兩大優(yōu)點：首先，層次聚類過程中的樹狀圖不僅包含各個類的組成信息，而且在每個類中還以節(jié)點高度的形成提供了每個類中元素的親密信息；其次，整個聚類過程不需要事先設定聚類的數(shù)目。

基于運動想象(motorimagery,mi)的bci相對于基于穩(wěn)態(tài)視覺誘發(fā)電位(stady-statevisualevokedpotential,ssvep)和p300電位的bci而言，它不依賴于屏幕的同步視覺刺激，可以實現(xiàn)異步bci，而且被試想象左右手運動意圖可以通過分析被試初級運動區(qū)c3和c4導聯(lián)的腦電信號加以識別，因而基于該類范式的bci更容易應用于便攜式bci設備中。近年，新興的便攜式bci、便攜式麻醉深度監(jiān)測系統(tǒng)、便攜式睡眠監(jiān)測、便攜式情緒狀態(tài)識別等應用領域通常針對在線單/少通道的腦電信號進行分析。因此，研究實用的在線少通道運動想象腦電偽跡剔除算法對于便攜式腦電設備的發(fā)展具有重要的意義。

技術實現(xiàn)要素：

本發(fā)明提供了一種結合dwt、ica和hc的實用化在線腦電偽跡剔除方法，即dwicac方法，以用于解決目前運動想象腦電預處理算法中，多種常規(guī)腦電偽跡無法在單通道的情況下，在線自動識別并剔除的問題。

本發(fā)明的技術方案是：一種實用化在線腦電偽跡剔除方法，所述方法的具體步驟如下：

a、對各通道的原始腦電信號進行降采樣、工頻陷波和線性漂移校正；

b、對校正之后各通道的腦電信號進行多尺度小波分解和小波系數(shù)單支重構；

c、分別對各通道的小波系數(shù)單支重構信號進行ica分解，計算各個分量的時域特征、頻譜特征以及序間相似度特征；

d、對各通道的獨立分量分別進行層次聚類，自動識別并剔除偽跡分量，并對剩余的獨立分量進行fastica逆變換重構和小波系數(shù)重構，得到干凈的腦電信號。

所述步驟a中的降采樣、工頻陷波和線性漂移校正具體方法如下：

a1、利用matlab信號處理工具箱中的downsample函數(shù)分別對c3、c4通道實時的原始腦電信號進行250hz降采樣；

a2、利用matlab信號處理工具箱中的dlsim函數(shù)對c3、c4通道降采樣之后的腦電數(shù)據(jù)進行50hz工頻陷波處理；

a3、利用matlab信號處理工具箱中的detrend函數(shù)對陷波處理之后的信號進行線性漂移校正；

經(jīng)過上述處理得到校正之后的腦電信號x(t)＝[x1(t),x2(t)，…，xm(t)]^t，其中：x(t)是一組m維隨機向量，m為腦電采集通道數(shù)，t為轉置運算，t＝1,2,…,n，n為腦電信號x(t)的采樣點個數(shù)。

所述步驟b中，多尺度小波分解和小波系數(shù)單支重構的具體步驟如下：

b1、選擇“db4”小波基，利用matlab軟件離散小波變換工具箱中的wavedec函數(shù)對某一通道的信號x1(t)進行7層離散小波分解，得到信號x1(t)的逼近系數(shù)分量ca7和細節(jié)系數(shù)分量cd7、cd6、cd5、cd4、cd3、cd2、cd1；

b2、采用離散小波變換工具箱中的wrcoef函數(shù)分別對步驟b1中得到的逼近系數(shù)分量和細節(jié)系數(shù)分量進行單支重構，得到對應于b1中ca7,cd7，cd6,…,cd1的小波重構信號r(t)＝[r1(t),r2(t),…,r8(t)]^t。

所述步驟c中，ica分解和各個分量的時域特征、頻譜特征和序間相似度計算方法如下：

c1、將步驟b2中的r(t)作為fastica的輸入，利用matlab軟件下fastica工具箱中的fastica函數(shù)對r(t)進行ica分解，得到分離矩陣w和8個獨立分量(s1,s2,…,s8)；

c2、對步驟c1得到的8個獨立分量計算峰度值作為各個獨立分量的時域特征，式中，si表示第i個獨立分量，i∈[1,8]，e()為統(tǒng)計期望函數(shù)，式中，l為獨立分量si的數(shù)據(jù)量，為獨立分量si中第l個數(shù)據(jù)點的m次方，和為獨立分量i的4階矩和2階矩；

c3、計算步驟c1得到的8個獨立分量在5個頻帶范圍內(nèi)的平均頻帶功率：fspectral＝[f(δ)f(θ)f(α)f(β)f(γ)]，式中，f()是平均頻帶功率；

c4、對步驟c1中8個獨立分量s1,s2,…,s8設置滑動窗時間長度為l的通道獨立分量數(shù)據(jù)作為一個子序列，總共得到nep個子序列，一個時序中通道的8獨立分量最后就分成了nep個子序列的8×nep個獨立分量，然后計算這8×nep個子序列分量的序間相似度，式中，nep是一個時序中的子序列個數(shù)，是分別從子序列j和k中計算出來的兩個獨立分量和之間的相關系數(shù),式中，和分別為兩個獨立分量和中n個數(shù)據(jù)點的均值。

所述步驟d中，層次聚類、偽跡的自動識別和剔除、fastica逆變換和小波系數(shù)重構的具體方法如下：

d1、利用matlab軟件信號處理工具箱中的zscore函數(shù)對特征向量f進行數(shù)據(jù)標準化；利用linkage函數(shù)計算特征向量中任意向量之間的相關距離，以距離最小準則為標準，利用linkage函數(shù)根據(jù)內(nèi)平方距離法計算聚類樹；

d2、根據(jù)步驟d1的聚類結果自動給8個獨立分量s1,s2,…,s8貼上相應的類別標簽，對偽跡分量進行置零處理，其余分量保持不變，得到剔除偽跡分量的獨立分量矩陣z(t)；

d3、利用fastica逆變換將z(t)進行投影映射，得到小波變換系數(shù)u(t)＝w^-1z(t)，利用離散小波變換工具箱中的waverec函數(shù)對u(t)進行重構，得到干凈的腦電信號。

本發(fā)明的有益效果是：

1、解決了現(xiàn)有運動想象腦電預處理方法需要訓練和偽跡參考電極，以及無法在線實現(xiàn)的問題。

2、通過離散小波變換可以將少通道(c3和c4通道)的運動想象腦電信號轉換為多通通信號用于后續(xù)fastica的分析。

3、小波變換系數(shù)比原始信號的超高斯性更強，峰度更大，從小波域進行ica在迭代算法的收斂速度、抗噪能力等方面具有顯著優(yōu)勢。

4、利用層次聚類算法基于獨立分量的時域、頻譜和序間相關性特征，可高效的自動識別和剔除多種常規(guī)偽跡。

附圖說明

圖1為本發(fā)明的在線腦電偽跡剔除方法流程圖；

圖2為腦電采集電極的放置示意圖；

圖3為被試進行左右手運動想象任務的實驗時序圖；

圖4為c3和c4通道原始腦電信號的時域圖和頻譜圖，其中，(a)為c3和c4通道原始腦電信號的時域圖，(b)為c3通道原始腦電信號的頻譜圖；

圖5為校正之后c3和c4通道腦電信號的時域圖和頻譜圖，其中，(a)為校正之后c3和c4通道腦電信號的時域圖，(b)為校正之后c3通道腦電信號的頻譜圖；

圖6為c3通道小波系數(shù)單支重構信號；

圖7為c3通道ica分解后各個獨立分量的波形圖；

圖8為c3通道ica分解后眼電分量的時域特征和頻譜特征，其中，(a)為眼電分量的時域特征，(b)為眼電分量的頻譜特征；

圖9為c3通道ica分解后肌電分量的時域波形圖和頻譜特征，其中，(a)為肌電分量的時域波形圖，(b)為肌電分量的頻譜特征；

圖10為c3通道各個分量進行層次聚類的樹狀圖；

圖11為預處理之后c3和c4通道的腦電信號波形圖。

具體實施方式

下面結合附圖和具體實施例對本發(fā)明作進一步說明。

實施例1：如圖1-11所示，一種實用化在線腦電偽跡剔除方法，實驗設計和數(shù)據(jù)說明如下：

本文提出的預處理方法一共在10名健康被試的左右手運動想象腦電中進行了測試，其中被試為6男4女；均為右利手，年齡范圍：19歲～25歲；本科及碩士學歷，均無影響腦功能的病史和腦電采集經(jīng)歷，所有被試在實驗前簽署了實驗研究知情同意書。

本次試驗的數(shù)據(jù)采集設備為16導聯(lián)eeg放大器(mipower-uc,eegcollectionv2，清華大學神經(jīng)工程實驗室)，0hz～250hz信號頻帶，采樣頻率為1000hz，24位a/d轉換器，無工頻陷波；采用國際標準導聯(lián)10-20系統(tǒng)定制的16導聯(lián)腦電帽(ag-agcl粉末電極，武漢格林泰克科技有限公司)。腦電采集電極的放置示意圖如圖2所示，實驗主要采集運動功能區(qū)c3和c4兩個電極的腦電信號，分別為x1(t)和x2(t)；左側乳突m1作為記錄參考電極，fpz為接地電極；實驗之前要確保電極與被試頭皮之間的阻抗低于5kω。

實驗在安靜的環(huán)境中進行，實驗時序圖如圖3所示，實驗過程中，被試坐在舒適的座椅上，眼睛距離電腦屏幕約1米，每次實驗開始時刻會有蜂鳴提示，同時屏幕正中央出現(xiàn)“+”，提示被試本次實驗已經(jīng)開始，該狀態(tài)持續(xù)2s；開始提示結束后“+”消失，同時屏幕正中央會隨機出現(xiàn)任務圖片提示，提示被試接下來進行哪種運動想象任務，該狀態(tài)持續(xù)1.5s；提示圖片消失之后屏幕正中央會出現(xiàn)“*”，提示被試開始進行之前圖片所提示的左手/右手運動想象任務，這個任務時間持續(xù)6s，任務結束之后會進入休息狀態(tài)，屏幕白屏，該狀態(tài)持續(xù)時間為3s。該刺激范式由廣泛應用于生物醫(yī)學工程領域的刺激軟件e-prime(版本version1.1)來完成。被試進行每次實驗之間的休息時間為5分鐘，每個被試進行5次實驗，每次實驗包括6個時序，10個被試的所有實驗均在一天內(nèi)完成。

本發(fā)明選取了被試5在第3次實驗中某一次想象右手運動時的腦電數(shù)據(jù)用于后續(xù)分析。

對于上述腦電數(shù)據(jù)結果進行在線腦電偽跡剔除的方法包括如下具體步驟：

a、對c3、c4通道的原始腦電信號進行降采樣、工頻陷波和線性漂移校正；

b、對校正之后兩個通道的腦電信號進行多尺度小波分解和小波系數(shù)單支重構；

c、分別對兩個通道的小波系數(shù)單支重構信號進行ica分解，計算各個分量的時域特征、頻譜特征以及序間相關性；

d、對兩個通道的獨立分量進行層次聚類，自動識別并剔除偽跡分量，并對剩余的獨立分量進行fastica逆變換重構和小波系數(shù)重構，得到干凈的腦電信號。

所述步驟a中的降采樣、工頻陷波和線性漂移校正具體方法如下：

a1、利用matlab信號處理工具箱中的downsample函數(shù)分別對c3、c4通道的原始腦電信號進行250hz降采樣，以降低數(shù)據(jù)總量，提高數(shù)據(jù)處理的實時性；

降采樣之后兩個通道的時域波形圖和頻譜圖如圖4所示，可以看出兩個通道的數(shù)據(jù)具有明顯的數(shù)據(jù)漂移和50hz工頻干擾現(xiàn)象；

a2、利用matlab信號處理工具箱中的dlsim函數(shù)對c3、c4通道降采樣之后的腦電數(shù)據(jù)進行50hz工頻陷波處理，剔除50hz工頻干擾；

a3、利用matlab信號處理工具箱中的detrend函數(shù)對陷波處理之后的信號進行線性校正，消除數(shù)據(jù)線性漂移帶來的偽跡；

經(jīng)過a2和a3兩步校正之后，兩個通道的時域波形圖和頻譜圖如圖5所示，可以看出兩個通道的數(shù)據(jù)已經(jīng)去除了50hz工頻干擾和數(shù)據(jù)漂移現(xiàn)象，但是數(shù)據(jù)中還存在嚴重的eog和emg等偽跡干擾；

經(jīng)過上述處理得到校正之后的腦電信號x(t)＝[x1(t),x2(t)]^t∈r^m×n，其中：x(t)是一組m維隨機向量，m為腦電采集通道數(shù)，這里m為2；t為轉置運算，x1(t)和x2(t)分別代表c3和c4通道的腦電信號；t＝1,2,…,n，為腦電信號x(t)的n個采樣點，這里n為1500；經(jīng)過步驟a處理之后得到兩個通道分別6s的腦電數(shù)據(jù)。以下分析均以c3通道為例，c4通道的處理方法和步驟同理。

所述步驟b中，多尺度小波分解和小波系數(shù)單支重構的具體步驟如下：

b1、選擇“db4”小波基，利用matlab軟件離散小波變換工具箱中的wavedec函數(shù)對c3通道的信號x1(t)進行7層離散小波分解，得到信號x1(t)的逼近系數(shù)分量ca7和細節(jié)系數(shù)分量cd1、cd2、cd3、cd4、cd5、cd6、cd7，此時，單通道的信號就變成了多通道的信號；

本發(fā)明選取了daubechies小波c3通道腦電信號的分解與重構。它通過最大頻響平坦頻率來刻畫，不同階的daubechies小波基被用在目前的一些研究中；然而，研究表明，本發(fā)明最終選取的4階daubechies小波(db4)能夠更好的正確表達eeg信號和尖峰。

b2、采用離散小波變換工具箱中的wrcoef函數(shù)對上一步分解得到的高低頻小波系數(shù)進行單支重構，得到對應于b1中小波系數(shù)ca7,cd7，cd6,…,cd1的重構信號r(t)＝[r1(t),r2(t),…,r8(t)]^t，其中，8個小波重構信號與x1(t)的數(shù)據(jù)點數(shù)相同，其波形如圖6所示。

所述步驟c中，ica分解和各個分量的時域特征、頻譜特征以及序間相關性計算方法如下：

c1、將步驟b2中的小波重構信號r(t)作為fastica的輸入，利用matlab軟件下fastica工具箱中的fastica函數(shù)，通過不斷迭代得到分離矩陣w，此處設定迭代精度為0.0001。當權值向量達到此精度時，就結束迭代過程，此時的權值即為最終的權值向量。經(jīng)過迭代，此處求出的分離矩陣w為：

r(t)經(jīng)過fastica分解后得到的8個獨立分量(s1,s2,…,s8)，其波形如圖7所示；

c2、分別對fastica分解后得到的8個獨立分量計算峰度值作為各個獨立分量的時域特征，式中，si表示第i個獨立分量，e()統(tǒng)計期望函數(shù),式中，l為獨立分量si的數(shù)據(jù)量，這里l＝1500，為獨立分量si中第l個數(shù)據(jù)點的m次方，當m為2和4時就得到和如果kurt(si)的值越大，則表示信號的動態(tài)分布越陡峭。各個獨立分量的峰度值如圖8(a)所示，從圖中可以看出，s1和s2的峰度值排在前兩位(kurt(s1)＝11.08,kurt(s2)＝8.878)，顯著大于其它分量的峰度值，s3的峰度值為-0.7788，是8個獨立分量中峰度值唯一為負值的分量，而以往的大量研究表明，腦電信號的偽跡成分中眼電偽跡往往具有最高的峰度值，噪聲偽跡作為亞高斯信號，峰度值往往為負值。

c3、分別計算fastica分解后8個獨立分量在5個頻帶范圍內(nèi)的平均頻帶功率：fspectral＝[f(δ)f(θ)f(α)f(β)f(γ)]，式中，f()是平均頻帶功率，它通過matlab軟件中的pwelch函數(shù)來計算，將平均頻帶功率作為各個獨立分量的頻譜特征；其中，s1的頻譜分布如圖8(b)所示，從圖中可以看出：s1的頻率較低，主要分布在0-12hz范圍內(nèi)，且具有較高的幅值，s2也表現(xiàn)出和s1同樣的頻譜特征；而腦電信號的偽跡成分中只有眼電偽跡往往頻率最低、幅值較大；s6的時域波形圖和頻譜分布分別如圖8(a)和8(b)所示，從圖8(b)中可以看出：s6的頻譜主要分布在30-80hz的頻帶范圍內(nèi)，且在這一范圍內(nèi)具有相對較高的幅值，而腦電信號的偽跡成分中，由于頭部轉動、手部晃動等身體運動所導致的肌電偽跡往往頻率最高(通?！?0hz)、幅值相對較大；

c4、在8個獨立分量s1,s2,…,s8中設置滑動窗時間長度為1s(即250個采樣點)的c3通道獨立分量數(shù)據(jù)作為一個子序列，一個時序中c3通道的8獨立分量最后就分成了6個子序列的48(6×8)個獨立分量，然后計算這48個子序列分量的序間相似度，式中，nep是一個時序中的子序列個數(shù)，這里nep為6，是分別從子序列j和k中計算出來的兩個獨立分量(和)之間的相關系數(shù),式中，和分別為兩個獨立分量和中n個數(shù)據(jù)點的均值，這里n為250；在這48個子序列分量中，當某些子序列分量和大多數(shù)子序列分量(這里設定為60％的子序列分量)都表現(xiàn)出極低相關系數(shù)時(其中：|·|為取絕對值運算)，取這些子序列分量中的最小相關系數(shù)作為原8個獨立分量s1,s2,…,s8中對應獨立分量的序間相似度；當某些子序列分量和大多數(shù)子序列分量(這里設定為60％的子序列分量)都表現(xiàn)出很高的相關系數(shù)時取這些子序列分量中的最大相關系數(shù)作為原8個獨立分量s1,s2,…,s8中對應獨立分量的序間相似度，從而得到一個時序中c3通道8個獨立分量s1,s2,…,s8的序間相似度特征向量。

腦電采集過程中偽跡是隨機出現(xiàn)的，難以預料，而且有時候在一些子序列中只會出現(xiàn)眼電偽跡或肌電等隨機偽跡。而包含偽跡成分的子序列分量和其它不含偽跡的子序列分量往往沒有相似的特征或相似度極低，另一方面，運動想象誘發(fā)的腦電子序列分量和其它有用腦電信號的子序列分量表現(xiàn)出極高的相似度。

通過上述計算，得到一個7維的特征向量f＝[f1,f2,f3,f4,f5,f6,f7]^t,其中f1為各個獨立分量的時域特征，f2,f3,f4,f5,f6為頻譜特征，f7為相似度特征；

所述步驟d中，對c1中得到的8個獨立分量s1,s2,…,s8進行層次聚類，自動識別并剔除偽跡，然后依次對剩余的獨立分量進行fastica逆變換重構和小波系數(shù)重構的具體方法如下：

d1、首先利用matlab軟件信號處理工具箱中的zscore函數(shù)對步驟c中計算出的特征向量f進行數(shù)據(jù)標準化，使得特征向量矩陣的平均值為零，標準差為1；然后利用linkage函數(shù)計算特征向量中任意向量之間的相關距離，這里實際計算的是兩個向量之間的歐氏距離，接下來以距離最小準則利用linkage函數(shù)根據(jù)內(nèi)平方距離法計算聚類樹，最后利用dendrogram函數(shù)根據(jù)linkage函數(shù)的輸出繪制出如圖10所示的層次聚類樹狀圖；

d2、如圖10所示，眼電分量s1和s2由于區(qū)別于其它獨立分量的時域、頻譜和序間相關性特征而被自動聚合為一類，根據(jù)以往的研究可以進一步看出：s1為眨眼偽跡，s2為眼動偽跡；肌電分量s6由于區(qū)別于其它獨立分量的頻譜和序間相似度特征而被單獨聚合為一類；正常的腦電信號分量s4、s5、s7、s8首先聚合為一類，雖然噪聲分量s3具有區(qū)別于獨特的時域特征，但是它的總體特征和腦電相似，所以最后和正常腦電聚為了一類；聚類的同時給8個獨立分量貼上相應的標簽，對偽跡分量s1、s2、s3和s6進行置零處理，其余分量保持不變，即得到無偽跡成分的獨立成分矩陣z(t)；

d3、利用fastica逆變換將z(t)進行投影映射，得到小波變換系數(shù)u(t)，即：

u(t)＝w^-1z(t)

選擇小波基函數(shù)“db4”，利用離散小波變換工具箱中的waverec函數(shù)對步驟d2中剔除了各種偽跡的小波系數(shù)u(t)重構腦電信號，如圖11所示，可以看出腦電信號中的eog偽跡、emg偽跡、50hz工頻干擾、噪聲干擾和數(shù)據(jù)漂移現(xiàn)象都已經(jīng)被很好的去除了。

為了進一步驗證本發(fā)明剔除偽跡的有益效果，分別將本發(fā)明中的方法和其它方法用于本實驗時的性能進行了對比，10名被試采用李明愛提出的dwt-ica、bogdanmijovie提出的emd-ica和本發(fā)明提出的dwicac偽跡剔除方法并計算預處理之后信號的頻帶能量特征，最后基于支持向量機(supportvectormachine,svm)分類器的分類準確率如表1所示，svm是基于統(tǒng)計學習理論的機器學習方法，適合小樣本分類問題，泛化能力強，能很好地適應運動想象腦電信號特征的單次試驗分類。由于本發(fā)明只使用了c3和cz通道的腦電數(shù)據(jù)用于分析，分類時只提取了腦電信號的頻帶能量特征，特征維數(shù)相對較低，因此本發(fā)明最后選擇了基于線性核的線性支持向量機，誤差懲罰因子取值為1。

表1采用emd-ica、dwt-ica、dwicac偽跡剔除方法并基于svm分類器的分類準確率

從表1中可以發(fā)現(xiàn)，本發(fā)明提出的dwicac方法相對于emd-ica方法和dwt-ica偽跡剔除方法而言，幾乎能自動識別并剔除腦電中的全部常規(guī)偽跡(eog、emg和噪聲等)；而在分類準確率方面該方法也顯著優(yōu)于其它兩種方法，其中：本發(fā)明對于左右手運動想象腦電的平均分類準確率為84.64％，分別高出emd-ica方法和dwt-ica方法8.23％和5.26％，最高分類準確率為89.27％，分別高出emd-ica方法和dwt-ica方法7.64％和3.55％；上述分析進一步體現(xiàn)了dwicac方法在少通道運動想象腦電在線偽跡剔除中的優(yōu)越性。

上面結合附圖對本發(fā)明的具體實施方式作了詳細說明，但是本發(fā)明并不限于上述實施方式，在本領域普通技術人員所具備的知識范圍內(nèi)，還可以在不脫離本發(fā)明宗旨的前提下作出各種變化。