第二講-EEG信號預(yù)處理與特征量

1、2021/3/261 第二講 EEG信號預(yù)處理及特征量 2021/3/262 信號預(yù)處理 v疊加平均降噪（1維信號） v信號濾波（1維信號） v小波分解（DWT,1維信號） v主成分分解（PCA,多維信號） v獨(dú)立分量分解（ICA,多維信號） vCommon Spatial Pattern（CSP,多維 信號） 2021/3/263 疊加平均降噪 x1 x2 xM M xavg xxxxN iiii N , : 12 1 vector x M xN avgi i M 1 1 1 : vector xsz iii ssss iiii N ,12 zzzz iiii N ,12 x M s M z

2、 sz avgi i M i i M avgavg 11 11 2021/3/264 疊加平均信號平均能量估計疊加平均信號平均能量估計: s,avg 2 2 11 2 1 2 2 2 1 N s N M M s avg n n N n n N s,*, 單次信號平均能量估計單次信號平均能量估計: sn n N N s 22 1 1 , 2 z 2 , 2 z ni zEn E z z iji jM iji jM i nj n, , , z 2 1 01 單次噪聲平均能量單次噪聲平均能量: 2021/3/265 No Image 信號與噪聲平均能量估計信號與噪聲平均能量估計 疊加平均噪聲平均能量

3、疊加平均噪聲平均能量: 2 , 2 avgz,navg zEn M i ni M i ni zE M z M E 1 2 , 2 2 1 , 11 11 2 z 2 z 2 M M M 單次信噪比單次信噪比:SNRi s 2 z 2 SNR M M SNR avgi s,avg 2 z,avg 2 s 2 z 2 1 疊加平均信噪比疊加平均信噪比: 2021/3/266 No Image 疊加平均在線計算 M i iavg x M x 1 1 疊加平均批量計算公式疊加平均批量計算公式: 疊加單次計算公式疊加單次計算公式: kavgavg x k kx k k kx 1 ) 1( 1 )( 20

4、21/3/267 信號濾波 v信號濾波涉及:低通、高通、帶通、陷波 v濾波器的比較: IIR滿足相同特性階數(shù)較低,只能近似線性相位,必須浮點運(yùn) 算 FIR滿足相同特性階數(shù)較高,可以做到嚴(yán)格線性相位,可以 采用整數(shù)運(yùn)算 v濾波器一般要結(jié)合實現(xiàn)時的計算效率和濾波器特性 等綜合考慮 階數(shù) 通帶、阻帶和過渡帶特性 延遲 2021/3/268 IIR濾波器 vIIR濾波器的模型為: vIIR濾波器的類型:貝塞爾、巴特沃斯、切比雪夫I型、 切比雪夫II型和橢圓型 v相位特性:貝塞爾巴特沃斯切比雪夫橢圓 v過渡帶寬度:貝塞爾巴特沃斯切比雪夫橢圓 N N M M zazaza zbzbzbb zH 2 2 1

5、 1 2 2 1 10 1 )( 2021/3/269 EEG信號波提取濾波器設(shè)計 v頻帶為8-13Hz,所以設(shè)計帶通濾波器,采樣頻 率為512,選擇4階橢圓濾波器,通帶邊緣頻率 為8.1, 12.8 B,A=ellip(4,0.5,20,8.1/256,12.8/256) h,f=freqz(B,A,1000,512); plot(f(1:100),20*log(abs(h(1:100) axis(6,15,-60,5) grid on 2021/3/2610 6789101112131415 -60 -50 -40 -30 -20 -10 0 Hz db 2021/3/2611 67891

6、01112131415 -3 -2 -1 0 1 2 3 2021/3/2612 濾波器參數(shù) vB =0.6968 -0.769432577410378 2.672998836022327 -5.3229 6.642781096781421 -5.3228 2.672998836022326 -0.769432577410377 0.6968 vA =1 -7.8686 27.0359 -53.6737 66.463174003175396 -52.8070 26.278657639283775 -7.4916 0.936634045474436 2021/3/2613 腦電信號波提取 020

7、0040006000800010000120001400016000 -2 -1 0 1 2 x 10 -4 0200040006000800010000120001400016000 -4 -2 0 2 4 x 10 -5 2021/3/2614 頻譜比較 051015202530 0 0.01 0.02 0.03 0.04 Hz 051015202530 0 0.005 0.01 0.015 Hz 2021/3/2615 連續(xù)小波變換 說明 va為尺度因子(對應(yīng)頻率),較小的a對應(yīng)高頻,較大的a對應(yīng)低頻; vb為位移因子(對應(yīng)時間); v(t)為小波母函數(shù),一般取具有單位能量的窗函數(shù); v

8、小波變換的值表示了信號f(t)與小波函數(shù)匹配的程度,例如若 對某個a和b的取值信號f(t)與小波完全相同,則小波變換為1。 dt a bt tfabafW)()(|),)( 2 1 0 2 ),)( 2 )( a da dbbafW C tf 2021/3/2616 小波波形隨尺度因子和位移因子的變化 2021/3/2617 常用連續(xù)小波函數(shù) -4-3-2-101234 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 00.050.10.150.20.25 0 10 20 30 40 50 60 70 80 -5-4-3-2-1012345 -0.4 -0

9、.2 0 0.2 0.4 0.6 0.8 1 -8-6-4-202468 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 )cos()( 0 2 2 tet t 2 2 0 2 )1)( t etct （ 82 22 2 1 )( tt eet 00.511.522.533.544.55 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 00.511.522.533.544.55 0 0.5 1 1.5 2021/3/2618 頻率突變信號的墨西哥草帽小波變換 2021/3/2619 離散正交小波變換 v離散小波變換相當(dāng)于在連續(xù)

10、小波變換公式中取a = 2- j/2,b = k 2-j/2; v小波反變換成為小波級數(shù); v小波函數(shù)必須滿足一定的條件,才能使以上變換公式和 反變換公式成立; v離散小波變換需要計算積分,不利于實際應(yīng)用,需要更 高效的算法。 dtkttfd jj j k )2()(2 22 jk j j k ktdtf)2()( 2 2021/3/2620 多分辨率分析和金子塔算法 小波函數(shù)由尺度函數(shù)確定,尺度函數(shù)一般滿足 k k ktpt)2(2)( 小波函數(shù)則可表示為 k k ktqt)2(2)( 設(shè)信號f(t)可以表示為 k MM k ktctf)2()( i j iki j k cpc 2 1 i

11、j iki j k cqd 2 1 i j iik j iik j k dqcpc 1 2 1 2 分解算法 重構(gòu)算法 注意:此算法與尺度函數(shù)和小波函數(shù)的形式無關(guān) 令:qk = (-1)k-1p-k+1 2021/3/2621 分解算法和重構(gòu)算法的含意 1、金字塔算法是對信號按頻 帶逐層分解,一直達(dá)到需要的 頻帶為止; 2、尺度函數(shù)分量為低通分量, 小波分量為各個頻帶的帶通 分量。 k j j k j ktctf)2()( 2 k jj k j ktdtg)2()( 記: 則有:fM(t)=fM-1(t)+gM-1(t) =fM-2(t)+gM-2(t)+gM-1(t)= =fM-N(t)+g

12、M-N(t)+gM-1(t) 2021/3/2622 HAAR正交小波 尺度系數(shù):p0=1, p1=1 小波系數(shù):q0=1, q1= -1 特點: 1、非零尺度系數(shù)和小波系數(shù)個數(shù)有限; 2、尺度函數(shù)和小波函數(shù)的非零區(qū)域為0,1(緊支撐); 3、尺度函數(shù)和小波函數(shù)不連續(xù),頻率窗太寬。 2021/3/2623 Db2緊支撐正交小波 00.511.522.53 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 00.511.522.53 -1.5 -1 -0.5 0 0.5 1 1.5 2 尺度系數(shù): p0=0.4829629131445341, p1=0.83651230

13、37378077 p2=0.2241438680420134, p4=-0.12943 尺度函數(shù) 小波函數(shù) 特點:1、非零尺度系數(shù)和小波系數(shù)個數(shù)有限; 2、尺度函數(shù)和小波函數(shù)的非零區(qū)域為-4,4(緊支撐); 3、尺度函數(shù)和小波函數(shù)連續(xù)。 2021/3/2624 常用小波 vHaar小波(可以看作為Daubechies小波的特 例) vDaubechies正交緊支撐小波(波形不對稱) v半正交小波(波形具有對稱性) v緊支撐雙正交小波(波形可以具有對稱性) 2021/3/2625 采用Db3對sin函數(shù)的和構(gòu)成信號的分解 2021/3/2626 采用Db2對頻率突變信號的分解 2021/3/26

14、27 采用Db5對頻率突變信號的分解 2021/3/2628 采用Db9對頻率突變信號的分解 2021/3/2629 采用Db3對用電曲線的分解 2021/3/2630 小波分解與重構(gòu)法去除基線漂移 原腦電信號 加入基線漂移后的腦電信號 去除緩慢基線漂移后的腦電信號 2021/3/2631 腦電信號的7層分解 2021/3/2632 主成分分析 v主成分分析(或稱主分量分析,principal component analysis)由皮爾遜(Pearson,1901)首先引入,后來被 霍特林(Hotelling,1933)發(fā)展了。 v主成分分析是一種通過降維技術(shù)把多個變量化為少 數(shù)幾個主成分(

15、即綜合變量)的統(tǒng)計分析方法。這些 主成分能夠反映原始變量的絕大部分信息,它們通常 表示為原始變量的某種線性組合。 2021/3/2633 尋找主成分的正交旋轉(zhuǎn) v旋轉(zhuǎn)公式: 112 212 cossin sincos yxx yxx 2021/3/2634 主成分的定義及導(dǎo)出 v設(shè) 為一個 維隨機(jī)向量,E(X)=0,其協(xié) 方差矩陣為 該矩陣為實對稱矩陣,且特征值都是非負(fù)實數(shù),設(shè) 為 ,則存在正交矩陣P,使得 p T p xxxX),( 21 2 21 2 2 212 121 2 1 ppp p p T xExxExxE xxExExxE xxExxExE XXE )( 2 1 T p T PX

16、PXEPP 0 21 p 2021/3/2635 令: ppppp p p p x x x uuu uuu uuu PX y y y Y 2 1 21 22221 11211 2 1 則有: , ,因此 Y 的任意兩個分量 不相關(guān)。Y的分量稱為X的主分量。由于 ii yVar)(0),cov( ji yy p i i T p i i p i i xVartracePPtraceyVar 111 )()()()( 總方差中屬于第 主成分 yi 的比例為 稱為主成分 yi 的貢獻(xiàn)率。 1 p ii i 2021/3/2636 v前 m 個主成分的貢獻(xiàn)率之和 稱為主成分 的累計貢獻(xiàn)率,它表明 解釋

17、的能力。 v通常取(相對于p)較小的 m,使得累計貢獻(xiàn)達(dá)到一個較 高的百分比(如8090)。此時, 可用來代替 ,從而達(dá)到降維的目的,而信息 的損失卻不多。 11 pm ii ii 12 , m y yy 12 , m y yy 12 , p x xx pm 12 , m y yy 12 , p x xx 2021/3/2637 主成分分析的步驟 v對原始信號數(shù)據(jù) pNpp N N xxx xxx xxx X 21 22221 11211 v進(jìn)行標(biāo)準(zhǔn)化處理,即 其中 i iij ij x x * N j iji x N 1 1 2 1 1 2 )( 1 ( N j iiji x N 2021/

18、3/2638 v計算相關(guān)系數(shù)矩陣計算相關(guān)系數(shù)矩陣 R=(rij) N k jkikij xx N r 1 * 1 v計算矩陣計算矩陣 R 的特征值和對應(yīng)單位正交特征向量的特征值和對應(yīng)單位正交特征向量,并并 按從大到小排列按從大到小排列: 0 21 p p QQQ, 21 v取取 ,則則 Y = P X 的每一個行向量的每一個行向量 即為主成分分量。即為主成分分量。 T p QQQP),( 21 2021/3/2639 |計算主成分貢獻(xiàn)率及累計貢獻(xiàn)率計算主成分貢獻(xiàn)率及累計貢獻(xiàn)率 貢獻(xiàn)率貢獻(xiàn)率: : 累計貢獻(xiàn)率累計貢獻(xiàn)率: : 一般取累計貢獻(xiàn)率達(dá)一般取累計貢獻(xiàn)率達(dá)8595%8595%的特征值的特征

19、值 所對應(yīng)的第一、第二、所對應(yīng)的第一、第二、第、第mm（mpmp）個主成分）個主成分 m , 21 ), 2 , 1( 1 pk p i i k k j p i ij 11 / 2021/3/2640 獨(dú)立分量分析 v主成分分析的局限性:在主成分分解 Y=PX 中,當(dāng)X不服從正態(tài) 分布時,Y的各個分量是不相關(guān)的,但不能保證是獨(dú)立的。當(dāng) X 是獨(dú)立信號的混合時,即 X = A S,主成分分析得不到 S。 v獨(dú)立分量分析的目的是:當(dāng) X = A S 時,求矩陣 W,使得 Y=WX 的各個分量獨(dú)立,此時W可能不是A的逆,但是 WA 是置換矩 陣。 v由于生物信號一般具有非平穩(wěn)、非正態(tài)等性質(zhì),因此IC

20、A比 PCA更有優(yōu)勢。 v獨(dú)立分量分解的局限性:求解ICA的計算復(fù)雜度比PCA高,理論 深奧算法復(fù)雜,各個分量需要解釋判讀。 2021/3/2641 各類ICA算法 v批處理算法:指依據(jù)一批已經(jīng)取得的數(shù)據(jù)X來進(jìn)行處理,而不是隨著數(shù)據(jù)的 不斷輸入做遞歸式處理。已有算法: 成對數(shù)據(jù)旋轉(zhuǎn)法(Jacobi法)及極大峰度法(Maxkurt法) 特征矩陣的聯(lián)合近似對角化法(JADE法) 四階盲辨識(FOBI) JADE法和Maxkurt法的混合 v自適應(yīng)算法:根據(jù)數(shù)據(jù)陸續(xù)得到而逐步更新參數(shù),使處理所得逐步趨近于期 望結(jié)果,即各分量獨(dú)立。已有算法: 常規(guī)的隨機(jī)梯度法 自然梯度與相對梯度 串行矩陣更新及其自適

21、應(yīng)算法 擴(kuò)展的Infomax法 非線性PCA自適應(yīng)法 2021/3/2642 各類ICA算法 v探查性投影追蹤:按照一定次序把各獨(dú)立分量一個一 個的逐次提取出來,每提取一個,就將該分量從原始 數(shù)據(jù)中去掉,對剩下的部分提取下一個分量。已有算 法: 梯度算法 旋轉(zhuǎn)因子乘積法 固定點算法(fastICA)-最常用算法 胎兒心電提取胎兒心電提取 8通道原始波形通道原始波形 2021/3/2644 2021/3/2645 ICA分解分解 8個分量波形個分量波形 2021/3/2646 01002003004005006007008009001000 -6 -4 -2 0 2 4 6 2021/3/264

22、7 PCA分解分解 8個分量波形個分量波形 2021/3/2648 2021/3/2649 Common Spatial Pattern (CSP) v設(shè) 和 為代表兩 個類的兩個p 維隨機(jī)向量,E(X)=0, E(Y)=0,其協(xié)方差 矩陣分別為 T p xxxX),( 21 T X XXE T p yyyY),( 21 T Y YYE 若這兩個矩陣都是正定矩陣,則存在矩陣 Q 滿足: ),( 21nX T diagQQ )1 ,1 ,1 ( 21pY T diagIQQ 實際上,Q的列向量和 為廣義特 征值問題 0 21 p 2021/3/2650 的特征向量和特征值。實際上,存在矩陣 G 滿足: T YX GG 于是 為實對稱正定矩陣,因此存在正交矩 陣 P 和非負(fù)特征值 使得: ),( 21 1 p T X T diagPGGP T XG G 1 0 21 p )( YXX 令 Q = G-TP , 則有 QQ X T IPGGGGPQQ TTT YX T 1 )( 另一方面 QQQQQQQQ Y T Y T X T YX T )( 2021/3/2651 因此有 ),( 21nX T diagQQ )1 ,1 ,1 ( 21pY T diagIQQ