Lecture-7：Matlab生物醫(yī)學(xué)信號處理

Lecture-7：Matlab生物醫(yī)學(xué)信號處理第一頁，共41頁。6.1多元信號分析[主要內(nèi)容]一、多元分析概述二、主成分分析1、概述2、工作目的及基本性質(zhì)3、Matlab算法4、信號處理舉例

三、獨(dú)立成分分析1、概述2、基本推到及性質(zhì)3、matlab算法4、信號處理舉例第二頁，共41頁。多元分析概述一、多元分析基本概念研究客觀事物中多個變量（或多個因素，多個測量值）之間相互依賴的統(tǒng)計(jì)規(guī)律性。它的重要基礎(chǔ)之一是多元正態(tài)分析，又稱多元分析。

在多元分析中，多個因素或變量通常用一個多值的矢量變量表達(dá)：X包含M個變量，每一個變量都有N個觀察。X表達(dá)的多元數(shù)據(jù)可以看做位于M維空間的變量，每一維都包含一個信號或者（圖像）。多元分析研究的是變量自身及多個變量之間的相互關(guān)系。多元分析的一個主要的任務(wù)就是尋找一種變換能夠讓復(fù)雜的多元數(shù)據(jù)尺寸減少或者更便于理解。在多元信號變量里包含的一些相關(guān)信息能否由較低維數(shù)或較少的變量表達(dá)？較少維數(shù)的多元數(shù)據(jù)變量集是不是比原始的數(shù)據(jù)集更有意義？EEG信號分析,大量大腦皮層的信號可能來自于極少的區(qū)域的神經(jīng)源。第三頁，共41頁。二、變換產(chǎn)生新數(shù)據(jù)變量集的變換可以是線性的也可以是非線性的。通常我們使用線性變換，因?yàn)楸阌谟?jì)算和解釋。線性變換類似于旋轉(zhuǎn)變換第四頁，共41頁。請看散點(diǎn)圖含兩個變量的數(shù)據(jù)集：線性變換前（left）var(x1)=0.34,var(x2)=0.20;變換后Howaboutthistransformation?第五頁，共41頁。6.2主成份分析一、問題引出美國的統(tǒng)計(jì)學(xué)家斯通(stone)在1947年關(guān)于國民經(jīng)濟(jì)的研究。他曾利用美國1929一1938年各年的數(shù)據(jù)，得到了17個反映國民收入與支出的變量要素，例如雇主補(bǔ)貼、消費(fèi)資料和生產(chǎn)資料、純公共支出、凈增庫存、股息、利息外貿(mào)平衡等等。在進(jìn)行主成分分析后，竟以％的精度，用三新變量就取代了原17個變量。根據(jù)經(jīng)濟(jì)學(xué)知識，斯通給這三個新變量分別命名為總收入F1、總收入變化率F2和經(jīng)濟(jì)發(fā)展或衰退的趨勢F3。更有意思的是，這三個變量其實(shí)都是可以直接測量的。第六頁，共41頁。人類基因組中的DNA全序列是由4個堿基A，T，C，G按一定順序排成的長約30億的序列，毫無疑問，這是一本記錄著人類自身生老病死及遺傳進(jìn)化的全部信息的“天書”。但是，除了這四種堿基外，人們對它所包含的內(nèi)容知之甚少，如何破譯這部“天書”是二十一世紀(jì)最重要的任務(wù)之一。在這個目標(biāo)中，研究DNA全序列具有什么結(jié)構(gòu)，由這4個字符排成的看似隨機(jī)的序列中隱藏著什么規(guī)律，又是解讀這部天書的基礎(chǔ)，是生物信息學(xué)（Bioinformatics）最重要的課題之一。雖然人類對這部“天書”知之甚少，但也發(fā)現(xiàn)了DNA序列中的一些規(guī)律性和結(jié)構(gòu)。例如，在全序列中有一些是用于編碼蛋白質(zhì)的序列片段，即由這4個字符組成的64種不同的3字符串，其中大多數(shù)用于編碼構(gòu)成蛋白質(zhì)的20種氨基酸。又例如，在不用于編碼蛋白質(zhì)的序列片段中，A和T的含量特別多些，于是以某些堿基特別豐富作為特征去研究DNA序列的結(jié)構(gòu)也取得了一些結(jié)果。利用統(tǒng)計(jì)的方法還發(fā)現(xiàn)序列的某些片段之間具有相關(guān)性，等等。這些發(fā)現(xiàn)讓人們相信，DNA序列中存在著局部的和全局性的結(jié)構(gòu)，充分發(fā)掘序列的結(jié)構(gòu)對理解DNA全序列是十分有意義的。第七頁，共41頁。問題：下面有20個已知類別的人工制造的序列，長度為10000，其中序列標(biāo)號1—10為A類，11-20為B類。請從中提取特征，使精度在90%時，長度序列可為100或以下；此外構(gòu)造分類方法，并用這些已知類別的序列。（2000年“網(wǎng)易杯”全國大學(xué)生數(shù)學(xué)建模競賽，DNA序列試題）二、主成分分析的工作目標(biāo)及基本原理亦稱主分量分析,是把各變量之間互相關(guān)聯(lián)的復(fù)雜關(guān)系進(jìn)行簡化分析的方法。主成分分析的工作目標(biāo)數(shù)學(xué)本質(zhì)——利用降維技術(shù)，將多個變量（指標(biāo)）歸結(jié)為線性無關(guān)的幾個（少數(shù)）主成分（綜合指標(biāo)）.工作目標(biāo)——化簡多指標(biāo)系統(tǒng)，構(gòu)造方便系統(tǒng)分析的少數(shù)綜合指標(biāo).在力求數(shù)據(jù)信息丟失最少的原則下，對高維的變量空間降維，即研究指標(biāo)體系的少數(shù)幾個線性組合，并且這幾個線性組合所構(gòu)成的綜合指標(biāo)將盡可能多地保留原來指標(biāo)變異方面的信息。這些綜合指標(biāo)就稱為主成分第八頁，共41頁。例：兩變量（兩維）數(shù)據(jù)集：每一變量由兩個正弦信號疊加構(gòu)成。最后加入少量噪聲。兩個變量是高度相關(guān)的。第九頁，共41頁。使用PCA解相關(guān)：尋找到某一套坐標(biāo)系（旋轉(zhuǎn)原來的坐標(biāo)系），讓數(shù)據(jù)點(diǎn)集在其上分布沿坐標(biāo)軸方向方差最大。這樣使得數(shù)據(jù)點(diǎn)集應(yīng)該具有針對均值處于對稱形狀。兩個變量的主成分；經(jīng)過旋轉(zhuǎn)后的新的成份量是不相關(guān)的（對稱形狀）思考：統(tǒng)計(jì)不相關(guān)是否意味獨(dú)立？第十頁，共41頁。F

2F

1X

1X

2三、主成分的幾何解釋

以最簡單的二元正態(tài)變量來說明主成分的幾何意義

.其

n個樣本點(diǎn)的散布大致為一個橢圓.n個點(diǎn)的在平面上作一個坐標(biāo)變換，即按坐標(biāo)

X1

和

X2

呈現(xiàn)某種線性相關(guān)性

.逆時針方向旋轉(zhuǎn)θ角度

.在坐標(biāo)系

X1

OX2

中，取新坐標(biāo)軸，在橢圓長軸方向取F1

,

短軸方向取F2

.第十一頁，共41頁。根據(jù)旋軸變換公式新老坐標(biāo)之間有關(guān)系

n個點(diǎn)的坐標(biāo)

F1

和

F2

幾乎不相關(guān).在坐標(biāo)系

F1

OF2

中，F(xiàn)

2F

1X

1X

2在

F1

軸上的方差達(dá)到最大，在此方向上所含的有關(guān)

n個樣品間差異的信息是最多的

，故稱

F1

為

第一主成分

.在和

F1

正交的軸

F2

上方差較小，稱

F2

為

第二主成分

.因此，用一維空間代替二維空間時，選用F1

可使信息的損失降到最小.這種系統(tǒng)簡化方法體現(xiàn)了抓事物主要矛盾的哲學(xué)思維.θ=？第十二頁，共41頁。四、主成份的推到及性質(zhì)一、兩個線性代數(shù)的結(jié)論

1、若A是p階實(shí)對稱陣，則一定可以找到正交陣U，使其中

是A的特征根。第十三頁，共41頁。

2、若上述矩陣的特征根所對應(yīng)的單位特征向量為

則實(shí)對稱陣屬于不同特征根所對應(yīng)的特征向量是正交的，即有令第十四頁，共41頁。（一）第一主成分設(shè)X的協(xié)方差陣為由于Σx為非負(fù)定的對稱陣，則有利用線性代數(shù)的知識可得，必存在正交陣U，使得五、主成分的推導(dǎo)第十五頁，共41頁。

其中1，2，…，p為Σx的特征根，不妨假設(shè)12

…p

。而U恰好是由特征根相對應(yīng)的特征向量所組成的正交陣。

下面我們來看，是否由U的第一列元素所構(gòu)成為原始變量的線性組合是否有最大的方差。第十六頁，共41頁。設(shè)有P維正交向量第十七頁，共41頁。第十八頁，共41頁。

當(dāng)且僅當(dāng)a1=u1時，即時，有最大的方差1。因?yàn)閂ar(F1)=U’1xU1=1。

如果第一主成分的信息不夠，則需要尋找第二主成分。第十九頁，共41頁。（二）第二主成分在約束條件下，尋找第二主成分因?yàn)樗?/p>

則，對p維向量，有第二十頁，共41頁。

所以如果取線性變換：

則的方差次大。

類推

第二十一頁，共41頁。寫為矩陣形式：第二十二頁，共41頁。通常我們在matlab里使用奇異值分解SVD來取主成分（特征向量）X是數(shù)據(jù)矩陣，可分解成D（特征根的平方根），U(主成分）第二十三頁，共41頁。六、Matlab實(shí)例一、旋轉(zhuǎn)一個兩周期的正弦波，旋轉(zhuǎn)角45度。%Exampleofdatarotation%Createatwovariabledatasety=sin(x),%Thenrotatethedatasetbyangleof45degclearall,closeall;N=100;x(1,:)=(1:N)/10;x(2,:)=sin(x(1,:)*4*pi/10);plot(x(1,:),x(2,:),'*k');xlabel('x1');ylabel('x2');phi=45*(2*pi/360);y=rotation(x,phi);holdon;plot(y(1,:),y(2,:),'xk');%Functionrotationfunctionout=rotation(input,phi)[rc]=size(input);ifr<cinput=input';transpose_flag='y';endR=[cos(phi),sin(phi);-sin(phi),cos(phi)];out=input*R;iftranspose_flag=='y'out=out';end第二十四頁，共41頁。二、根據(jù)2個信號源和噪聲，產(chǎn)生一個包含5變量的數(shù)據(jù)集。使用主成分分析，求取主成分并畫出重要主成分，并畫出特征根比值圖。%ExampleofPCAanalysisclearall,closeall;N=1000;fs=500;w=(1:N)*2*pi/fs;t=1:N;x=0.75*sin(w*5);y=sawtooth(w*7,0.5);D(1,:)=.5*y+.5*x+.1*rand(1,N);D(2,:)=.2*y+.7*x+.1*rand(1,N);D(3,:)=.7*y+.2*x+.1*rand(1,N);D(4,:)=-.6*y+-.24*x+.2*rand(1,N);D(5,:)=.6*rand(1,N);plot(t,D(1,:)+0,t,D(2,:)+2,t,D(3,:)+4,t,D(4,:)+6,t,D(5,:)+8);第二十五頁，共41頁。figure;fori=1:5D(i,:)=D(i,:)-mean(D(i,:));end[U,S,pc]=svd(D,0);eigen=diag(S).^2;pc=pc(:,1:5);fori=1:5pc(:,i)=pc(:,i)*sqrt(eigen(i));endeigen=eigen/N;plot(eigen);total_eigen=sum(eigen);fori=1:5pct(i)=sum(eigen(i:5))/total_eigen;enddisp(pct*100);S=cov(pc);figure;subplot(1,2,1);plot(t,pc(:,1)-2,t,pc(:,2)+2);subplot(1,2,2);plot(t,x-2,'k',t,y+2,'k');第二十六頁，共41頁。第二十七頁，共41頁。6.2獨(dú)立成份分析在多元統(tǒng)計(jì)分析中，獨(dú)立成分分析或獨(dú)立分量分析（Independentcomponentsanalysis(ICA)）是一種利用統(tǒng)計(jì)原理進(jìn)行計(jì)算的方法。它是一個線性變換。這個變換把數(shù)據(jù)或信號分離成統(tǒng)計(jì)獨(dú)立的非高斯的信號源的線性組合。獨(dú)立成分分析是盲信號分離（blindsourceseparation(BSS)）的一種特例。一、問題的引出在上述主成份分析中，我們可以注意到，把數(shù)據(jù)之間的相關(guān)性去除，并不足以使變量獨(dú)立，特別是當(dāng)變量是非高斯分布的時候。獨(dú)立主成份分析是要尋找到一種變換把原始數(shù)據(jù)轉(zhuǎn)換成若干個獨(dú)立的變量。獨(dú)立主成份分析的主要目的是要揭示數(shù)據(jù)當(dāng)中一些更有意義的變量而不是減少數(shù)據(jù)的維數(shù)。第二十八頁，共41頁。

獨(dú)立成分分析的最重要的假設(shè)就是信號源統(tǒng)計(jì)獨(dú)立。這個假設(shè)在大多數(shù)盲信號分離的情況中符合實(shí)際情況。即使當(dāng)該假設(shè)不滿足時，仍然可以用獨(dú)立成分分析來把觀察信號統(tǒng)計(jì)獨(dú)立化，從而進(jìn)一步分析數(shù)據(jù)的特性。獨(dú)立成分分析的經(jīng)典問題是“雞尾酒會問題”（cocktailpartyproblem）。該問題描述的是給定混合信號，如何分離出雞尾酒會中同時說話的每個人的獨(dú)立信號。當(dāng)有N個信號源時，通常假設(shè)觀察信號也有N個（例如N個MiC或者錄音機(jī)）。該假設(shè)意味著混合矩陣是個方陣，即J=D，其中D是輸入數(shù)據(jù)的維數(shù)，J是系統(tǒng)模型的維數(shù)。對于J<D和J>D的情況，也分別有不同研究。這種問題非常類似于EEG信號分析。EEG信號由放置在頭部的許多電極記錄產(chǎn)生，而這些信號實(shí)際上是隱含的神經(jīng)信號源組成生產(chǎn)的。第二十九頁，共41頁。mixedICASeparation第三十頁，共41頁。ICA和PCA信號分析的不同在于PCA僅使用的是信號的二階統(tǒng)計(jì)量，而ICA使用的是更高階的統(tǒng)計(jì)量。易知，高斯分布的信號二階以上的統(tǒng)計(jì)距為0，而很多信號是非高斯分布的，這樣就會有更高階的統(tǒng)計(jì)距存在。而ICA正好很好的利用了這些高階統(tǒng)計(jì)特性。二、ICA的系統(tǒng)模型ICA通常假設(shè)被測信號是由一組獨(dú)立信號源瞬時線性組成產(chǎn)生的。注意有N個變量，就有N個方程。通常我們考慮信號的順序與時間無關(guān):s與x是無關(guān)的時間函數(shù)。第三十一頁，共41頁。s是信號源矢量，A是混合矩陣，x是被測信號矢量。這一模型也稱為隱變量模型。ICA就是要求解混合矩陣A，使得但很多時候我們不能準(zhǔn)確獲得混合矩陣A，這時我們通常采用估計(jì)的策略。此外，可以看出使用這種ICA模型并不能完全恢復(fù)信號源的具體數(shù)值，也不能解出信號源的正負(fù)符號、信號的階數(shù)或者信號的數(shù)值范圍。第三十二頁，共41頁。三、ICA的計(jì)算步驟1.中心化數(shù)據(jù)2.白化（whiten)數(shù)據(jù)；去相關(guān)，各方向方差縮放為1。上面兩步可通過PCA的方法完成第三十三頁，共41頁。3.獲得原信號獨(dú)立成份的估計(jì)b應(yīng)是一個合適的矢量估計(jì)能夠重建獨(dú)立成份。為估計(jì)b，我們要建立與獨(dú)立變量有關(guān)的目標(biāo)函數(shù)，然后通過最優(yōu)化算法來估計(jì)出b。一個直觀的方法：利用中心極限定理：大量的獨(dú)立分布的信號的混合最后會趨向于成為高斯分布的信號。Gaussdistribution（rand）Singlesinusoidat100HZ第三十四頁，共41頁。Twosinusoidsmixed（100and30HZ)Foursinusoidsmixed（100,70,30,25HZ)我們可以找到一種計(jì)算方法來測試數(shù)據(jù)變量的高斯或者非高斯性。在這種方法下，就可以找到某個b使得被測的數(shù)據(jù)集ici的非高斯性達(dá)到最大。常用的這樣的計(jì)算方法包括：變量的四階積累矩、負(fù)熵，互信息等。第三十五頁，共41頁。ICA程序工作量大。目前流行的兩個比較好的ICA程序FastICA和Jade可從下面地址下載：FastICA，是一種交互式的程序。而我們使用Jade算法用在下面的例子。第三十六頁，共41頁。%ExampleofICA%createamixtureusingthreedifferentsignalsmixedfivewaysplusnoiseclearall,closeall;N=1000;fs=500;w=(1:N)*2*pi/fs;t=1:N;%generatethethreesignalsplusnoises1=.75*sin(w*12)+.1*randn(1,