EEG信号处理基础和特征讲解

1、EEG信号处理基础和特征讲解EEG信号模型EEG信号特征EEG信号的特点随机信号非平稳非线性非高斯过程EEG信号模型基于神经元的生化物理模型Hodgkin and Huxley 模型MorrisLecar 模型基于随机信号的动态模型线性模型：AR模型和ARMA模型非线性模型：GARCH模型EEG信号线性模型AR模型ARMA模型AIC准则确定模型阶数：即最小化如下目标函数多维AR模型：对应多通道诱发电位的Prony方法首先建立输入为脉冲函数的AR模型：假设有样本值y(1),y(2),y(N), 令解得AR模型系数：求解特征方程： 得到其p个根方程的通解为：求解其中EEG信号非线性模型AR模型AR

2、MA模型AIC准则确定模型阶数：即最小化如下目标函数多维AR模型：对应多通道由模型化神经元活动产生EEG信号AR模型ARMA模型由模型化神经元活动产生EEG信号AR模型EEG信号常用统计特征随机变量或向量统计量均值、方差、协方差与相关系数矩阵偏度(skewness)：峭度或峰度(kurtosis)：随机过程或随机信号统计量自相关函数与功率谱双谱等高阶谱EEG信号常用统计特征设m个电极通道采集的脑电信号表示为均值其N点的采样数据为EEG信号常用统计特征方差偏度峭度(峰度)协方差矩阵EEG信号常用统计特征相关系数矩阵注意：以上统计量只涉及随机变量取值，不涉及时序计算时要考虑实时性，如均值的计算，可

3、以考虑采用算法平稳EEG信号时序特征自相关函数功率谱：信号的功率谱是其自相关函数的傅里叶变换，即多通道信号的自相关函数矩阵平稳EEG信号时序特征功率谱的计算：周期图法与线性时序模型法周期图法计算频谱得到，简单、功率谱不光滑、误差大；假设EEG信号的时序模型为AR(p)，即其功率谱为平稳EEG信号时序特征EEG信号节律波能量提取： 设信号x(n)的傅里叶变换为X()，则有对有限个采样值x(0),x(1),x(N-1), 其DFT为X(k), 对应频率kfs/N, 则有：则EEG信号节律波平均能量提取算法可以考虑两种方案。平稳EEG信号时序特征设节律波频带范围为fL,fH，脑电信号的N个采样值x(

4、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);第二步：EEG信号滤波信号滤波涉及：低通、高通、带通、陷波滤波器的比较：IIR满足相同特性阶数较低，只能近似线性相位，必须浮点运算FIR满足相同特性阶数较高，可以做到严格线性相位，可以采用整数运算滤波器一般要结合实现时的计算效率和滤波器特性等综合考虑阶数通带、阻带和过渡带特性延迟IIR滤波器IIR滤波器的模型为：IIR滤波器的类型：贝塞尔、巴特沃

5、斯、切比雪夫I型、切比雪夫II型和椭圆型相位特性：贝塞尔巴特沃斯切比雪夫椭圆过渡带宽度：贝塞尔巴特沃斯切比雪夫0发散，a=0为常数，a0 时，系统具有混沌特征。时间序列Lyapunov指数的计算 在实际时间序列混沌识别中，通常只估计最大Lyapunov指数，下面介绍一种算法: 小数据量法。 设时间序列x1,x2,xN, 嵌入维数m，时间延迟，重构相空间后有： 其中kN-(m-1) ,选取一距离阈值，对初始点x(0)，选取与其距离最近的点不妨设为x(k0)=z0(0)，当经过t0次迭代后，z0(t0)与x(t0)的距离大于，令时间序列Lyapunov指数的计算 取与x(t0)的距离小于的点x(k

6、1)=z1(0)，且z1(0)-x(t0)与z0(t0)-x(t0)的夹角最小，如图所示重复以上过程，假设直到M-1步结束，此时tM-1=N-(m-1)，则有 Kolmogorov熵 Kolmogorov熵K在混沌的度量中有着相当重要的应用。对于规则运动，K=0; 随机系统K=无穷，若系统表现确定性混沌 ， 则 Kolmogorov熵是大于0的常数。Kolmogorov熵越大，那么信息的损失速度越大，系统的混沌程度越大。混沌相空间重构理论相空间重构是对一维时间序列，按照延迟时间和嵌入维数重构一个与原动力系统等价的相空间;相空间重构的理论基础是Takens相空间嵌入定理，即通过对一个m维流形上的

7、连续流的观测值h(t)的延迟(h(t),h(t+ ),h(t+2m)，可以将该连续流光滑嵌入到2m+1维空间(相空间)中；相空间重构过程中有两个参数选取特别重要: 延迟时间 和嵌入维数.延迟时间间隔 的计算主要方法 线性自相关函数法平均互信息法延迟时间间隔 的选取线性自相关函数法定义自相关函数为 选择使得自相关函数C()第一次为零时的 的值。缺点：对嵌入维数大于2时不是最优的。延迟时间间隔 的选取平均互信息法选择使I( ) 为第一个局部极小的为延迟时间间隔。缺点：要对概率密度进行估计，需要的样本点多，计算复杂，且对嵌入维数大于2维不是最优的。嵌入维数m的计算主要方法Lyapunov指数法关联积

8、分法嵌入维数m的计算Lyapunov指数法（1）首先确定延迟 ；（2）依次递增取嵌入维数m，计算嵌入相空间(x(n),x(n- ),x(n-(m-1)的Lyapunov指数，选取指数趋于常数的m作为嵌入相空间维数。缺点：当样本点数不多时可信度较低，当嵌入维数较高时面临“维数灾”。嵌入维数m的计算关联积分: 定义：设t 为时间序列的延时，m是嵌入维数，N是样本数据集的大小，数据点个数M=N-(m-1)t ，则重构相空间中嵌入时间序列Y(i)的关联积分定义为： 关联积分是一个累积分布函数，表示相空间中任意两点之间距离小于半径r的概率，这里点与点之间的距离用矢量之差的无穷范数表示。嵌入维数m的计算关

9、联积分法的主要步骤(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的增加而在一定误差范围内不变为止。缺点：当样本点数不多时可信度较低，当嵌入维数较高时面临“维数灾”。总 结 关于混沌判定，一般应用最大Lyapunov指数，或者Kolmogorov熵，或者结合两者

10、判定。在计算延迟时间方法上，常用的有自相关函数法和互信息法，计算嵌入维常用的方法有G-P算法，这种方法是先计算出混沌时间序列的关联维，然后再计算出嵌入维数。自适应滤波噪声消除自适应滤波是一种采用参数可在线调整的有限脉冲响应滤波器进行滤波；相对于固定参数滤波器，自适应滤波器有可能达到更好的去除噪声的效果；相对于主成分和独立分量分析，自适应滤波算法简单；自适应滤波需要有欲滤除噪声的参考信号。自适应滤波算法一可以考虑由每一次的采样值修正参数w，由为减小该式误差，可以按照负梯度方向修正w，即考虑到稳定性，一般应比较小，满足自适应滤波算法二假设有n个样本，考虑遗忘因子加权的误差令其关于w的梯度为零得：得

11、到自适应滤波算法二令有自适应滤波总结自适应滤波需要确定参考信号；需要确定参考信号的FIR滤波器阶数M；计算复杂度低；有针对单样本的在线算法；在EEG信号处理中的应用，如：对滤除眨眼的眼电干扰信号，参考信号可以考虑从FP1和FP2上取；ERP中，参考信号可以考虑各段的平均值；在去除心电干扰中，可以取心电信号作为参考信号。主成分分析(K-L变换)信号的DFT和DCT是将信号投影或表示到固定正交函数系上的变换主成分分析或称K-L变换是将信号表示到一组标准正交信号系上的变换，该信号系一般由信号的统计特性得到该正交信号系张成的空间称为信号子空间，其正交补空间称为噪声子空间主成分分析或K-L变换的应用用于

12、分离信号与噪声对信号的维数进行压缩降维主成分的定义及导出设 x(n)=(x1(n),x2(n),xp(n)T 为一个 p 维平稳随机信号，E(x(n)=0，其协方差矩阵为 该矩阵为实对称矩阵，其非负特征值设为 ，则存在正交矩阵U，使得令：则有： ， ，因此 y(n) 的任意两个分量不相关。y(n)的分量称为X的主分量。由于总方差(平均能量)中属于第 i 主成分 yi 的比例为 称为主成分 yi 的贡献率。主成分的计算对原始信号数据由其奇异值分解，知存在 p 阶正交矩阵U和 N 阶正交矩阵V，使得 其中样本协方差矩阵因此，矩阵R的特征值为主分量为V的前p行为正交信号，由其张成的空间为信号子空间主

13、成分分析在降维中的应用首先对样本协方差矩阵R计算正交变换矩阵U和特征值计算前 m 个主成分的贡献率之和对给定的阈值(一般在0.80.95), 选取累计贡献达到最小的m。由 y(n)=Umx(n) (Um为矩阵U的前m行子矩阵)计算得到m维向量，可用来代替p维向量x(n)，从而达到降维的目的，而信息的损失却不多。主成分分析在降噪中的应用首先对样本数据矩阵X计算正交矩阵U和V对实际采样信号xi(n)计算分解也可以采用在线最小二乘算法(LMS)计算分解，即由主成分分析总结要求信号服从高斯分布各个主分量之间不相关，但不保证独立独立分量分析 独立分量分析的目的是：当 X = A S 时，求矩阵 W，使得

14、 Y=WX 的各个分量独立，此时W可能不是A的逆，但是 WA 是置换矩阵。典型例子： “鸡尾酒会”问题，从酒会嘈杂人声中提取所关心对象的语音, 人的大脑可以很快辨出或集中听某种需要关注声音。麦克风1麦克风2麦克风3独立分量分析主要研究结构:(1) 美国加州大学生物系，计算神经生物实验室，提出信息极大化( infomax )算法。 .(2) 日本Riken的数量神经科学实验室，互信息极小化（minimization of mutual information MMI)采用人工神经网络优化。 (3) 芬兰赫尔辛基工业大学神经网络研究中心，提出了立足于逐次提取独立分量的固定点算 法（fixed po

15、int algorithm) ：fastICA。 (4) 法国学者：，提出了JADE算法、批数据处理算法、近年来引人注意的稀疏分量分析。 . ICA相关的基本概念n阶矩(moment)：mn=E(xn) 特征函数第二特征函数n阶累计量(cumulant) 对高斯型信号，二阶以上的累计量都为0，因此可由一、二阶统计特征来完整描述。k40 超高斯, k40 亚高斯, 用k4大小作为衡量信号距离高斯型程度的度量。均值 方差 峭度ICA相关的基本概念联合矩中最常用的是协方差；联合累计量最常用的是4阶累积量，一般用cum(x1,x2,x3,x4)表示联合累积量性质：当x各分量相互独立时，其互累计量恒等于

16、0。比例性：cum(w1x1,w2x2,w3x3,w4x4)=w1w2w3w4cum(x1,x2,x3,x4)当两个随机变量x, y独立时，有k4(x+y)=k4(x)+k4(y)随机向量x=(x1,x2,xn)各个分量独立当且仅当其联合概率分布(密度)等于各个分量概率分布(密度)的乘积。ICA计算数据的预处理解混系统B白化W正交系统UZ(t)混合系统Ax(t)y(t)s(t)系统简图 一般在ICA求解之前，先对数据进行白化处理，这样可以使得各个分量之间是不相关的，便于进行解混计算，如下图所示样本数据的白化对原始信号数据假设数据已经零均值化，先计算数据的协方差矩阵计算其特征值和特征向量，得到对

17、角化分解：Rx=UTU计算白化矩阵W： ，计算白化数据：Z=WX， 此时z的协方差矩阵是单位阵，即：Rz=I。fastICA-基于非高斯极大化的算法假设数据经过白化处理，即m维随机向量z的均值向量为零向量，协方差矩阵为单位阵；采用一次提取一个分量的方法，以极大化非高斯性作为优化的目标函数，可以考虑以下两种度量极大化峭度(Kurtosis)的绝度值极小化负熵(Negentropy)由于峭度：Kurt(y)=Ey4-Ey22 的鲁棒性不好，一般 fastICA 算法采用负熵作为目标函数。fastICA算法也称为固定点算法或投影寻踪算法。牛顿迭代法牛顿法迭代法是用于求解方程 f(x)=0 的解。设在

18、某一点Pk，对应x坐标为xk，其切线方程为与x轴的交点为xk+1，则有迭代公式例.用牛顿迭代法求方程的根:解:由牛顿迭代法x0 =0.5;迭代四次精度达10-8 负熵信息熵： 当随机变量的均值和方差一定时，其 pdf 为高斯分布时其信息熵最大。负熵：任意 pdf 的 p(x) 和具有相同协方差矩阵的高斯分布pG(x)的差作为该 pdf 非高斯程度的度量，定义为负熵 J(x)0，p(x)=pG(x)：当且仅当Jp(x)=0代表高斯分布。 由以上公式可知，计算负熵需要随机变量的概率密度，在实际计算中是不方便的，因此需要给出其计算的近似公式。负熵的近似高阶统计量形式：设x零均值，方差为1(白化数据)

19、Edgeworth展开Gram-Charlier展开缺点：大值野点会引起较大误差负熵的近似非多项式函数形式： 其中G(1)(y)为奇函数，G(2)(y)为偶函数，v为与y有相同均值与方差的高斯随机变量。可以考虑如下更简洁的形式fastICA算法就以上式为基础。设y=wTz，其关于w的梯度为所以其极值点使 gradw(EG(y)=EG(wTz)z=0 基于负熵的牛顿迭代公式考虑到约束wTw=1，由拉格朗日乘子法得极值点满足： 考虑映射： 由线性近似：得到迭代公式：基于负熵的牛顿迭代公式于是考虑到单位化，可以简化为由于G(y)为偶函数，其导数g(y)为奇函数，一般可选 基于负熵的fastICA基于

20、负熵提取一个分量的固定点算法（2001）（1）X去均值，进行白化得Z；（2）任选 w1 的初值，使 w1的2-范数为1，k=1；（3）令 ，均值用样本均值；（4）归一化：wk+1=wk+1/|wk+1|;（5）若 不接近1，返回步骤（3），否则迭代结束（6）提取分量该方法收敛具有三阶收敛速度基于负熵的fastICA基于负熵逐个提取多个分量的固定点算法（1）X去均值，进行白化得Z；（2）m为待提取独立分量数目，p=1；（3）任选 wp,1 的初值，使 wp,1 的2-范数为1，k=1；（4）令 ，均值用样本均值；（5）正交化：（6）归一化：wp,k+1=wp,k+1/|wp,k+1|;（7）若 不接近1，k=k+1,返回步骤（4）;（8）若p0，取gi(yi)=-2tanh(yi)，否则g(yi)=tanh(yi)-yi（5）更新分离矩阵BB+I+g(y)yTB（6）如果尚未收敛