EEG信號處理基礎(chǔ)和特征講解

1、EEG信號處理基礎(chǔ)和特征講解EEG信號模型EEG信號特征EEG信號的特點隨機信號非平穩(wěn)非線性非高斯過程EEG信號模型基于神經(jīng)元的生化物理模型Hodgkin and Huxley 模型MorrisLecar 模型基于隨機信號的動態(tài)模型線性模型：AR模型和ARMA模型非線性模型：GARCH模型EEG信號線性模型AR模型ARMA模型AIC準則確定模型階數(shù)：即最小化如下目標函數(shù)多維AR模型：對應(yīng)多通道誘發(fā)電位的Prony方法首先建立輸入為脈沖函數(shù)的AR模型：假設(shè)有樣本值y(1),y(2),y(N), 令解得AR模型系數(shù)：求解特征方程： 得到其p個根方程的通解為：求解其中EEG信號非線性模型AR模型AR

2、MA模型AIC準則確定模型階數(shù)：即最小化如下目標函數(shù)多維AR模型：對應(yīng)多通道由模型化神經(jīng)元活動產(chǎn)生EEG信號AR模型ARMA模型由模型化神經(jīng)元活動產(chǎn)生EEG信號AR模型EEG信號常用統(tǒng)計特征隨機變量或向量統(tǒng)計量均值、方差、協(xié)方差與相關(guān)系數(shù)矩陣偏度(skewness)：峭度或峰度(kurtosis)：隨機過程或隨機信號統(tǒng)計量自相關(guān)函數(shù)與功率譜雙譜等高階譜EEG信號常用統(tǒng)計特征設(shè)m個電極通道采集的腦電信號表示為均值其N點的采樣數(shù)據(jù)為EEG信號常用統(tǒng)計特征方差偏度峭度(峰度)協(xié)方差矩陣EEG信號常用統(tǒng)計特征相關(guān)系數(shù)矩陣注意：以上統(tǒng)計量只涉及隨機變量取值，不涉及時序計算時要考慮實時性，如均值的計算，可

3、以考慮采用算法平穩(wěn)EEG信號時序特征自相關(guān)函數(shù)功率譜：信號的功率譜是其自相關(guān)函數(shù)的傅里葉變換，即多通道信號的自相關(guān)函數(shù)矩陣平穩(wěn)EEG信號時序特征功率譜的計算：周期圖法與線性時序模型法周期圖法計算頻譜得到，簡單、功率譜不光滑、誤差大；假設(shè)EEG信號的時序模型為AR(p)，即其功率譜為平穩(wěn)EEG信號時序特征EEG信號節(jié)律波能量提?。?設(shè)信號x(n)的傅里葉變換為X()，則有對有限個采樣值x(0),x(1),x(N-1), 其DFT為X(k), 對應(yīng)頻率kfs/N, 則有：則EEG信號節(jié)律波平均能量提取算法可以考慮兩種方案。平穩(wěn)EEG信號時序特征設(shè)節(jié)律波頻帶范圍為fL,fH，腦電信號的N個采樣值x(

4、0), x(1), , x(N-1)方案一：第一步：采用FFT算法計算DFTX(0),X(1),X(N-1); 第二步：計算NL=fL*N/fs, NH=fH*N/fs第三步：方案二第一步：對腦電信號x(n)進行帶通fL,fH濾波得到輸出信號y(n);第二步：EEG信號濾波信號濾波涉及：低通、高通、帶通、陷波濾波器的比較：IIR滿足相同特性階數(shù)較低，只能近似線性相位，必須浮點運算FIR滿足相同特性階數(shù)較高，可以做到嚴格線性相位，可以采用整數(shù)運算濾波器一般要結(jié)合實現(xiàn)時的計算效率和濾波器特性等綜合考慮階數(shù)通帶、阻帶和過渡帶特性延遲IIR濾波器IIR濾波器的模型為：IIR濾波器的類型：貝塞爾、巴特沃

5、斯、切比雪夫I型、切比雪夫II型和橢圓型相位特性：貝塞爾巴特沃斯切比雪夫橢圓過渡帶寬度：貝塞爾巴特沃斯切比雪夫0發(fā)散，a=0為常數(shù)，a0 時，系統(tǒng)具有混沌特征。時間序列Lyapunov指數(shù)的計算 在實際時間序列混沌識別中，通常只估計最大Lyapunov指數(shù)，下面介紹一種算法: 小數(shù)據(jù)量法。 設(shè)時間序列x1,x2,xN, 嵌入維數(shù)m，時間延遲，重構(gòu)相空間后有： 其中kN-(m-1) ,選取一距離閾值，對初始點x(0)，選取與其距離最近的點不妨設(shè)為x(k0)=z0(0)，當(dāng)經(jīng)過t0次迭代后，z0(t0)與x(t0)的距離大于，令時間序列Lyapunov指數(shù)的計算 取與x(t0)的距離小于的點x(k

6、1)=z1(0)，且z1(0)-x(t0)與z0(t0)-x(t0)的夾角最小，如圖所示重復(fù)以上過程，假設(shè)直到M-1步結(jié)束，此時tM-1=N-(m-1)，則有 Kolmogorov熵 Kolmogorov熵K在混沌的度量中有著相當(dāng)重要的應(yīng)用。對于規(guī)則運動，K=0; 隨機系統(tǒng)K=無窮，若系統(tǒng)表現(xiàn)確定性混沌 ， 則 Kolmogorov熵是大于0的常數(shù)。Kolmogorov熵越大，那么信息的損失速度越大，系統(tǒng)的混沌程度越大。混沌相空間重構(gòu)理論相空間重構(gòu)是對一維時間序列，按照延遲時間和嵌入維數(shù)重構(gòu)一個與原動力系統(tǒng)等價的相空間;相空間重構(gòu)的理論基礎(chǔ)是Takens相空間嵌入定理，即通過對一個m維流形上的

7、連續(xù)流的觀測值h(t)的延遲(h(t),h(t+ ),h(t+2m)，可以將該連續(xù)流光滑嵌入到2m+1維空間(相空間)中；相空間重構(gòu)過程中有兩個參數(shù)選取特別重要: 延遲時間 和嵌入維數(shù).延遲時間間隔 的計算主要方法 線性自相關(guān)函數(shù)法平均互信息法延遲時間間隔 的選取線性自相關(guān)函數(shù)法定義自相關(guān)函數(shù)為 選擇使得自相關(guān)函數(shù)C()第一次為零時的 的值。缺點：對嵌入維數(shù)大于2時不是最優(yōu)的。延遲時間間隔 的選取平均互信息法選擇使I( ) 為第一個局部極小的為延遲時間間隔。缺點：要對概率密度進行估計，需要的樣本點多，計算復(fù)雜，且對嵌入維數(shù)大于2維不是最優(yōu)的。嵌入維數(shù)m的計算主要方法Lyapunov指數(shù)法關(guān)聯(lián)積

8、分法嵌入維數(shù)m的計算Lyapunov指數(shù)法（1）首先確定延遲 ；（2）依次遞增取嵌入維數(shù)m，計算嵌入相空間(x(n),x(n- ),x(n-(m-1)的Lyapunov指數(shù)，選取指數(shù)趨于常數(shù)的m作為嵌入相空間維數(shù)。缺點：當(dāng)樣本點數(shù)不多時可信度較低，當(dāng)嵌入維數(shù)較高時面臨“維數(shù)災(zāi)”。嵌入維數(shù)m的計算關(guān)聯(lián)積分: 定義：設(shè)t 為時間序列的延時，m是嵌入維數(shù)，N是樣本數(shù)據(jù)集的大小，數(shù)據(jù)點個數(shù)M=N-(m-1)t ，則重構(gòu)相空間中嵌入時間序列Y(i)的關(guān)聯(lián)積分定義為： 關(guān)聯(lián)積分是一個累積分布函數(shù)，表示相空間中任意兩點之間距離小于半徑r的概率，這里點與點之間的距離用矢量之差的無窮范數(shù)表示。嵌入維數(shù)m的計算關(guān)

9、聯(lián)積分法的主要步驟(1)利用時間序列X1,X2,XN,先給定一個較小的嵌入維數(shù)m0,重構(gòu)相空間，得到新的序列Yi(2)計算關(guān)聯(lián)積分C(r)(3)對于r的某個取值范圍，吸引子的維數(shù)d與累積分布函數(shù)C(r)應(yīng)滿足對數(shù)線性關(guān)系，即d(m)=LnC(r)/Lnr,從而可用最小二乘擬合得到對應(yīng)于m0的關(guān)聯(lián)維數(shù)估計d(m0)(4)增加嵌入維數(shù)m0,重新計算步驟(2)和(3),直到相應(yīng)的維數(shù)估計值d(m)不再隨著m的增加而在一定誤差范圍內(nèi)不變?yōu)橹?。缺點：當(dāng)樣本點數(shù)不多時可信度較低，當(dāng)嵌入維數(shù)較高時面臨“維數(shù)災(zāi)”。總 結(jié) 關(guān)于混沌判定，一般應(yīng)用最大Lyapunov指數(shù)，或者Kolmogorov熵，或者結(jié)合兩者

10、判定。在計算延遲時間方法上，常用的有自相關(guān)函數(shù)法和互信息法，計算嵌入維常用的方法有G-P算法，這種方法是先計算出混沌時間序列的關(guān)聯(lián)維，然后再計算出嵌入維數(shù)。自適應(yīng)濾波噪聲消除自適應(yīng)濾波是一種采用參數(shù)可在線調(diào)整的有限脈沖響應(yīng)濾波器進行濾波；相對于固定參數(shù)濾波器，自適應(yīng)濾波器有可能達到更好的去除噪聲的效果；相對于主成分和獨立分量分析，自適應(yīng)濾波算法簡單；自適應(yīng)濾波需要有欲濾除噪聲的參考信號。自適應(yīng)濾波算法一可以考慮由每一次的采樣值修正參數(shù)w，由為減小該式誤差，可以按照負梯度方向修正w，即考慮到穩(wěn)定性，一般應(yīng)比較小，滿足自適應(yīng)濾波算法二假設(shè)有n個樣本，考慮遺忘因子加權(quán)的誤差令其關(guān)于w的梯度為零得：得

11、到自適應(yīng)濾波算法二令有自適應(yīng)濾波總結(jié)自適應(yīng)濾波需要確定參考信號；需要確定參考信號的FIR濾波器階數(shù)M；計算復(fù)雜度低；有針對單樣本的在線算法；在EEG信號處理中的應(yīng)用，如：對濾除眨眼的眼電干擾信號，參考信號可以考慮從FP1和FP2上??；ERP中，參考信號可以考慮各段的平均值；在去除心電干擾中，可以取心電信號作為參考信號。主成分分析(K-L變換)信號的DFT和DCT是將信號投影或表示到固定正交函數(shù)系上的變換主成分分析或稱K-L變換是將信號表示到一組標準正交信號系上的變換，該信號系一般由信號的統(tǒng)計特性得到該正交信號系張成的空間稱為信號子空間，其正交補空間稱為噪聲子空間主成分分析或K-L變換的應(yīng)用用于

12、分離信號與噪聲對信號的維數(shù)進行壓縮降維主成分的定義及導(dǎo)出設(shè) x(n)=(x1(n),x2(n),xp(n)T 為一個 p 維平穩(wěn)隨機信號，E(x(n)=0，其協(xié)方差矩陣為 該矩陣為實對稱矩陣，其非負特征值設(shè)為 ，則存在正交矩陣U，使得令：則有： ， ，因此 y(n) 的任意兩個分量不相關(guān)。y(n)的分量稱為X的主分量。由于總方差(平均能量)中屬于第 i 主成分 yi 的比例為 稱為主成分 yi 的貢獻率。主成分的計算對原始信號數(shù)據(jù)由其奇異值分解，知存在 p 階正交矩陣U和 N 階正交矩陣V，使得 其中樣本協(xié)方差矩陣因此，矩陣R的特征值為主分量為V的前p行為正交信號，由其張成的空間為信號子空間主

13、成分分析在降維中的應(yīng)用首先對樣本協(xié)方差矩陣R計算正交變換矩陣U和特征值計算前 m 個主成分的貢獻率之和對給定的閾值(一般在0.80.95), 選取累計貢獻達到最小的m。由 y(n)=Umx(n) (Um為矩陣U的前m行子矩陣)計算得到m維向量，可用來代替p維向量x(n)，從而達到降維的目的，而信息的損失卻不多。主成分分析在降噪中的應(yīng)用首先對樣本數(shù)據(jù)矩陣X計算正交矩陣U和V對實際采樣信號xi(n)計算分解也可以采用在線最小二乘算法(LMS)計算分解，即由主成分分析總結(jié)要求信號服從高斯分布各個主分量之間不相關(guān)，但不保證獨立獨立分量分析 獨立分量分析的目的是：當(dāng) X = A S 時，求矩陣 W，使得

14、 Y=WX 的各個分量獨立，此時W可能不是A的逆，但是 WA 是置換矩陣。典型例子： “雞尾酒會”問題，從酒會嘈雜人聲中提取所關(guān)心對象的語音, 人的大腦可以很快辨出或集中聽某種需要關(guān)注聲音。麥克風(fēng)1麥克風(fēng)2麥克風(fēng)3獨立分量分析主要研究結(jié)構(gòu):(1) 美國加州大學(xué)生物系，計算神經(jīng)生物實驗室，提出信息極大化( infomax )算法。 .(2) 日本Riken的數(shù)量神經(jīng)科學(xué)實驗室，互信息極小化（minimization of mutual information MMI)采用人工神經(jīng)網(wǎng)絡(luò)優(yōu)化。 (3) 芬蘭赫爾辛基工業(yè)大學(xué)神經(jīng)網(wǎng)絡(luò)研究中心，提出了立足于逐次提取獨立分量的固定點算 法（fixed po

15、int algorithm) ：fastICA。 (4) 法國學(xué)者：，提出了JADE算法、批數(shù)據(jù)處理算法、近年來引人注意的稀疏分量分析。 . ICA相關(guān)的基本概念n階矩(moment)：mn=E(xn) 特征函數(shù)第二特征函數(shù)n階累計量(cumulant) 對高斯型信號，二階以上的累計量都為0，因此可由一、二階統(tǒng)計特征來完整描述。k40 超高斯, k40 亞高斯, 用k4大小作為衡量信號距離高斯型程度的度量。均值 方差 峭度ICA相關(guān)的基本概念聯(lián)合矩中最常用的是協(xié)方差；聯(lián)合累計量最常用的是4階累積量，一般用cum(x1,x2,x3,x4)表示聯(lián)合累積量性質(zhì)：當(dāng)x各分量相互獨立時，其互累計量恒等于

16、0。比例性：cum(w1x1,w2x2,w3x3,w4x4)=w1w2w3w4cum(x1,x2,x3,x4)當(dāng)兩個隨機變量x, y獨立時，有k4(x+y)=k4(x)+k4(y)隨機向量x=(x1,x2,xn)各個分量獨立當(dāng)且僅當(dāng)其聯(lián)合概率分布(密度)等于各個分量概率分布(密度)的乘積。ICA計算數(shù)據(jù)的預(yù)處理解混系統(tǒng)B白化W正交系統(tǒng)UZ(t)混合系統(tǒng)Ax(t)y(t)s(t)系統(tǒng)簡圖 一般在ICA求解之前，先對數(shù)據(jù)進行白化處理，這樣可以使得各個分量之間是不相關(guān)的，便于進行解混計算，如下圖所示樣本數(shù)據(jù)的白化對原始信號數(shù)據(jù)假設(shè)數(shù)據(jù)已經(jīng)零均值化，先計算數(shù)據(jù)的協(xié)方差矩陣計算其特征值和特征向量，得到對

17、角化分解：Rx=UTU計算白化矩陣W： ，計算白化數(shù)據(jù)：Z=WX， 此時z的協(xié)方差矩陣是單位陣，即：Rz=I。fastICA-基于非高斯極大化的算法假設(shè)數(shù)據(jù)經(jīng)過白化處理，即m維隨機向量z的均值向量為零向量，協(xié)方差矩陣為單位陣；采用一次提取一個分量的方法，以極大化非高斯性作為優(yōu)化的目標函數(shù)，可以考慮以下兩種度量極大化峭度(Kurtosis)的絕度值極小化負熵(Negentropy)由于峭度：Kurt(y)=Ey4-Ey22 的魯棒性不好，一般 fastICA 算法采用負熵作為目標函數(shù)。fastICA算法也稱為固定點算法或投影尋蹤算法。牛頓迭代法牛頓法迭代法是用于求解方程 f(x)=0 的解。設(shè)在

18、某一點Pk，對應(yīng)x坐標為xk，其切線方程為與x軸的交點為xk+1，則有迭代公式例.用牛頓迭代法求方程的根:解:由牛頓迭代法x0 =0.5;迭代四次精度達10-8 負熵信息熵： 當(dāng)隨機變量的均值和方差一定時，其 pdf 為高斯分布時其信息熵最大。負熵：任意 pdf 的 p(x) 和具有相同協(xié)方差矩陣的高斯分布pG(x)的差作為該 pdf 非高斯程度的度量，定義為負熵 J(x)0，p(x)=pG(x)：當(dāng)且僅當(dāng)Jp(x)=0代表高斯分布。 由以上公式可知，計算負熵需要隨機變量的概率密度，在實際計算中是不方便的，因此需要給出其計算的近似公式。負熵的近似高階統(tǒng)計量形式：設(shè)x零均值，方差為1(白化數(shù)據(jù))

19、Edgeworth展開Gram-Charlier展開缺點：大值野點會引起較大誤差負熵的近似非多項式函數(shù)形式： 其中G(1)(y)為奇函數(shù)，G(2)(y)為偶函數(shù)，v為與y有相同均值與方差的高斯隨機變量?？梢钥紤]如下更簡潔的形式fastICA算法就以上式為基礎(chǔ)。設(shè)y=wTz，其關(guān)于w的梯度為所以其極值點使 gradw(EG(y)=EG(wTz)z=0 基于負熵的牛頓迭代公式考慮到約束wTw=1，由拉格朗日乘子法得極值點滿足： 考慮映射： 由線性近似：得到迭代公式：基于負熵的牛頓迭代公式于是考慮到單位化，可以簡化為由于G(y)為偶函數(shù)，其導(dǎo)數(shù)g(y)為奇函數(shù)，一般可選 基于負熵的fastICA基于

20、負熵提取一個分量的固定點算法（2001）（1）X去均值，進行白化得Z；（2）任選 w1 的初值，使 w1的2-范數(shù)為1，k=1；（3）令 ，均值用樣本均值；（4）歸一化：wk+1=wk+1/|wk+1|;（5）若 不接近1，返回步驟（3），否則迭代結(jié)束（6）提取分量該方法收斂具有三階收斂速度基于負熵的fastICA基于負熵逐個提取多個分量的固定點算法（1）X去均值，進行白化得Z；（2）m為待提取獨立分量數(shù)目，p=1；（3）任選 wp,1 的初值，使 wp,1 的2-范數(shù)為1，k=1；（4）令 ，均值用樣本均值；（5）正交化：（6）歸一化：wp,k+1=wp,k+1/|wp,k+1|;（7）若 不接近1，k=k+1,返回步驟（4）;（8）若p0，取gi(yi)=-2tanh(yi)，否則g(yi)=tanh(yi)-yi（5）更新分離矩陣BB+I+g(y)yTB（6）如果尚未收斂