第二讲-EEG信号处理基础PPT课件

.,1,第二讲EEG信号处理基础,主讲人：谢宏信息工程学院,.,2,EEG信号模型EEG信号特征,.,3,EEG信号的特点,随机信号非平稳非线性非高斯过程,.,4,EEG信号模型,基于神经元的生化物理模型HodgkinandHuxley模型MorrisLecar模型基于随机信号的动态模型线性模型：AR模型和ARMA模型非线性模型：GARCH模型,.,5,EEG信号线性模型,AR模型,ARMA模型,AIC准则确定模型阶数：即最小化如下目标函数,多维AR模型：对应多通道,.,6,诱发电位的Prony方法,首先建立输入为脉冲函数的AR模型：,假设有样本值y(1),y(2),y(N),令,.,7,解得AR模型系数：,求解特征方程：,得到其p个根,方程的通解为：,.,8,求解,其中,.,9,EEG信号非线性模型,AR模型,ARMA模型,AIC准则确定模型阶数：即最小化如下目标函数,多维AR模型：对应多通道,.,10,由模型化神经元活动产生EEG信号,AR模型,ARMA模型,.,11,由模型化神经元活动产生EEG信号,AR模型,.,12,EEG信号常用统计特征,随机变量或向量统计量均值、方差、协方差与相关系数矩阵偏度(skewness)：峭度或峰度(kurtosis)：随机过程或随机信号统计量自相关函数与功率谱双谱等高阶谱,.,13,EEG信号常用统计特征,设m个电极通道采集的脑电信号表示为,均值,其N点的采样数据为,.,14,EEG信号常用统计特征,方差,偏度,峭度(峰度),协方差矩阵,.,15,EEG信号常用统计特征,相关系数矩阵,注意：以上统计量只涉及随机变量取值，不涉及时序计算时要考虑实时性，如均值的计算，可以考虑采用算法,.,16,平稳EEG信号时序特征,自相关函数,功率谱：信号的功率谱是其自相关函数的傅里叶变换，即,多通道信号的自相关函数矩阵,.,17,平稳EEG信号时序特征,功率谱的计算：周期图法与线性时序模型法周期图法计算频谱得到，简单、功率谱不光滑、误差大；假设EEG信号的时序模型为AR(p)，即,其功率谱为,.,18,.,19,平稳EEG信号时序特征,EEG信号节律波能量提取：设信号x(n)的傅里叶变换为X()，则有,对有限个采样值x(0),x(1),x(N-1),其DFT为X(k),对应频率kfs/N,则有：,则EEG信号节律波平均能量提取算法可以考虑两种方案。,.,20,平稳EEG信号时序特征,设节律波频带范围为fL,fH，脑电信号的N个采样值x(0),x(1),x(N-1)方案一：第一步：采用FFT算法计算DFTX(0),X(1),X(N-1);第二步：计算NL=fL*N/fs,NH=fH*N/fs第三步：,方案二第一步：对脑电信号x(n)进行带通fL,fH滤波得到输出信号y(n);第二步：,.,21,EEG信号滤波,信号滤波涉及：低通、高通、带通、陷波滤波器的比较：IIR满足相同特性阶数较低，只能近似线性相位，必须浮点运算FIR满足相同特性阶数较高，可以做到严格线性相位，可以采用整数运算滤波器一般要结合实现时的计算效率和滤波器特性等综合考虑阶数通带、阻带和过渡带特性延迟,.,22,IIR滤波器,IIR滤波器的模型为：IIR滤波器的类型：贝塞尔、巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆型相位特性：贝塞尔巴特沃斯切比雪夫椭圆过渡带宽度：贝塞尔巴特沃斯0时，系统具有混沌特征。,.,69,时间序列Lyapunov指数的计算,在实际时间序列混沌识别中，通常只估计最大Lyapunov指数，下面介绍一种算法:小数据量法。设时间序列x1,x2,xN,嵌入维数m，时间延迟，重构相空间后有：其中kN-(m-1),选取一距离阈值，对初始点x(0)，选取与其距离最近的点不妨设为x(k0)=z0(0)，当经过t0次迭代后，z0(t0)与x(t0)的距离大于，令,.,70,时间序列Lyapunov指数的计算,取与x(t0)的距离小于的点x(k1)=z1(0)，且z1(0)-x(t0)与z0(t0)-x(t0)的夹角最小，如图所示重复以上过程，假设直到M-1步结束，此时tM-1=N-(m-1)，则有,.,71,Kolmogorov熵,Kolmogorov熵K在混沌的度量中有着相当重要的应用。对于规则运动，K=0;随机系统K=无穷，若系统表现确定性混沌，则Kolmogorov熵是大于0的常数。Kolmogorov熵越大，那么信息的损失速度越大，系统的混沌程度越大。,.,72,混沌相空间重构理论,相空间重构是对一维时间序列，按照延迟时间和嵌入维数重构一个与原动力系统等价的相空间;相空间重构的理论基础是Takens相空间嵌入定理，即通过对一个m维流形上的连续流的观测值h(t)的延迟(h(t),h(t+),h(t+2m)，可以将该连续流光滑嵌入到2m+1维空间(相空间)中；相空间重构过程中有两个参数选取特别重要:延迟时间和嵌入维数.,.,73,延迟时间间隔的计算,主要方法线性自相关函数法平均互信息法,.,74,延迟时间间隔的选取线性自相关函数法,定义自相关函数为,选择使得自相关函数C()第一次为零时的的值。缺点：对嵌入维数大于2时不是最优的。,.,75,延迟时间间隔的选取平均互信息法,选择使I()为第一个局部极小的为延迟时间间隔。缺点：要对概率密度进行估计，需要的样本点多，计算复杂，且对嵌入维数大于2维不是最优的。,.,76,嵌入维数m的计算,主要方法Lyapunov指数法关联积分法,.,77,嵌入维数m的计算,Lyapunov指数法（1）首先确定延迟；（2）依次递增取嵌入维数m，计算嵌入相空间(x(n),x(n-),x(n-(m-1)的Lyapunov指数，选取指数趋于常数的m作为嵌入相空间维数。缺点：当样本点数不多时可信度较低，当嵌入维数较高时面临“维数灾”。,.,78,嵌入维数m的计算,关联积分:定义：设t为时间序列的延时，m是嵌入维数，N是样本数据集的大小，数据点个数M=N-(m-1)t，则重构相空间中嵌入时间序列Y(i)的关联积分定义为：关联积分是一个累积分布函数，表示相空间中任意两点之间距离小于半径r的概率，这里点与点之间的距离用矢量之差的无穷范数表示。,.,79,嵌入维数m的计算,关联积分法的主要步骤(1)利用时间序列X1,X2,XN,先给定一个较小的嵌入维数m0,重构相空间，得到新的序列Yi(2)计算关联积分C(r)(3)对于r的某个取值范围，吸引子的维数d与累积分布函数C(r)应满足对数线性关系，即d(m)=LnC(r)/Lnr,从而可用最小二乘拟合得到对应于m0的关联维数估计d(m0)(4)增加嵌入维数m0,重新计算步骤(2)和(3),直到相应的维数估计值d(m)不再随着m的增加而在一定误差范围内不变为止。缺点：当样本点数不多时可信度较低，当嵌入维数较高时面临“维数灾”。,.,80,.,81,总结,关于混沌判定，一般应用最大Lyapunov指数，或者Kolmogorov熵，或者结合两者判定。在计算延迟时间方法上，常用的有自相关函数法和互信息法，计算嵌入维常用的方法有G-P算法，这种方法是先计算出混沌时间序列的关联维，然后再计算出嵌入维数。,.,82,自适应滤波噪声消除,自适应滤波是一种采用参数可在线调整的有限脉冲响应滤波器进行滤波；相对于固定参数滤波器，自适应滤波器有可能达到更好的去除噪声的效果；相对于主成分和独立分量分析，自适应滤波算法简单；自适应滤波需要有欲滤除噪声的参考信号。,.,83,.,84,自适应滤波算法一,可以考虑由每一次的采样值修正参数w，由,为减小该式误差，可以按照负梯度方向修正w，即,考虑到稳定性，一般应比较小，满足,.,85,自适应滤波算法二,假设有n个样本，考虑遗忘因子加权的误差,令其关于w的梯度为零得：,得到,.,86,自适应滤波算法二,令,有,.,87,自适应滤波总结,自适应滤波需要确定参考信号；需要确定参考信号的FIR滤波器阶数M；计算复杂度低；有针对单样本的在线算法；在EEG信号处理中的应用，如：对滤除眨眼的眼电干扰信号，参考信号可以考虑从FP1和FP2上取；ERP中，参考信号可以考虑各段的平均值；在去除心电干扰中，可以取心电信号作为参考信号。,.,88,主成分分析(K-L变换),信号的DFT和DCT是将信号投影或表示到固定正交函数系上的变换主成分分析或称K-L变换是将信号表示到一组标准正交信号系上的变换，该信号系一般由信号的统计特性得到该正交信号系张成的空间称为信号子空间，其正交补空间称为噪声子空间主成分分析或K-L变换的应用用于分离信号与噪声对信号的维数进行压缩降维,.,89,主成分的定义及导出,设x(n)=(x1(n),x2(n),xp(n)T为一个p维平稳随机信号，E(x(n)=0，其协方差矩阵为该矩阵为实对称矩阵，其非负特征值设为，则存在正交矩阵U，使得,.,90,令：,则有：，因此y(n)的任意两个分量不相关。y(n)的分量称为X的主分量。由于,总方差(平均能量)中属于第i主成分yi的比例为称为主成分yi的贡献率。,.,91,主成分的计算,对原始信号数据,由其奇异值分解，知存在p阶正交矩阵U和N阶正交矩阵V，使得其中,.,92,样本协方差矩阵,因此，矩阵R的特征值为,主分量为,V的前p行为正交信号，由其张成的空间为信号子空间,.,93,主成分分析在降维中的应用,首先对样本协方差矩阵R计算正交变换矩阵U和特征值计算前m个主成分的贡献率之和对给定的阈值(一般在0.80.95),选取累计贡献达到最小的m。由y(n)=Umx(n)(Um为矩阵U的前m行子矩阵)计算得到m维向量，可用来代替p维向量x(n)，从而达到降维的目的，而信息的损失却不多。,.,94,主成分分析在降噪中的应用,首先对样本数据矩阵X计算正交矩阵U和V对实际采样信号xi(n)计算分解,也可以采用在线最小二乘算法(LMS)计算分解，即由,.,95,主成分分析总结,要求信号服从高斯分布各个主分量之间不相关，但不保证独立,.,96,独立分量分析,独立分量分析的目的是：当X=AS时，求矩阵W，使得Y=WX的各个分量独立，此时W可能不是A的逆，但是WA是置换矩阵。典型例子：“鸡尾酒会”问题，从酒会嘈杂人声中提取所关心对象的语音,人的大脑可以很快辨出或集中听某种需要关注声音。,.,97,独立分量分析,主要研究结构:(1)美国加州大学生物系，计算神经生物实验室，提出信息极大化(infomax)算法。.(2)日本Riken的数量神经科学实验室，互信息极小化（minimizationofmutualinformationMMI)采用人工神经网络优化。http:/www.brain.riken.jp/lab/mns/amari(3)芬兰赫尔辛基工业大学神经网络研究中心，提出了立足于逐次提取独立分量的固定点算法（fixedpointalgorithm)：fastICA。www.cis.hut.fi/ojawww.cs.helsinki.fi/aapo.hyvarinen(4)法国学者：J.F.Cardoso，提出了JADE算法、批数据处理算法、近年来引人注意的稀疏分量分析。http:/tsi.enst.fr/cardoso.,.,98,ICA相关的基本概念,n阶矩(moment)：mn=E(xn)特征函数第二特征函数n阶累计量(cumulant),对高斯型信号，二阶以上的累计量都为0，因此可由一、二阶统计特征来完整描述。k40超高斯,k4loadfoetal_ecg.datasciiX=foetal_ecg(:,2:9);plot(X(1:1000,1);,.,121,.,122,ICA分解Y=fastica(X);plot(Y(1:1000,1);8个分量波形,.,123,.,124,PCA分解m,n=size(X);X=X-ones(m,1)*mean(X);A=X*X/m;P,D=eig(A);d=diag(D);d=d(n:-1:1);P=P(:,n:-1:1);Y=X*P;,.,125,.,126,CommonSpatialPattern(CSP)