1、第二講 EEG信號處理基礎(chǔ)主講人：謝宏主講人：謝宏信息工程學(xué)院信息工程學(xué)院vEEG信號模型vEEG信號特征EEG信號的特點(diǎn)v隨機(jī)信號v非平穩(wěn)v非線性v非高斯過程EEG信號模型v基于神經(jīng)元的生化物理模型Hodgkin and Huxley 模型MorrisLecar 模型v基于隨機(jī)信號的動態(tài)模型線性模型：AR模型和ARMA模型非線性模型：GARCH模型EEG信號線性模型vAR模型)()()(1nxknyanypkk), 0()(2NnxvARMA模型qkpkkknxknyany01)()()(), 0()(2NnxvAIC準(zhǔn)則確定模型階數(shù)：即最小化如下目標(biāo)函數(shù))(2)ln(),(2,qpNqpA

2、ICqpv多維AR模型：對應(yīng)多通道)()()()(111nxknyaknyanyimijjpkjjkpkiiki誘發(fā)電位的Prony方法方法v首先建立輸入為脈沖函數(shù)的AR模型：)()() 1()(1nApnanyanypv假設(shè)有樣本值y(1),y(2),y(N), 令)()2() 1()2()(1() 1 () 1()(pNYNyNyypypyypypyY)()2() 1(NYpYpyyv解得AR模型系數(shù)：yYaaaap21v求解特征方程：0111ppppaaa 得到其p個(gè)根pjpjjeAeAeA,2121v方程的通解為：njnppnjnnjnpeAweAweAwny212211)(v求解其中

3、kjkkeAr1) 1() 1 ()0(111211121121pyyywwwrrrrrrRppppppEEG信號非線性模型vAR模型)()()(1nxknyanypkk), 0()(2NnxvARMA模型qkpkkknxknyany01)()()(), 0()(2NnxvAIC準(zhǔn)則確定模型階數(shù)：即最小化如下目標(biāo)函數(shù))(2)ln(),(2,qpNqpAICqpv多維AR模型：對應(yīng)多通道)()()()(111nxknyaknyanyimijjpkjjkpkiiki由模型化神經(jīng)元活動產(chǎn)生EEG信號vAR模型vARMA模型由模型化神經(jīng)元活動產(chǎn)生EEG信號vAR模型EEG信號常用統(tǒng)計(jì)特征v隨機(jī)變量或向

4、量統(tǒng)計(jì)量均值、方差、協(xié)方差與相關(guān)系數(shù)矩陣偏度(skewness)：峭度或峰度(kurtosis)：v隨機(jī)過程或隨機(jī)信號統(tǒng)計(jì)量自相關(guān)函數(shù)與功率譜雙譜等高階譜EEG信號常用統(tǒng)計(jì)特征v設(shè)m個(gè)電極通道采集的腦電信號表示為Tmnxnxnxnx)()()()(21) 1() 1 ()0() 1() 1 ()0() 1() 1 ()0(222111NxxxNxxxNxxxXmmm11)(1)(NkiiikxNnxEv均值v其N點(diǎn)的采樣數(shù)據(jù)為EEG信號常用統(tǒng)計(jì)特征v方差11222)(1)(NkiiiiikxNnxEv偏度113333/)(1)()(NkiiiiiiikxNnxEnxSkewv峭度(峰度)3)(

5、)(24iiiinxEnxKurtv協(xié)方差矩陣)()()()(TnxEnxnxEnxEA1)()()(1)()()(NkjjiijjiiijkxkxNnxnxEaEEG信號常用統(tǒng)計(jì)特征v相關(guān)系數(shù)矩陣v注意：以上統(tǒng)計(jì)量只涉及隨機(jī)變量取值，不涉及時(shí)序計(jì)算時(shí)要考慮實(shí)時(shí)性，如均值的計(jì)算，可以考慮采用算法)(ijrR 1)()(1)()(NkjjjiiijjjiiiijkxkxNnxnxEr) 1(1) )(11(1)(1)(2010NxNkxNNNkxNnxEiNkiNkiii)(1) 1(1)(nxnnnnniii平穩(wěn)EEG信號時(shí)序特征v自相關(guān)函數(shù)v功率譜：信號的功率譜是其自相關(guān)函數(shù)的傅里葉變換，即

6、10)()(1)()()(kNniiiiiknxnxNknxnxEkv多通道信號的自相關(guān)函數(shù)矩陣10)()(1)()()(kNnTTknxnxNknxnxEkRkkjiikkjiieknxnxEekS)()()()()()()()(injikkjiiXenxEeknxnxE2*| )(|)()()()(iiinnjiiXXXenxX平穩(wěn)EEG信號時(shí)序特征v功率譜的計(jì)算：周期圖法與線性時(shí)序模型法v周期圖法計(jì)算頻譜得到，簡單、功率譜不光滑、誤差大；v假設(shè)EEG信號的時(shí)序模型為AR(p)，即)()()(1nknxanxpkkv其功率譜為212|1 |)(pkkjkpeaS平穩(wěn)EEG信號時(shí)序特征vEE

7、G信號節(jié)律波能量提?。?設(shè)信號x(n)的傅里葉變換為X()，則有對有限個(gè)采樣值x(0),x(1),x(N-1), 其DFT為X(k), 對應(yīng)頻率kfs/N, 則有：dXnxin2022| )(|21| )(|sff2則EEG信號節(jié)律波平均能量提取算法可以考慮兩種方案。102102| )(|1| )(|NkNnkXNnx平穩(wěn)EEG信號時(shí)序特征設(shè)節(jié)律波頻帶范圍為fL,fH，腦電信號的N個(gè)采樣值x(0), x(1), , x(N-1)v方案一：v第一步：采用FFT算法計(jì)算DFTX(0),X(1),X(N-1); v第二步：計(jì)算NL=fL*N/fs, NH=fH*N/fsv第三步：v方案二v第一步：對

8、腦電信號x(n)進(jìn)行帶通fL,fH濾波得到輸出信號y(n);v第二步：HLNNkavgkXNE2| )(|1102| )(|1NnavgnyNEEEG信號濾波v信號濾波涉及：低通、高通、帶通、陷波v濾波器的比較：IIR滿足相同特性階數(shù)較低，只能近似線性相位，必須浮點(diǎn)運(yùn)算FIR滿足相同特性階數(shù)較高，可以做到嚴(yán)格線性相位，可以采用整數(shù)運(yùn)算v濾波器一般要結(jié)合實(shí)現(xiàn)時(shí)的計(jì)算效率和濾波器特性等綜合考慮階數(shù)通帶、阻帶和過渡帶特性延遲IIR濾波器vIIR濾波器的模型為：vIIR濾波器的類型：貝塞爾、巴特沃斯、切比雪夫I型、切比雪夫II型和橢圓型v相位特性：貝塞爾巴特沃斯切比雪夫橢圓v過渡帶寬度：貝塞爾巴特沃斯

9、切比雪夫0發(fā)散，a=0為常數(shù)，a3.6時(shí)進(jìn)入混沌狀態(tài)時(shí)進(jìn)入混沌狀態(tài)64正 文v分析論證過程65羅淪茲方程 其中其中x x、y y、z z為無量綱量為無量綱量，分別表征對流強(qiáng)度，對流中，分別表征對流強(qiáng)度，對流中升流與降流間的溫差和豎直方升流與降流間的溫差和豎直方向溫度分布的非線性度。任意向溫度分布的非線性度。任意給定初值，系統(tǒng)最終都會回到給定初值，系統(tǒng)最終都會回到狀態(tài)空間的特定區(qū)域內(nèi)，其吸狀態(tài)空間的特定區(qū)域內(nèi)，其吸引子具有精巧而奇特的結(jié)構(gòu)，引子具有精巧而奇特的結(jié)構(gòu)，如圖所示，表明系統(tǒng)進(jìn)入了混如圖所示，表明系統(tǒng)進(jìn)入了混沌狀態(tài)。沌狀態(tài)。( (圖為圖為LorenzLorenz混沌吸混沌吸引子引子) )

10、混沌識別 混沌識別主要包括定性和定量兩種方法，定性方法主要通過揭示混沌信號在時(shí)域或頻域中表現(xiàn)出的特殊空間結(jié)構(gòu)或頻率特性來判別，這種方法簡單直觀，但是過于籠統(tǒng)。 定量方法通過計(jì)算混沌信號奇異吸引子的特性參數(shù)來辨別混沌行為的方法。主要有兩個(gè)： (1)描述鄰近軌道發(fā)散率的Laypunov指數(shù) (2)描述吸引子維數(shù)的關(guān)聯(lián)維數(shù)和反映信息產(chǎn)生頻率的Kolmogorov熵Lyapunov指數(shù) 混沌系統(tǒng)初值敏感性是指相空間中初始距離很近的兩條軌跡會以指數(shù)速率發(fā)散，Lyapunov指數(shù)即是根據(jù)相軌跡有無擴(kuò)散運(yùn)動特征來判別系統(tǒng)的混沌特性。在相空間中，軌跡間的距離分別表現(xiàn)為線度、面積和體積。 對一維映射xn+1 =

11、 fxn,假設(shè)初始位置有兩個(gè)相鄰的狀態(tài)x0和y0=x0+, 則經(jīng)過n次迭代后，有 所以)()(11nnnnnxfyfxyE11111)()(nnnnnExfxyxfnnkknnnnnexfEEEEEEEE10012110)(Lyapunov指數(shù)vLyapunov指數(shù)定義：設(shè)相軌跡上兩點(diǎn)之間的初始距離為|E0|=，用|En|表示經(jīng)過n次迭代后該兩點(diǎn)之間的距離，則稱 為系統(tǒng)Lyapunov指數(shù)。 當(dāng) 0 時(shí)，系統(tǒng)具有混沌特征。10110)(ln1limln1limnkknnkkknxfnEEn時(shí)間序列Lyapunov指數(shù)的計(jì)算 在實(shí)際時(shí)間序列混沌識別中，通常只估計(jì)最大Lyapunov指數(shù)，下面介紹

12、一種算法: 小數(shù)據(jù)量法。 設(shè)時(shí)間序列x1,x2,xN, 嵌入維數(shù)m，時(shí)間延遲，重構(gòu)相空間后有： 其中kN-(m-1) ,選取一距離閾值，對初始點(diǎn)x(0)，選取與其距離最近的點(diǎn)不妨設(shè)為x(k0)=z0(0)，當(dāng)經(jīng)過t0次迭代后，z0(t0)與x(t0)的距離大于，令mTmkkkRxxxkx),()()1(|)0()0(|)()(|ln1000000 xztxtztl時(shí)間序列Lyapunov指數(shù)的計(jì)算 取與x(t0)的距離小于的點(diǎn)x(k1)=z1(0)，且z1(0)-x(t0)與z0(t0)-x(t0)的夾角最小，如圖所示重復(fù)以上過程，假設(shè)直到M-1步結(jié)束，此時(shí)tM-1=N-(m-1)，則有 10

13、1MkklMKolmogorov熵 Kolmogorov熵K在混沌的度量中有著相當(dāng)重要的應(yīng)用。對于規(guī)則運(yùn)動，K=0; 隨機(jī)系統(tǒng)K=無窮，若系統(tǒng)表現(xiàn)確定性混沌 ， 則 Kolmogorov熵是大于0的常數(shù)。Kolmogorov熵越大，那么信息的損失速度越大，系統(tǒng)的混沌程度越大?；煦缦嗫臻g重構(gòu)理論v相空間重構(gòu)是對一維時(shí)間序列，按照延遲時(shí)間和嵌入維數(shù)重構(gòu)一個(gè)與原動力系統(tǒng)等價(jià)的相空間;v相空間重構(gòu)的理論基礎(chǔ)是Takens相空間嵌入定理，即通過對一個(gè)m維流形上的連續(xù)流的觀測值h(t)的延遲(h(t),h(t+ ),h(t+2m)，可以將該連續(xù)流光滑嵌入到2m+1維空間(相空間)中；v相空間重構(gòu)過程中有兩

14、個(gè)參數(shù)選取特別重要: 延遲時(shí)間 和嵌入維數(shù).延遲時(shí)間間隔延遲時(shí)間間隔 的計(jì)算的計(jì)算主要方法v 線性自相關(guān)函數(shù)法v平均互信息法延遲時(shí)間間隔延遲時(shí)間間隔 的選取的選取線性自相關(guān)函數(shù)法NnnNnnNnnnxNxxxNxxxxNC11211,)(1)(1)(其中定義自相關(guān)函數(shù)為 選擇使得自相關(guān)函數(shù)C()第一次為零時(shí)的 的值。缺點(diǎn)：對嵌入維數(shù)大于2時(shí)不是最優(yōu)的。延遲時(shí)間間隔延遲時(shí)間間隔 的選取的選取平均互信息法為概率，其中定義平均互信息為對于時(shí)間序列),()(,)()(),(log),()(,21nnnnnnnnnNnnxxPxPxPxPxxPxxPIx選擇使I( ) 為第一個(gè)局部極小的為延遲時(shí)間間隔。

15、缺點(diǎn)：要對概率密度進(jìn)行估計(jì)，需要的樣本點(diǎn)多，計(jì)算復(fù)雜，且對嵌入維數(shù)大于2維不是最優(yōu)的。嵌入維數(shù)嵌入維數(shù)m的計(jì)算的計(jì)算主要方法vLyapunov指數(shù)法v關(guān)聯(lián)積分法嵌入維數(shù)嵌入維數(shù)m的計(jì)算的計(jì)算Lyapunov指數(shù)法（1）首先確定延遲 ；（2）依次遞增取嵌入維數(shù)m，計(jì)算嵌入相空間(x(n),x(n- ),x(n-(m-1)的Lyapunov指數(shù)，選取指數(shù)趨于常數(shù)的m作為嵌入相空間維數(shù)。缺點(diǎn)：當(dāng)樣本點(diǎn)數(shù)不多時(shí)可信度較低，當(dāng)嵌入維數(shù)較高時(shí)面臨“維數(shù)災(zāi)”。嵌入維數(shù)嵌入維數(shù)m的的計(jì)算v關(guān)聯(lián)積分: 定義：設(shè) 為時(shí)間序列的延時(shí)，m是嵌入維數(shù)，N是樣本數(shù)據(jù)集的大小，數(shù)據(jù)點(diǎn)個(gè)數(shù)M=N-(m-1) ，則重構(gòu)相空間中

16、嵌入時(shí)間序列Y(i)的關(guān)聯(lián)積分定義為： 關(guān)聯(lián)積分是一個(gè)累積分布函數(shù)，表示相空間中任意兩點(diǎn)之間距離小于半徑r的概率，這里點(diǎn)與點(diǎn)之間的距離用矢量之差的無窮范數(shù)表示。0,01,00,| , ( )ZijijZrdYYZMjijidrMMrC1)() 1(2)(嵌入維數(shù)m的計(jì)算關(guān)聯(lián)積分法的主要步驟(1)利用時(shí)間序列X1,X2,XN,先給定一個(gè)較小的嵌入維數(shù)m0,重構(gòu)相空間，得到新的序列Yi(2)計(jì)算關(guān)聯(lián)積分C(r)(3)對于r的某個(gè)取值范圍，吸引子的維數(shù)d與累積分布函數(shù)C(r)應(yīng)滿足對數(shù)線性關(guān)系，即d(m)=LnC(r)/Lnr,從而可用最小二乘擬合得到對應(yīng)于m0的關(guān)聯(lián)維數(shù)估計(jì)d(m0)(4)增加嵌入

17、維數(shù)m0,重新計(jì)算步驟(2)和(3),直到相應(yīng)的維數(shù)估計(jì)值d(m)不再隨著m的增加而在一定誤差范圍內(nèi)不變?yōu)橹?。缺點(diǎn)：當(dāng)樣本點(diǎn)數(shù)不多時(shí)可信度較低，當(dāng)嵌入維數(shù)較高時(shí)面臨“維數(shù)災(zāi)”。總 結(jié) 關(guān)于混沌判定，一般應(yīng)用最大Lyapunov指數(shù)，或者Kolmogorov熵，或者結(jié)合兩者判定。在計(jì)算延遲時(shí)間方法上，常用的有自相關(guān)函數(shù)法和互信息法，計(jì)算嵌入維常用的方法有G-P算法，這種方法是先計(jì)算出混沌時(shí)間序列的關(guān)聯(lián)維，然后再計(jì)算出嵌入維數(shù)。自適應(yīng)濾波噪聲消除v自適應(yīng)濾波是一種采用參數(shù)可在線調(diào)整的有限脈沖響應(yīng)濾波器進(jìn)行濾波；v相對于固定參數(shù)濾波器，自適應(yīng)濾波器有可能達(dá)到更好的去除噪聲的效果；v相對于主成分和獨(dú)立

18、分量分析，自適應(yīng)濾波算法簡單；v自適應(yīng)濾波需要有欲濾除噪聲的參考信號。)()()( )()(nnxnvnxneTrw)()()()(2)(2wrrwwrnnnnxnxETTT)()()(22nnxEneETrwwrrwwr)()()()(2)(2nnEnnxEnxETTT) 1() 1()()(Mnrnrnrnr)()()()(1nnxEnnETrrrw自適應(yīng)濾波算法一v可以考慮由每一次的采樣值修正參數(shù)w，由22)()()(nnxneTrwv為減小該式誤差，可以按照負(fù)梯度方向修正w，即)()(2) 1()() 1()(2nnendndennrwwww)()() 1()(2) 1(nnnnxnw

19、rrwv考慮到穩(wěn)定性，一般應(yīng)比較小，滿足)()(max10nnETrr自適應(yīng)濾波算法二v假設(shè)有n個(gè)樣本，考慮遺忘因子加權(quán)的誤差niinie02)()(wv令其關(guān)于w的梯度為零得：v得到niiniiedd0)()(2)(rwwniTiniiix0)()()(2rrw0wrrrniniTininiiiix11)()( 2)()( 2niinniTiniixii111)()( )()(rrrw自適應(yīng)濾波算法二v令v有)()() 1()()()(1nnniinTniTinrrRrrR)()()() 1()() 1()(1nnnnnxnnTrRrwww11)()() 1()(nnnnTrrRR)() 1

20、()() 1()()() 1() 1(1111nnnnnnnnTTrRrRrrRR)()()()(nnnxneTrw自適應(yīng)濾波總結(jié)v自適應(yīng)濾波需要確定參考信號；v需要確定參考信號的FIR濾波器階數(shù)M；v計(jì)算復(fù)雜度低；v有針對單樣本的在線算法；v在EEG信號處理中的應(yīng)用，如：對濾除眨眼的眼電干擾信號，參考信號可以考慮從FP1和FP2上??；ERP中，參考信號可以考慮各段的平均值；在去除心電干擾中，可以取心電信號作為參考信號。主成分分析(K-L變換)v信號的DFT和DCT是將信號投影或表示到固定正交函數(shù)系上的變換v主成分分析或稱K-L變換是將信號表示到一組標(biāo)準(zhǔn)正交信號系上的變換，該信號系一般由信號的

21、統(tǒng)計(jì)特性得到v該正交信號系張成的空間稱為信號子空間，其正交補(bǔ)空間稱為噪聲子空間v主成分分析或K-L變換的應(yīng)用用于分離信號與噪聲對信號的維數(shù)進(jìn)行壓縮降維主成分的定義及導(dǎo)出v設(shè) x(n)=(x1(n),x2(n),xp(n)T 為一個(gè) p 維平穩(wěn)隨機(jī)信號，E(x(n)=0，其協(xié)方差矩陣為 該矩陣為實(shí)對稱矩陣，其非負(fù)特征值設(shè)為 ，則存在正交矩陣U，使得)()()()()()()()()()()()()()()()()(2212221212121nxEnxnxEnxnxEnxnxEnxEnxnxEnxnxEnxnxEnxEnxnxERpppppT)()(21TpTnUxnUxEURU021p令：)()

22、()()()()()(2121222211121121nxnxnxuuuuuuuuunUxnyynynypppppppp則有： ， ，因此 y(n) 的任意兩個(gè)分量不相關(guān)。y(n)的分量稱為X的主分量。由于iinyVar)(0)(),(cov(nynyjipiiTpiipiixVartracePPtraceyVar111)()()()(總方差(平均能量)中屬于第 i 主成分 yi 的比例為 稱為主成分 yi 的貢獻(xiàn)率。1piii主成分的計(jì)算v對原始信號數(shù)據(jù))()2() 1 ()()2() 1 ()()2() 1 (222111NxxxNxxxNxxxXpppv由其奇異值分解，知存在 p 階正交

23、矩陣U和 N 階正交矩陣V，使得 其中VUXH0S p1Sv樣本協(xié)方差矩陣v因此，矩陣R的特征值為UUNUVVUNXXNRTTTHT2111SNii2v主分量為pSVUXYvV的前p行為正交信號，由其張成的空間為信號子空間)()2() 1 ()()2() 1 ()()2() 1 (222111NvvvNvvvNvvvVpppp主成分分析在降維中的應(yīng)用v首先對樣本協(xié)方差矩陣R計(jì)算正交變換矩陣U和特征值v計(jì)算前 m 個(gè)主成分的貢獻(xiàn)率之和v對給定的閾值(一般在0.80.95), 選取累計(jì)貢獻(xiàn)達(dá)到最小的m。v由 y(n)=Umx(n) (Um為矩陣U的前m行子矩陣)計(jì)算得到m維向量，可用來代替p維向量

24、x(n)，從而達(dá)到降維的目的，而信息的損失卻不多。11pmiiii主成分分析在降噪中的應(yīng)用v首先對樣本數(shù)據(jù)矩陣X計(jì)算正交矩陣U和Vv對實(shí)際采樣信號xi(n)計(jì)算分解)()()(1nnvwnxipkkikiNnkikiiknvnxnvnxnw1)()()(),()(v也可以采用在線最小二乘算法(LMS)計(jì)算分解，即由212)()()()()(iTipkkikiiwnvnxnvwnxniiidwndnwnw)() 1()(2)() 1()()(2) 1(nvnwnvnxnwTiTii主成分分析總結(jié)v要求信號服從高斯分布v各個(gè)主分量之間不相關(guān)，但不保證獨(dú)立獨(dú)立分量分析 v獨(dú)立分量分析的目的是：當(dāng) X

25、 = A S 時(shí)，求矩陣 W，使得 Y=WX 的各個(gè)分量獨(dú)立，此時(shí)W可能不是A的逆，但是 WA 是置換矩陣。v典型例子： “雞尾酒會”問題，從酒會嘈雜人聲中提取所關(guān)心對象的語音, 人的大腦可以很快辨出或集中聽某種需要關(guān)注聲音。)()()()()()()()()()()()(333232131332322212123132121111tsatsatsatxtsatsatsatxtsatsatsatx麥克風(fēng)1麥克風(fēng)2麥克風(fēng)3)(1tx) (2tx)(3tx11a12a13a21a22a)(1ts)(2ts)(3ts23a31a32a33a獨(dú)立分量分析主要研究結(jié)構(gòu):(1) 美國加州大學(xué)生物系，計(jì)算神

26、經(jīng)生物實(shí)驗(yàn)室，提出信息極大化( infomax )算法。 .(2) 日本Riken的數(shù)量神經(jīng)科學(xué)實(shí)驗(yàn)室，互信息極小化（minimization of mutual information MMI)采用人工神經(jīng)網(wǎng)絡(luò)優(yōu)化。 http:/www.brain.riken.jp/lab/mns/amari(3) 芬蘭赫爾辛基工業(yè)大學(xué)神經(jīng)網(wǎng)絡(luò)研究中心，提出了立足于逐次提取獨(dú)立分量的固定點(diǎn)算 法（fixed point algorithm) ：fastICA。 www.cis.hut.fi/oja www.cs.helsinki.fi/aapo.hyvarinen(4) 法國

27、學(xué)者：J.F.Cardoso，提出了JADE算法、批數(shù)據(jù)處理算法、近年來引人注意的稀疏分量分析。 http:/tsi.enst.fr/cardoso. ICA相關(guān)的基本概念vn階矩(moment)：mn=E(xn) v特征函數(shù)v第二特征函數(shù)vn階累計(jì)量(cumulant) 0)()()(nnnnddjdxxpx0)()(nnnnddjk)(log)(11mk 2122mmkdxexpxj)()(31123323mmmmk4121213224461243mmmmmmmkv對高斯型信號，二階以上的累計(jì)量都為0，因此可由一、二階統(tǒng)計(jì)特征來完整描述。vk40 超高斯, k40 亞高斯, 用k4大小作為

28、衡量信號距離高斯型程度的度量。均值 方差 峭度ICA相關(guān)的基本概念v聯(lián)合矩中最常用的是協(xié)方差；v聯(lián)合累計(jì)量最常用的是4階累積量，一般用cum(x1,x2,x3,x4)表示v聯(lián)合累積量性質(zhì)：當(dāng)x各分量相互獨(dú)立時(shí)，其互累計(jì)量恒等于0。比例性：cum(w1x1,w2x2,w3x3,w4x4)=w1w2w3w4cum(x1,x2,x3,x4)當(dāng)兩個(gè)隨機(jī)變量x, y獨(dú)立時(shí)，有k4(x+y)=k4(x)+k4(y)v隨機(jī)向量x=(x1,x2,xn)各個(gè)分量獨(dú)立當(dāng)且僅當(dāng)其聯(lián)合概率分布(密度)等于各個(gè)分量概率分布(密度)的乘積。ICA計(jì)算數(shù)據(jù)的預(yù)處理解混系統(tǒng)解混系統(tǒng)B白化白化W正交系統(tǒng)正交系統(tǒng)UZ(t)混合系

29、統(tǒng)混合系統(tǒng)Ax(t)y(t)s(t)系統(tǒng)簡圖系統(tǒng)簡圖 一般在一般在ICA求解之前，先對數(shù)據(jù)進(jìn)行白化處理，這樣可以求解之前，先對數(shù)據(jù)進(jìn)行白化處理，這樣可以使得各個(gè)分量之間是不相關(guān)的，便于進(jìn)行解混計(jì)算，如下圖使得各個(gè)分量之間是不相關(guān)的，便于進(jìn)行解混計(jì)算，如下圖所示所示樣本數(shù)據(jù)的白化v對原始信號數(shù)據(jù))()2() 1 ()()2() 1 ()()2() 1 (222111NxxxNxxxNxxxXmmmv假設(shè)數(shù)據(jù)已經(jīng)零均值化，先計(jì)算數(shù)據(jù)的協(xié)方差矩陣TxXXNR1v計(jì)算其特征值和特征向量，得到對角化分解：Rx=UTUv計(jì)算白化矩陣W： ，計(jì)算白化數(shù)據(jù)：Z=WX， 此時(shí)z的協(xié)方差矩陣是單位陣，即：Rz=I

30、。UW21fastICA-基于非高斯極大化的算法v假設(shè)數(shù)據(jù)經(jīng)過白化處理，即m維隨機(jī)向量z的均值向量為零向量，協(xié)方差矩陣為單位陣；v采用一次提取一個(gè)分量的方法，以極大化非高斯性作為優(yōu)化的目標(biāo)函數(shù)，可以考慮以下兩種度量極大化峭度(Kurtosis)的絕度值極小化負(fù)熵(Negentropy)v由于峭度：Kurt(y)=Ey4-Ey22 的魯棒性不好，一般 fastICA 算法采用負(fù)熵作為目標(biāo)函數(shù)。vfastICA算法也稱為固定點(diǎn)算法或投影尋蹤算法。牛頓迭代法v牛頓法迭代法是用于求解方程 f(x)=0 的解。v設(shè)在某一點(diǎn)Pk，對應(yīng)x坐標(biāo)為xk，其切線方程為v與x軸的交點(diǎn)為xk+1，則有迭代公式)()(

31、kkkxfxxxfy)()(1kkkkxfxfxx例.用牛頓迭代法求方程的根:0133xx解:13)(3xxxf設(shè)33)(2xxf由牛頓迭代法)()(1kkkkxfxfxx得取初值,5 .00 xx0 =0.5;x1 =0.3333333333x2 =0.3472222222x3 =0.3472963532x4 =0.3472963553331323kkkkxxxx迭代四次精度達(dá)10-8 1kx*x)(xfy kx負(fù)熵負(fù)熵v信息熵： 當(dāng)隨機(jī)變量的均值和方差一定時(shí)，其 pdf 為高斯分布時(shí)其信息熵最大。v負(fù)熵：任意 pdf 的 p(x) 和具有相同協(xié)方差矩陣的高斯分布pG(x)的差作為該 pdf

32、 非高斯程度的度量，定義為負(fù)熵dxxpxpXH)(log)()()()()(xHxHxJG J(x)0，p(x)=pG(x)：當(dāng)且僅當(dāng)Jp(x)=0代表高斯分布。v 由以上公式可知，計(jì)算負(fù)熵需要隨機(jī)變量的概率密度，在實(shí)際計(jì)算中是不方便的，因此需要給出其計(jì)算的近似公式。負(fù)熵的近似負(fù)熵的近似v高階統(tǒng)計(jì)量形式：設(shè)x零均值，方差為1(白化數(shù)據(jù))Edgeworth展開Gram-Charlier展開3434( )1( )( )( )3!4!Gp xkkHxHxpx ( ),nH x Hermite多項(xiàng)式2232344341 ( )47648J p xkkkk k2242343341 ( )431848J

33、p xkkkk k缺點(diǎn)：缺點(diǎn)：大值野點(diǎn)會引大值野點(diǎn)會引起較大誤差起較大誤差負(fù)熵的近似負(fù)熵的近似v非多項(xiàng)式函數(shù)形式： 其中G(1)(y)為奇函數(shù)，G(2)(y)為偶函數(shù)，v為與y有相同均值與方差的高斯隨機(jī)變量?？梢钥紤]如下更簡潔的形式vfastICA算法就以上式為基礎(chǔ)。設(shè)y=wTz，其關(guān)于w的梯度為v所以其極值點(diǎn)使 gradw(EG(y)=EG(wTz)z=0 2)2()2(22)1(1)()()()(vGEyGEkyGEkyJ2)()()(vGEyGExJ)(grad)()(2)(gradyGEvGEyGEyJww基于負(fù)熵的牛頓迭代公式基于負(fù)熵的牛頓迭代公式v考慮到約束wTw=1，由拉格朗日乘

34、子法得極值點(diǎn)滿足： v考慮映射： v由線性近似：v得到迭代公式：wzzwgEwFT)()()()(yGyg)()()(kkwTwFwwdwwdFFk)()(11kwTkkwFdwwdFwwkIzzzwgEdwwdFTTT)()(IzwgT)(0 wzzwgET)(IzzEzwgTT)(基于負(fù)熵的牛頓迭代公式基于負(fù)熵的牛頓迭代公式v于是v考慮到單位化，可以簡化為v由于G(y)為偶函數(shù)，其導(dǎo)數(shù)g(y)為奇函數(shù)，一般可選 )()(11kTkTkkkwzzwgEzwgEww)()()(1zzwgEwzwgEzwgETkkTkTkkTkTkkwzwgEzzwgEw)()(1332211)()2/exp(

35、)()tanh()(yygyyygyayg2322212113)()2/exp()1 ()()(tanh1 ()(yygyyygyaayg基于負(fù)熵的基于負(fù)熵的fastICA基于負(fù)熵提取一個(gè)分量的固定點(diǎn)算法（基于負(fù)熵提取一個(gè)分量的固定點(diǎn)算法（2001）（1）X去均值，進(jìn)行白化得去均值，進(jìn)行白化得Z；（2）任選）任選 w1 的初值，使的初值，使 w1的的2-范數(shù)為范數(shù)為1，k=1；（3）令）令 ，均值用樣本，均值用樣本均值；均值；（4）歸一化：）歸一化：wk+1=wk+1/|wk+1|;（5）若）若 不接近不接近1，返回步驟（，返回步驟（3），否則迭代結(jié)束），否則迭代結(jié)束（6）提取分量）提取分量該方法收斂具有三階收斂速度該方法收斂具有三階收斂速度kTkTkkwzwgEzzwgEw)()(1|1kTkwwZwyT基于負(fù)熵的基于負(fù)熵的fastICA基于負(fù)熵逐個(gè)提取多個(gè)分量的固定點(diǎn)算法基于負(fù)熵逐個(gè)提取多個(gè)分量的固定點(diǎn)算法（1）X去均值，進(jìn)行白化得去均值，進(jìn)行白化得Z；（2）m為待提取獨(dú)立分量數(shù)目，為待提取獨(dú)立分量數(shù)目，p=1；（3）任選）任選 wp,1 的初值，使的初值，使 wp,1 的的2-范數(shù)為范數(shù)為1，k=1；（4）令）令 ，均值用樣本均值；，均值用樣本均值；（5）正交化：）正交