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

1、计算机仿真,Lecture 7：Matlab生物医学信号分析,刘 恒 博士 ,西南科技大学信息工程学院,6.1 多元信号分析,主要内容,一、多元分析概述 二、主成分分析 1、概述 2、工作目的及基本性质 3、Matlab算法 4、信号处理举例 三、独立成分分析 1、概述 2、基本推到及性质 3、matlab算法 4、信号处理举例,6.1 多元分析概述,一、多元分析基本概念,研究客观事物中多个变量（或多个因素，多个测量值）之间相互依赖的 统计规律性。它的重要基础之一是多元正态分析，又称多元分析。,在多元分析中，多个因素或变量通常用一个多值的矢量变量表达：,X包含M个变量，每一个变量都有N个观察。

2、X表达的多元数据可以看做位于M维空间的变量，每一维都包含一个信号或者（图像）。多元分析研究的是变量自身及多个变量之间的相互关系。多元分析的一个主要的任务就是寻找一种变换能够让复杂的多元数据尺寸减少或者更便于理解。,在多元信号变量里包含的一些相关信息能否由较低维数或较少的变量表达？较少维数的多元数据变量集是不是比原始的数据集更有意义？,EEG信号分析,大量大脑皮层的信号可能来自于极少的区域的神经源。,二、变换,产生新数据变量集的变换可以是线性的也可以是非线性的。通常我们使用线性变换，因为便于计算和解释。,线性变换类似于旋转变换,请看散点图,含两个变量的数据集：线性变换前（left）var(x1)

3、=0.34, var(x2)=0.20; 变换后(right) ar(x1)=0.50,var(x2)=0.005,How about this transformation?,6.2 主成份分析,一、问题引出,美国的统计学家斯通(stone)在1947年关于国民经济的研究。他曾利用美国1929一1938年各年的数据，得到了17个反映国民收入与支出的变量要素，例如雇主补贴、消费资料和生产资料、纯公共支出、净增库存、股息、利息外贸平衡等等。,在进行主成分分析后，竟以97.4的精度，用三新变量就取代了原17个变量。根据经济学知识，斯通给这三个新变量分别命名为总收入F1、总收入变化率F2和经济发展或

4、衰退的趋势F3。更有意思的是，这三个变量其实都是可以直接测量的。,人类基因组中的DNA全序列是由4个碱基A，T，C，G按一定顺序排成的长约30亿的序列，毫无疑问，这是一本记录着人类自身生老病死及遗传进化的全部信息的“天书”。但是，除了这四种碱基外，人们对它所包含的内容知之甚少，如何破译这部“天书”是二十一世纪最重要的任务之一。在这个目标中，研究DNA全序列具有什么结构， 由这4个字符排成的看似随机的序列中隐藏着什么规律，又是解读这部天书的基础，是生物信息学（Bioinformatics）最重要的课题之一。,虽然人类对这部“天书”知之甚少，但也发现了DNA序列中的一些规律性和结构。例如，在全序列

5、中有一些是用于编码蛋白质的序列片段，即由这4个字符组成的64种不同的3字符串，其中大多数用于编码构成蛋白质的20种氨基酸。又例如，在不用于编码蛋白质的序列片段中，A和T的含量特别多些，于是以某些碱基特别丰富作为特征去研究DNA序列的结构也取得了一些结果。利用统计的方法还发现序列的某些片段之间具有相关性，等等。这些发现让人们相信，DNA序列中存在着局部的和全局性的结构，充分发掘序列的结构对理解DNA全序列是十分有意义的。,问题：下面有20个已知类别的人工制造的序列，长度为10000，其中序列标号110 为A类，11-20为B类。请从中提取特征，使精度在90%时，长度序列可为100或以下；此外构造

6、分类方法，并用这些已知类别的序列。（2000年“网易杯”全国大学生数学建模竞赛，DNA序列试题）,二、主成分分析的工作目标及基本原理,亦称主分量分析,是把各变量之间互相关联的复杂关系进行简化分析的方法。, 主成分分析的工作目标,数学本质利用降维技术，将多个变量（指标）归结为线性无关的几个（少数）主成分（综合指标）.,工作目标化简多指标系统，构造方便系统分析的少数综合指标.在力求数据信息丢失最少的原则下，对高维的变量空间降维，即研究指标体系的少数几个线性组合，并且这几个线性组合所构成的综合指标将尽可能多地保留原来指标变异方面的信息。这些综合指标就称为主成分,例：,两变量（两维）数据集：每一变量由

7、两个正弦信号叠加构成。最后加入少量噪声。两个变量是高度相关的。,使用PCA解相关：寻找到某一套坐标系（旋转原来的坐标系），让数据点集在其上分布沿坐标轴方向方差最大。这样使得数据点集应该具有针对均值处于对称形状。,两个变量的主成分；经过旋转后的新的成份量是不相关的（对称形状）,思考：统计不相关是否意味独立？,三、 主成分的几何解释,以最简单的二元正态变量来说明主成分的几何意义 .,其 n 个样本点,的散布大致为一个椭圆.,n 个点的,在平面上作一个坐标变换，即按,坐标 X1 和 X2 呈现某种线性相关性 .,逆时针方向旋转角度 .,在坐标系 X1 OX2 中，,取新坐标轴，在椭圆长轴方向取F1

8、, 短轴方向取F2 .,根据旋轴变换公式新老坐标之间有关系,n 个点的坐标 F1 和 F2 几乎不相关.,在坐标系 F1 OF2 中，,在 F1 轴上的方差达到最大，在此方向上所含的有关 n 个,样品间差异的信息是最多的 ，故,称 F1 为 第一主成分 .,在和 F1 正交的轴 F2 上方差较,小，称 F2 为 第二主成分 .,因此，用一维空间代替二维空,间时，选用 F1 可使信息的损失降到最小.,这种系统简化方法体现了抓事物主要矛盾的哲学思维.,=？,四、主成份的推到及性质,一、两个线性代数的结论,1、若A是p阶实对称阵，则一定可以找到正交阵U，使,其中 是A的特征根。,2、若上述矩阵的特征

9、根所对应的单位特征向量为,则实对称阵 属于不同特征根所对应的特征向量是正交的，即有,令,（一） 第一主成分,设X的协方差阵为,由于x为非负定的对称阵，则有利用线性代数的知识可得，必存在正交阵U，使得,五、主成分的推导,其中1， 2， p为x的特征根，不妨假设1 2 p 。而U恰好是由特征根相对应的特征向量所组成的正交阵。,下面我们来看，是否由U的第一列元素所构成为原始 变量的线性组合是否有最大的方差。,设有P维正交向量,当且仅当a1 =u1时，即 时，有最大的方差1。因为Var(F1)=U1xU1=1。 如果第一主成分的信息不够，则需要寻找第二主成分。,（二） 第二主成分,在约束条件 下，寻找

10、第二主成分,因为 所以,则，对p维向量 ，有,所以如果取线性变换：,则 的方差次大。,类推,写为矩阵形式：,通常我们在matlab里使用奇异值分解SVD来取主成分（特征向量）,X是数据矩阵，可分解成D（特征根的平方根），U(主成分）,六、Matlab实例,一、旋转一个两周期的正弦波，旋转角45度。,%Example of data rotation %Create a two variable data set y=sin(x), % Then rotate the data set by angle of 45 deg clear all, close all; N=100; x(1,:)=

11、(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 rotation function out=rotation(input,phi) r c=size (input); if rc input=input; transpose_flag=y; end R=cos(phi),sin(phi);-sin(phi),c

12、os(phi); out=input*R; if transpose_flag=y out=out; end,二、根据2个信号源和噪声，产生一个包含5变量的数据集。使用主成分分析，求取主成分并画出重要主成分，并画出特征根比值图。,%Example of PCA analysis clear all, close all; 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)

13、; 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);,figure; for i=1:5 D(i,:)=D(i,:)-mean(D(i,:); end U,S,pc=svd(D,0); eigen=diag(S).2; pc=pc(:,1:5); for i=1:5 pc(:,i)=pc(:,i)*sqrt(eigen(i); end eigen=ei

14、gen/N; plot(eigen); total_eigen=sum(eigen); for i=1:5 pct(i)=sum(eigen(i:5)/total_eigen; end disp(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);,6.2 独立成份分析,在多元统计分析中，独立成分分析或独立分量分析（Independent components analysis (ICA)） 是一种利用统计原理进行计算的方

15、法。它是一个线性变换。这个变换把数据或信号分离成统计独立的非高斯的信号源的线性组合。独立成分分析是盲信号分离（blind source separation (BSS)）的一种特例。,一、问题的引出,在上述主成份分析中，我们可以注意到，把数据之间的相关性去除，并不足以使变量独立，特别是当变量是非高斯分布的时候。独立主成份分析是要寻找到一种变换把原始数据转换成若干个独立的变量。独立主成份分析的主要目的是要揭示数据当中一些更有意义的变量而不是减少数据的维数。,独立成分分析的最重要的假设就是信号源统计独立。这个假设在大多数盲信号分离的情况中符合实际情况。即使当该假设不满足时，仍然可以用独立成分分析来

16、把观察信号统计独立化，从而进一步分析数据的特性。独立成分分析的经典问题是“鸡尾酒会问题”（cocktail party problem）。该问题描述的是给定混合信号，如何分离出鸡尾酒会中同时说话的每个人的独立信号。当有N个信号源时，通常假设观察信号也有N个（例如N个MiC或者录音机）。该假设意味着混合矩阵是个方阵，即J = D，其中D是输入数据的维数，J是系统模型的维数。对于J D的情况，也分别有不同研究。,这种问题非常类似于EEG信号分析。EEG信号由放置在头部的许多电极记录产生，而这些信号实际上是隐含的神经信号源组成生产的。,mixed,ICA Separation,ICA和PCA信号分析

17、的不同在于PCA仅使用的是信号的二阶统计量，而ICA使用的是更高阶的统计量。易知，高斯分布的信号二阶以上的统计距为0，而很多信号是非高斯分布的，这样就会有更高阶的统计距存在。而ICA正好很好的利用了这些高阶统计特性。,二、ICA的系统模型,ICA通常假设被测信号是由一组独立信号源瞬时线性组成产生的。,注意有N个变量，就有N个方程。通常我们考虑信号的顺序与时间无关:s 与 x 是无关的时间函数。,s是信号源矢量，A是混合矩阵，x是被测信号矢量。,这一模型也称为隐变量模型。,ICA就是要求解混合矩阵A，使得,但很多时候我们不能准确获得混合矩阵A，这时我们通常采用估计的策略。,此外，可以看出使用这种

18、ICA模型并不能完全恢复信号源的具体数值，也不能解出信号源的正负符号、信号的阶数或者信号的数值范围。,三、ICA的计算步骤,1.中心化数据 2.白化（whiten)数据；去相关，各方向方差缩放为1。,上面两步可通过PCA的方法完成,3.获得原信号独立成份的估计,b应是一个合适的矢量估计能够重建独立成份。为估计b，我们要建立与独立变量有关的目标函数，然后通过最优化算法来估计出b。,一个直观的方法：利用中心极限定理：大量的独立分布的信号的混合最后会趋向于成为高斯分布的信号。,Gauss distribution （rand）,Single sinusoid at 100HZ,Two sinusoi

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

20、stICA，是一种交互式的程序。而我们使用Jade算法用在下面的例子。,%Example of ICA %create a mixture using three different signals mixed five ways plus noise clear all, close all; N=1000; fs=500; w=(1:N)*2*pi/fs; t=1:N; %generate the three signals plus noise s1=.75*sin(w*12)+.1*randn(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,