Lecture-7：Matlab生物医学信号处理.ppt

《Lecture-7：Matlab生物医学信号处理.ppt》由会员分享，可在线阅读，更多相关《Lecture-7：Matlab生物医学信号处理.ppt（42页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、Biomedical Engineering6.1 多元信号分析多元信号分析主要内容主要内容一、多元分析概述一、多元分析概述二、主成分分析二、主成分分析 1 1、概述、概述 2 2、工作目的及基本性质、工作目的及基本性质 3 3、MatlabMatlab算法算法 4 4、信号处理举例、信号处理举例 三、独立成分分析三、独立成分分析 1 1、概述、概述 2 2、基本推到及性质、基本推到及性质 3 3、matlabmatlab算法算法 4 4、信号处理举例、信号处理举例Biomedical Engineering6.1 多元分析概述多元分析概述一、多元分析基本概念一、多元分析基本概念研究客研究客观

2、观事物中多个事物中多个变变量（或多个因素量（或多个因素，多个，多个测测量量值值）之）之间间相互依相互依赖赖的的统计规统计规律性。它的重要基律性。它的重要基础础之一是多元正之一是多元正态态分析分析， ，又称多元分析又称多元分析。 。 在多元分析中，多个因素或在多元分析中，多个因素或变变量通常用一个多量通常用一个多值值的矢量的矢量变变量表量表达： Mpfortxtxtxp1),(21xX包含包含M个个变变量，每一个量，每一个变变量都有量都有N个个观观察。察。X表达的多元数据可以表达的多元数据可以看做位于看做位于M维维空空间间的的变变量，每一量，每一维维都包含一个信号或者（都包含一个信号或者（图图像

3、）。多像）。多元分析研究的是元分析研究的是变变量自身及多个量自身及多个变变量之量之间间的相互关系。多元分析的的相互关系。多元分析的一个主要的任一个主要的任务务就是就是寻寻找一种找一种变换变换能能够让够让复复杂杂的多元数据尺寸减的多元数据尺寸减少或者更便于理解。少或者更便于理解。在多元信号变量里包含的一些相关信在多元信号变量里包含的一些相关信息能否由较低维数或较少的变量表达？息能否由较低维数或较少的变量表达？较少维数的多元数据变量集是不是比较少维数的多元数据变量集是不是比原始的数据集更有意义？原始的数据集更有意义？EEG信号分析信号分析,大量大量大脑皮层的信号可大脑皮层的信号可能来自于极少的区能

4、来自于极少的区域的神经源。域的神经源。Biomedical Engineering二、变换二、变换产产生新数据生新数据变变量集的量集的变换变换可以是可以是线线性的也可以是非性的也可以是非线线性的。性的。通常我通常我们们使用使用线线性性变换变换，因，因为为便于便于计计算和解算和解释释。 。 MixwtyiNjiji.11 txtxtxWtytytyMM.2121线线性性变换类变换类似于旋似于旋转变换转变换Biomedical Engineering请请看散点看散点图图How about this transformation?Biomedical Engineering6.2 主成份分析主成份分

5、析一、问题引出一、问题引出Biomedical Engineering人类基因组中的人类基因组中的DNADNA全序列是由全序列是由4 4个碱基个碱基A A，T T，C C，G G按一定顺序排成按一定顺序排成的长约的长约3030亿的序列，毫无疑问，这是一本记录着人类自身生老病死亿的序列，毫无疑问，这是一本记录着人类自身生老病死及遗传进化的全部信息的及遗传进化的全部信息的“天书天书”。但是，除了这四种碱基外，人。但是，除了这四种碱基外，人们对它所包含的内容知之甚少，如何破译这部们对它所包含的内容知之甚少，如何破译这部“天书天书”是二十一世是二十一世纪最重要的任务之一。在这个目标中，研究纪最重要的任

6、务之一。在这个目标中，研究DNADNA全序列具有什么结全序列具有什么结构，构，由这由这4 4个字符排成的看似随机的序列中隐藏着什么规律，又是解读个字符排成的看似随机的序列中隐藏着什么规律，又是解读这部天书的基础，是生物信息学（这部天书的基础，是生物信息学（BioinformaticsBioinformatics）最重要的课题）最重要的课题之一。之一。虽然人类对这部虽然人类对这部“天书天书”知之甚少，但也发现了知之甚少，但也发现了DNADNA序列中的一些序列中的一些规律性和结构。例如，在全序列中有一些是用于编码蛋白质的序列规律性和结构。例如，在全序列中有一些是用于编码蛋白质的序列片段，即由这片段

7、，即由这4 4个字符组成的个字符组成的6464种不同的种不同的3 3字符串，其中大多数用于字符串，其中大多数用于编码构成蛋白质的编码构成蛋白质的2020种氨基酸。又例如，在不用于编码蛋白质的序种氨基酸。又例如，在不用于编码蛋白质的序列片段中，列片段中，A A和和T T的含量特别多些，于是以某些碱基特别丰富作为特的含量特别多些，于是以某些碱基特别丰富作为特征去研究征去研究DNADNA序列的结构也取得了一些结果。利用统计的方法还发序列的结构也取得了一些结果。利用统计的方法还发现序列的某些片段之间具有相关性，等等。这些发现让人们相信，现序列的某些片段之间具有相关性，等等。这些发现让人们相信，DNAD

8、NA序列中存在着局部的和全局性的结构，充分发掘序列的结构对序列中存在着局部的和全局性的结构，充分发掘序列的结构对理解理解DNADNA全序列是十分有意义的。全序列是十分有意义的。 Biomedical Engineering问题：下面有20个已知类别的人工制造的序列，长度为10000，其中序列标号110 为A类，11-20为B类。请从中提取特征，使精度在90%时，长度序列可为100或以下；此外构造分类方法，并用这些已知类别的序列。（2000年“网易杯”全国大学生数学建模竞赛，DNA序列试题）http:/ Engineering例：两变量（两维）数据集：每一变量由两个正两变量（两维）数据集：每一变

9、量由两个正弦信号叠加构成。最后加入少量噪声。两个弦信号叠加构成。最后加入少量噪声。两个变量是高度相关的。变量是高度相关的。Biomedical Engineering使用使用PCAPCA解相关：寻找到某一套坐标系（旋转原来的坐标系），解相关：寻找到某一套坐标系（旋转原来的坐标系），让数据点集在其上分布沿坐标轴方向方差最大。这样使得数据让数据点集在其上分布沿坐标轴方向方差最大。这样使得数据点集应该具有针对均值处于对称形状。点集应该具有针对均值处于对称形状。两个变量的主成分；经过旋转后的新的成份量是不相关的（对称形状）思考：统计不相关是否意味独立？Biomedical Engineering ,

10、) , (),( 221 NXXX设设, , 2 , 1 ),(21nixxxiii ）（Biomedical Engineering . cossin , sincos 212211 XXFXXFBiomedical Engineering四、主成份的推到及性质四、主成份的推到及性质一、两个线性代数的结论一、两个线性代数的结论 1、若A是p阶实对称阵，则一定可以找到正交阵U，使ppp00000021AUU1pii. 2 . 1, 其中 是A A的特征根。Biomedical Engineering 2、若上述矩阵的特征根所对应的单位特征向量为 ppppppuuuuuuuuu212222111

11、211),(p1uuU 则实对称阵 属于不同特征根所对应的特征向量是正交的，即有p1uu,令AIUUUUBiomedical Engineering（一）（一） 第一主成分第一主成分设X的协方差阵为2212222111221pppppx由于x为非负定的对称阵，则有利用线性代数的知识可得，必存在正交阵U，使得p001UUX五、主成分的推导五、主成分的推导Biomedical Engineering 其中1， 2， p为x的特征根，不妨假设1 2 p 。而U恰好是由特征根相对应的特征向量所组成的正交阵。ppppppuuuuuuuuu212222111211),(p1uuUpiiiuuu，21iUi

12、Pi, 2 , 1 下面我们来看，是否由U的第一列元素所构成为原始变量的线性组合是否有最大的方差。Biomedical Engineering设有P维正交向量11111ppFa Xa X a X1211111)(aUUaaapFV121111,paaaa12p 12112p1puuau ,u ,uauBiomedical Engineeringpii121)( ua piii11auuaaUUa1aa1 1 1piiiia u u a21()piiia uBiomedical Engineering 当且仅当a1 =u1时，即 时，有最大的方差1。因为Var(F1)=U1xU1=1。 如果第一

13、主成分的信息不够，则需要寻找第二主成分。ppXuXuF11111Biomedical Engineering（二）（二） 第二主成分第二主成分在约束条件 下，寻找第二主成分 0),cov(21FFppXuXuF21122因为所以0),cov(),cov(121122121uuuuxuxuFF 则，对p维向量 ，有012uupiiipiiiiuuFV122122222)()(uuuuuu pii222)(uu22uBiomedical Engineeringpiii122uuuu2 22uUUu2 222uu 2 ppXuXuXuF22221122 所以如果取线性变换： 则 的方差次大。2F 类

14、推 ppppppppppXuXuXuFXuXuXuFXuXuXuF22112222112212211111Biomedical Engineering写为矩阵形式：XUFppppppuuuuuuuuu212222111211),(p1uuU),(21pXXXXBiomedical Engineering通常我们在matlab里使用奇异值分解SVD来取主成分（特征向量）X是数据矩阵，可分解成D（特征根的平方根），U(主成分）Biomedical Engineering六、六、Matlab实例实例一、旋转一个两周期的正弦波，旋转角45度。%Example of data rotation%Crea

15、te a two variable data set y=sin(x), % Then rotate the data set by angle of 45 degclear all, close all;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);hold on;plot(y(1,:),y(2,:), xk);%Function rotationfunction out=rot

16、ation(input,phi)r c=size (input);if rcinput=input; transpose_flag=y;endR=cos(phi),sin(phi);-sin(phi),cos(phi);out=input*R;if transpose_flag=yout=out;endBiomedical Engineering二、根据2个信号源和噪声，产生一个包含5变量的数据集。使用主成分分析，求取主成分并画出重要主成分，并画出特征根比值图。%Example of PCA analysisclear all, close all;N=1000;fs=500;w=(1:N)*

17、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);Biomedical Engineeringfigure;for i=1:5D(i,:)=D

18、(i,:)-mean(D(i,:);endU,S,pc=svd(D,0);eigen=diag(S).2;pc=pc(:,1:5);for i=1:5pc(:,i)=pc(:,i)*sqrt(eigen(i);endeigen=eigen/N;plot(eigen);total_eigen=sum(eigen);for i=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(

19、t,x-2,k,t,y+2,k);Biomedical EngineeringBiomedical Engineering6.2 独立成份分析独立成份分析 在多元统计分析中，独立成分分析或独立分量分析（Independent components analysis (ICA)） 是一种利用统计原理进行计算的方法。它是一个线性变换。这个变换把数据或信号分离成统计独立的非高斯的信号源的线性组合。独立成分分析是盲信号分离（blind source separation (BSS)）的一种特例。一、问题的引出一、问题的引出 在上述主成份分析中，我们可以注意到，把数据之间的相关性去除，并不足以使变量独立

20、，特别是当变量是非高斯分布的时候。独立主成份分析是要寻找到一种变换把原始数据转换成若干个独立的变量。独立主成份分析的主要目的是要揭示数据当中一些更有意义的变量而不是减少数据的维数。Biomedical Engineering 独立成分分析的最重要的假设就是信号源统计独立。这个假设在大多数盲信号分离的情况中符合实际情况。即使当该假设不满足时，仍然可以用独立成分分析来把观察信号统计独立化，从而进一步分析数据的特性。独立成分分析的经典问题是“鸡尾酒会问题”（cocktail party problem）。该问题描述的是给定混合信号，如何分离出鸡尾酒会中同时说话的每个人的独立信号。当有N个信号源时，通

21、常假设观察信号也有N个（例如N个MiC或者录音机）。该假设意味着混合矩阵是个方阵，即J = D，其中D是输入数据的维数，J是系统模型的维数。对于J D的情况，也分别有不同研究。这种问题非常类似于EEG信号分析。EEG信号由放置在头部的许多电极记录产生，而这些信号实际上是隐含的神经信号源组成生产的。Biomedical EngineeringmixedICA SeparationBiomedical EngineeringICA和PCA信号分析的不同在于PCA仅使用的是信号的二阶统计量，而ICA使用的是更高阶的统计量。易知，高斯分布的信号二阶以上的统计距为0，而很多信号是非高斯分布的，这样就会有

22、更高阶的统计距存在。而ICA正好很好的利用了这些高阶统计特性。二、二、ICAICA的系统模型的系统模型ICA通常假设被测信号是由一组独立信号源瞬时线性组成产生的。注意有N个变量，就有N个方程。通常我们考虑信号的顺序与时间无关:s 与 x 是无关的时间函数。Asx Biomedical Engineerings是信号源矢量，A是混合矩阵，x是被测信号矢量。这一模型也称为隐变量模型。ICA就是要求解混合矩阵A，使得xAs1但很多时候我们不能准确获得混合矩阵A，这时我们通常采用估计的策略。此外，可以看出使用这种ICA模型并不能完全恢复信号源的具体数值，也不能解出信号源的正负符号、信号的阶数或者信号的

23、数值范围。Biomedical Engineering三、三、ICAICA的计算步骤的计算步骤1.中心化数据2.白化（whiten)数据；去相关，各方向方差缩放为1。上面两步可通过PCA的方法完成Biomedical Engineering3.获得原信号独立成份的估计b应是一个合适的矢量估计能够重建独立成份。为估计b，我们要建立与独立变量有关的目标函数，然后通过最优化算法来估计出b。一个直观的方法：利用中心极限定理：大量的独立分布的信号的混合最后会趋向于成为高斯分布的信号。Gauss distribution （rand）Single sinusoid at 100HZBiomedical E

24、ngineeringTwo sinusoids mixed （100 and 30 HZ)Four sinusoids mixed （100,70, 30,25 HZ)我们可以找到一种计算方法来测试数据变量的高斯或者非高斯性。在这种方法下，就可以找到某个b使得被测的数据集ici的非高斯性达到最大。常用的这样的计算方法包括：变量的四阶积累矩、负熵，互信息等。Biomedical EngineeringICA程序工作量大。目前流行的两个比较好的ICA程序FastICA和Jade可从下面地址下载：http:/www.cis.hut.fi/projects/ica/fastical/fp.htmlht

25、tp:/sig.enst.fr/cardoso/stuff.htmlFastICA，是一种交互式的程序。而我们使用Jade算法用在下面的例子。Biomedical Engineering%Example of ICA%create a mixture using three different signals mixed five ways plus noiseclear all, close all;N=1000;fs=500;w=(1:N)*2*pi/fs;t=1:N;%generate the three signals plus noises1=.75*sin(w*12)+.1*ran

26、dn(1,N);s2=sawtooth(w*5,0.5)+.1*rand(1,N);s3=pulstran(0:999),(0:5)*180,kaiser(100,3)+0.07*randn(1,N);plot(t,s1-2,k,t,s2,k,t,s3+2,k);xlabel(Time(sec);ylabel(s(t);title(original components (before mixing);四、ICA Matlab例程Biomedical Engineering%combine the 3 source in 5 ways.Define mixing matrix;A=.5 .5

27、.5;.2 .7 .7;.7 .4 .2;-.5 .2 -.6; .7 -.5 -.4;s=s1;s2;s3;X=A*s;figure;%Center datafor i=1:5X(i,:)=X(i,:)-mean(X(i,:);plot(t,X(i,:)+2*(i-1),k);hold on;end%Do PCAfigure;U,S,pc=svd(X,0);eigen=diag(S).2;plot(eigen,k);%Compute ICAnu_ICA=input(Enter the number of indepentdent components);W=jadeR(X,nuICA);ic=(W*X);figure;plot(t,ic(:,1)-4,k,t,ic(:,2),k,t,ic(:,3)+4,k);Biomedical EngineeringBiomedical EngineeringBiomedical Engineering作业：根据3个信号源和噪声，产生一个包含7变量的数据集。使用主成分分析，求取主成分并画出重要主成分，并画出特征根比值图。实验：将课件里的cocktail里的声音信号用ICA的方法把它分离出来。并计算分离信号与原始信号的MSE值。