第二讲EEG信号处理基础.ppt

第二讲 EEG信号处理基础主讲人：谢宏主讲人：谢宏信息工程学院信息工程学院vEEG信号模型vEEG信号特征EEG信号的特点v随机信号v非平稳v非线性v非高斯过程EEG信号模型v基于神经元的生化物理模型™Hodgkin and Huxley 模型™Morris–Lecar 模型v基于随机信号的动态模型™线性模型：AR模型和ARMA模型™非线性模型：GARCH模型EEG信号线性模型vAR模型vARMA模型vAIC准则确定模型阶数：即最小化如下目标函数v多维AR模型：对应多通道诱发电位的Prony方法方法v首先建立输入为脉冲函数的AR模型：v假设有样本值y(1),y(2),…,y(N), 令v解得AR模型系数：v求解特征方程： 得到其p个根v方程的通解为：v求解其中EEG信号非线性模型vAR模型vARMA模型vAIC准则确定模型阶数：即最小化如下目标函数v多维AR模型：对应多通道由模型化神经元活动产生EEG信号vAR模型vARMA模型由模型化神经元活动产生EEG信号vAR模型EEG信号常用统计特征v随机变量或向量统计量™均值、方差、协方差与相关系数矩阵™偏度(skewness)：™峭度或峰度(kurtosis)：v随机过程或随机信号统计量™自相关函数与功率谱™双谱等高阶谱EEG信号常用统计特征v设m个电极通道采集的脑电信号表示为v均值v其N点的采样数据为EEG信号常用统计特征v方差v偏度v峭度(峰度)v协方差矩阵EEG信号常用统计特征v相关系数矩阵v注意：™以上统计量只涉及随机变量取值，不涉及时序™计算时要考虑实时性，如均值的计算，可以考虑采用算法平稳EEG信号时序特征v自相关函数v功率谱：信号的功率谱是其自相关函数的傅里叶变换，即v多通道信号的自相关函数矩阵平稳EEG信号时序特征v功率谱的计算：周期图法与线性时序模型法v周期图法计算频谱得到，简单、功率谱不光滑、误差大；v假设EEG信号的时序模型为AR(p)，即v其功率谱为平稳EEG信号时序特征vEEG信号节律波能量提取： 设信号x(n)的傅里叶变换为X(ω)，则有对有限个采样值x(0),x(1),…,x(N-1), 其DFT为X(k), 对应频率kfs/N, 则有：则EEG信号节律波平均能量提取算法可以考虑两种方案。

平稳EEG信号时序特征设节律波频带范围为[fL,fH]，脑电信号的N个采样值x(0), x(1), … , x(N-1)v方案一：v第一步：采用FFT算法计算DFTX(0),X(1),…,X(N-1); v第二步：计算NL=fL*N/fs, NH=fH*N/fsv第三步：v方案二v第一步：对脑电信号x(n)进行带通[fL,fH]滤波得到输出信号y(n);v第二步：EEG信号滤波v信号滤波涉及：低通、高通、带通、陷波v滤波器的比较：™IIR满足相同特性阶数较低，只能近似线性相位，必须浮点运算™FIR满足相同特性阶数较高，可以做到严格线性相位，可以采用整数运算v滤波器一般要结合实现时的计算效率和滤波器特性等综合考虑™阶数™通带、阻带和过渡带特性™延迟IIR滤波器vIIR滤波器的模型为：vIIR滤波器的类型：贝塞尔、巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆型v相位特性：贝塞尔>巴特沃斯>切比雪夫>椭圆v过渡带宽度：贝塞尔<巴特沃斯<切比雪夫<椭圆EEG信号α波提取滤波器设计v频带为8-13Hz，所以设计带通滤波器，采样频率为512，选择4阶椭圆滤波器，通带边缘频率为[8.1, 12.8][B,A]=ellip(4,0.5,20,[8.1/256,12.8/256])[h,f]=freqz(B,A,1000,512);plot(f(1:100),20*log(abs(h(1:100))))axis([6,15,-60,5])grid on滤波器参数脑电信号α波提取频谱比较非平稳信号分析v短时傅里叶变换v小波变换vWigner–Ville 分布2024/9/17连续小波变换说明va为尺度因子(对应频率)，较小的a对应高频，较大的a对应低频；vb为位移因子(对应时间)；vψ(t)为小波母函数，一般取具有单位能量的窗函数；v小波变换的值表示了信号f(t)与小波函数匹配的程度，例如若对某个a和b的取值信号f(t)与小波完全相同，则小波变换为1。

连续小波变换连续小波变换小波反变换小波反变换小波波形随尺度因子和位移因子的变化a a=1=1, , b b=0=0a a=1=1, , b b=6=6a a=3=3, , b b=6=6a a=1/3=1/3, b, b=6=6常用连续小波函数1 1、、Morlet Morlet 小波小波2 2、墨西哥草帽小波、墨西哥草帽小波 3 3、、DOGDOG小波小波 函函数数图图像像函函数数图图像像函函数数图图像像幅幅频频谱谱幅幅频频谱谱幅幅频频谱谱突变信号的墨西哥草帽变换vx=[sin(pi*t/200) sin(pi*t/30)];vC=cwt(x, 1:32, ‘mexh’, ‘plot’);离散正交小波变换v离散小波变换相当于在连续小波变换公式中取 a = 2-j/2，b = k 2-j/2；v小波反变换成为小波级数；v小波函数必须满足一定的条件，才能使以上变换公式和反变换公式成立；v离散小波变换需要计算积分，不利于实际应用，需要更高效的算法离散小波变换离散小波变换小波反变换小波反变换多分辨率分析和金子塔算法小波函数由尺度函数确定，尺度函数一般满足小波函数则可表示为设信号f(t)可以表示为分解算法 重构算法 注意：此算法与尺度函数和小波函数的形式无关令：qk = (-1)k-1p-k+1分解算法和重构算法的含意1、金字塔算法是对信号按频带逐层分解，一直达到需要的频带为止；2、尺度函数分量为低通分量，小波分量为各个频带的带通分量。

记：则有：fM(t)=fM-1(t)+gM-1(t) =fM-2(t)+gM-2(t)+gM-1(t)=… =fM-N(t)+gM-N(t)+…+gM-1(t)HAAR正交小波尺度系数：p0=1, p1=1小波系数：q0=1, q1= -1特点：1、非零尺度系数和小波系数个数有限；2、尺度函数和小波函数的非零区域为[0,1](紧支撑)；3、尺度函数和小波函数不连续，频率窗太宽Db2紧支撑正交小波尺度系数：p0=0.4829629131445341, p1p2p4尺度函数小波函数特点：1、非零尺度系数和小波系数个数有限； 2、尺度函数和小波函数的非零区域为[-4,4](紧支撑)； 3、尺度函数和小波函数连续常用小波vHaar小波(可以看作为Daubechies小波的特例)vDaubechies正交紧支撑小波(波形不对称)v半正交小波(波形具有对称性)v紧支撑双正交小波(波形可以具有对称性)采用Db3对sin函数的和构成信号的分解采用Db2对频率突变信号的分解采用Db5对频率突变信号的分解采用Db9对频率突变信号的分解小波分解与重构法去除基线漂移原脑电信号加入基线漂移后的脑电信号去除缓慢基线漂移后的脑电信号脑电信号的7层分解 多通道同步性和因果性v为了描述通道之间EEG信号的同步性或协同性，相关系数、相干性(coherence)等是常用的工具，由于相关系数会受信号延迟的影响，因此相干性更有优势。

v相干性定义：两个通道的相干性定义为 Xi(ω)和Xj(ω)分别是相应通道EEG信号的傅里叶变换GrangerGranger因果性因果性v定定义：对两个时间序列x(n)和y(n)，如果对x(n)、y(n)过去的了解比仅了解y(n)的过去减少了y(n)预测误差的方差，则称x(n)引起了y(n)v假设m个通道的EEG信号具有如下MVAR模型： 其中Lk是mⅹm阶方阵，对其两边取Fourier变换GrangerGranger因果性因果性v即￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿为MVAR模型的系统(传输)函数 若则称分量j与分量i有因果关系此外，也有考虑非线性模型的因果性的讨论动态系统与混沌分析v混沌的起源和发展 牛顿经典力学表明，力学系统服从确定的规律即当初始条件确定后，力学系统就将按确定的轨道运动在我们日常经验中，事物的一般重复现象很明显这些现象似乎表明我们的世界是规则、和谐、有序的，自然现象的变化是周期的、重复的 1963年美国气象学家洛伦兹发现大气变化的非周期性，天气预报难就难在天气变化不是周期性的动态系统与混沌分析v混沌的起源和发展 混沌理论把混沌描述为无序的，非周期性的现象。

自从科学家不懈地探索自然规律以来，无序、非周期性现象就一直被忽略自然界中不规则的方面、不连续和不稳定的方面，一直是科学的难题，是无法理解的怪物科学家常常把它们当做噪声或科学垃圾扔掉而混沌不同于噪声，有其自身的规律这正是科学家所要研究的.动态系统与混沌分析v连续动态系统：常微分方程组描述v离散动态系统：差分方程组描述v线性动态系统：当初值确定后，解就确定了，具有性质：™一般具有三种形态：收敛到一点(0维吸引子)、周期震荡、发散™当初值相差较小时，其轨道之间的差别也较小™具有确定性信号的特征线性连续动态系统v1维线性动态系统 其解为 ，a>0发散，a=0为常数，a<0收敛于0v2维线性动态系统 其轨道为x(0)eAt：当其系数矩阵A的特征值，有实部大于零时发散；实部都小于零时收敛到(0,0)，实部都等于零时为椭圆线性离散动态系统v2维线性动态系统 其轨道为x(0)An：当A的特征值模大于1时发散；模都小于1时收敛到(0,0)；模等于1时其约当标准型为以下几种形式：周期为2发散发散周期震荡非线性动态系统中的混沌v与线性系统不同，非线性系统可能出现非常复杂的行为(轨道)，混沌现象就是之一，其主要特征为：™对初态的敏感依赖性™确定论系统中的内随机性™局部不稳定而整体稳定™奇怪吸引子：吸引子具有分形结构，因此被称为奇怪吸引子。

™连续功率谱：与噪声(如高斯白噪声)具有相似的功率谱非线性动态系统中的混沌v与线性系统不同，非线性系统可能出现非常复杂的行为(轨道)，混沌现象就是之一，其主要特征为：™对初态的敏感依赖性™确定论系统中的内随机性™局部不稳定而整体稳定™奇怪吸引子：吸引子具有分形结构，因此被称为奇怪吸引子™连续功率谱：与噪声(如高斯白噪声)具有相似的功率谱非线性动态系统中的混沌v几种混沌的定义™非线性确定性系统中, 由于系统内部非线性相互作用而产生的一种非周期的行为™对初始状态敏感，表现似周期、非周期和不可预报性的过程 ™在确定性的非线性动态系统中出现的貌似随机的、不能预测的运动 ™对初始条件敏感的非线性确定性系统，具有正的李雅普诺夫指数 58离散动态系统v考虑如下描述人口或虫口问题的Logistic模型，其状态演化方程：xn=axn-1(1－－xn-1)xn=axn-1(1－－xn-1)v参数参数a=2.5，，x(0)=0.3，，迭代迭代200次，状态收敛于次，状态收敛于0.6xn=axn-1(1－－xn-1)参数参数a=3.1，，x(0)=0.3，在，在0.558和和0.7646两个状态变化两个状态变化xn=axn-1(1－－xn-1)参数参数a=3.5，，x(0)=0.3，在四个状态变化，周期为，在四个状态变化，周期为4xn=axn-1(1－－xn-1)参数参数a=3.6，，x(0)=0.3，状态进入混沌状态，状态进入混沌状态63xn=axn-1(1－－xn-1)v模型分析结论当当a>3.6时进入混沌状态时进入混沌状态64正 文v分析论证过程65罗沦兹方程 其中其中x x、、y y、、z z为无量纲量，为无量纲量，分别表征对流强度，对流中升分别表征对流强度，对流中升流与降流间的温差和竖直方向流与降流间的温差和竖直方向温度分布的非线性度。

任意给温度分布的非线性度任意给定初值，系统最终都会回到状定初值，系统最终都会回到状态空间的特定区域内，其吸引态空间的特定区域内，其吸引子具有精巧而奇特的结构，如子具有精巧而奇特的结构，如图所示，表明系统进入了混沌图所示，表明系统进入了混沌状态 (图为图为LorenzLorenz混沌吸引混沌吸引子子) )混沌识别 混沌识别主要包括定性和定量两种方法，定性方法主要通过揭示混沌信号在时域或频域中表现出的特殊空间结构或频率特性来判别，这种方法简单直观，但是过于笼统 定量方法通过计算混沌信号奇异吸引子的特性参数来辨别混沌行为的方法主要有两个： (1)描述邻近轨道发散率的Laypunov指数 (2)描述吸引子维数的关联维数和反映信息产生频率的Kolmogorov熵Lyapunov指数 混沌系统初值敏感性是指相空间中初始距离很近的两条轨迹会以指数速率发散，Lyapunov指数即是根据相轨迹有无扩散运动特征来判别系统的混沌特性在相空间中，轨迹间的距离分别表现为线度、面积和体积 对一维映射xn+1 = f[xn],假设初始位置有两个相邻的状态x0和y0=x0+ε, 则经过n次迭代后，有 所以Lyapunov指数vLyapunov指数定义：设相轨迹上两点之间的初始距离为|E0|=ε，用|En|表示经过n次迭代后该两点之间的距离，则称 为系统Lyapunov指数。

当 λ>0 时，系统具有混沌特征时间序列Lyapunov指数的计算 在实际时间序列混沌识别中，通常只估计最大Lyapunov指数，下面介绍一种算法: 小数据量法 设时间序列{x1,x2,…,xN}, 嵌入维数m，时间延迟τ，重构相空间后有： 其中k≤N-(m-1)τ ,选取一距离阈值δ，对初始点x(0)，选取与其距离最近的点不妨设为x(k0)=z0(0)，当经过t0次迭代后，z0(t0)与x(t0)的距离大于δ，令时间序列Lyapunov指数的计算 取与x(t0)的距离小于δ的点x(k1)=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熵越大，那么信息的损失速度越大，系统的混沌程度越大。

混沌相空间重构理论v相空间重构是对一维时间序列，按照延迟时间和嵌入维数重构一个与原动力系统等价的相空间;v相空间重构的理论基础是Takens相空间嵌入定理，即通过对一个m维流形上的连续流的观测值h(t)的延迟(h(t),h(t+τ ),…,h(t+2mτ))，可以将该连续流光滑嵌入到2m+1维空间(相空间)中；v相空间重构过程中有两个参数选取特别重要: 延迟时间τ 和嵌入维数.　延迟时间间隔　延迟时间间隔τ 的计算的计算主要方法v 线性自相关函数法v平均互信息法延迟时间间隔延迟时间间隔τ 的选取的选取线性自相关函数法定义自相关函数为 选择使得自相关函数C(τ)第一次为零时的 τ 的值缺点：对嵌入维数大于2时不是最优的延迟时间间隔延迟时间间隔τ 的选取的选取平均互信息法选择使I(τ ) 为第一个局部极小的τ为延迟时间间隔缺点：要对概率密度进行估计，需要的样本点多，计算复杂，且对嵌入维数大于2维不是最优的嵌入维数嵌入维数m的计算的计算主要方法vLyapunov指数法v关联积分法嵌入维数嵌入维数m的计算的计算Lyapunov指数法（1）首先确定延迟τ ；（2）依次递增取嵌入维数m，计算嵌入相空间(x(n),x(n-τ ),…,x(n-(m-1)τ))的Lyapunov指数，选取指数趋于常数的m作为嵌入相空间维数。

缺点：当样本点数不多时可信度较低，当嵌入维数较高时面临“维数灾”嵌入维数嵌入维数m的的计算v关联积分: 定义：设t 为时间序列的延时，m是嵌入维数，N是样本数据集的大小，数据点个数M=N-(m-1)t ，则重构相空间中嵌入时间序列Y(i)的关联积分定义为： 关联积分是一个累积分布函数，表示相空间中任意两点之间距离小于半径r的概率，这里点与点之间的距离用矢量之差的无穷范数表示嵌入维数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的增加而在一定误差范围内不变为止缺点：当样本点数不多时可信度较低，当嵌入维数较高时面临“维数灾”总 结 关于混沌判定，一般应用最大Lyapunov指数，或者Kolmogorov熵，或者结合两者判定。

在计算延迟时间方法上，常用的有自相关函数法和互信息法，计算嵌入维常用的方法有G-P算法，这种方法是先计算出混沌时间序列的关联维，然后再计算出嵌入维数自适应滤波噪声消除v自适应滤波是一种采用参数可调整的有限脉冲响应滤波器进行滤波；v相对于固定参数滤波器，自适应滤波器有可能达到更好的去除噪声的效果；v相对于主成分和独立分量分析，自适应滤波算法简单；v自适应滤波需要有欲滤除噪声的参考信号自适应滤波算法一v可以考虑由每一次的采样值修正参数w，由v为减小该式误差，可以按照负梯度方向修正w，即v考虑到稳定性，一般μ应比较小，满足自适应滤波算法二v假设有n个样本，考虑遗忘因子加权的误差v令其关于w的梯度为零得：v得到自适应滤波算法二v令v有自适应滤波总结v自适应滤波需要确定参考信号；v需要确定参考信号的FIR滤波器阶数M；v计算复杂度低；v有针对单样本的算法；v在EEG信号处理中的应用，如：™对滤除眨眼的眼电干扰信号，参考信号可以考虑从FP1和FP2上取；™ERP中，参考信号可以考虑各段的平均值；™在去除心电干扰中，可以取心电信号作为参考信号主成分分析(K-L变换)v信号的DFT和DCT是将信号投影或表示到固定正交函数系上的变换v主成分分析或称K-L变换是将信号表示到一组标准正交信号系上的变换，该信号系一般由信号的统计特性得到v该正交信号系张成的空间称为信号子空间，其正交补空间称为噪声子空间v主成分分析或K-L变换的应用™用于分离信号与噪声™对信号的维数进行压缩降维主成分的定义及导出v设 x(n)=(x1(n),x2(n),…,xp(n))T 为一个 p 维平稳随机信号，E(x(n))=0，其协方差矩阵为 该矩阵为实对称矩阵，其非负特征值设为 ，则存在正交矩阵U，使得令：则有： ， ，因此 y(n) 的任意两个分量不相关。

y(n)的分量称为X的主分量由于总方差(平均能量)中属于第 i 主成分 yi 的比例为 称为主成分 yi 的贡献率主成分的计算v对原始信号数据v由其奇异值分解，知存在 p 阶正交矩阵U和 N 阶正交矩阵V，使得 其中v样本协方差矩阵v因此，矩阵R的特征值为v主分量为vV的前p行为正交信号，由其张成的空间为信号子空间主成分分析在降维中的应用v首先对样本协方差矩阵R计算正交变换矩阵U和特征值v计算前 m 个主成分的贡献率之和v对给定的阈值α(一般在0.8~0.95), 选取累计贡献达到α最小的mv由 y(n)=Umx(n) (Um为矩阵U的前m行子矩阵)计算得到m维向量，可用来代替p维向量x(n)，从而达到降维的目的，而信息的损失却不多主成分分析在降噪中的应用v首先对样本数据矩阵X计算正交矩阵U和Vv对实际采样信号xi(n)计算分解v也可以采用最小二乘算法(LMS)计算分解，即由主成分分析总结v要求信号服从高斯分布v各个主分量之间不相关，但不保证独立独立分量分析 v独立分量分析的目的是：当 X = A S 时，求矩阵 W，使得 Y=WX 的各个分量独立，此时W可能不是A的逆，但是 WA 是置换矩阵。

v典型例子： “鸡尾酒会”问题，从酒会嘈杂人声中提取所关心对象的语音, 人的大脑可以很快辨出或集中听某种需要关注声音麦克风1麦克风2麦克风3独立分量分析主要研究结构:(1) 美国加州大学生物系，计算神经生物实验室，提出信息极大化( infomax )算法 .(2) 日本Riken的数量神经科学实验室，互信息极小化（（minimization of mutual information MMI)采用人工神经网络优化 (3) 芬兰赫尔辛基工业大学神经网络研究中心，提出了立足于逐次提取独立分量的固定点算 法（（fixed point algorithm) ：fastICA www.cis.hut.fi/~oja (4) 法国学者：，提出了JADE算法、批数据处理算法、近年来引人注意的稀疏分量分析 http://tsi.enst.fr/~cardoso. ICA相关的基本概念vn阶矩(moment)：mn=E(xn) v特征函数v第二特征函数vn阶累计量(cumulant) v对高斯型信号，二阶以上的累计量都为0，因此可由一、二阶统计特征来完整描述。

vk4>0 超高斯, k4<0 亚高斯, 用︱k4︱大小作为衡量信号距离高斯型程度的度量均值 方差 峭度ICA相关的基本概念v联合矩中最常用的是协方差；v联合累计量最常用的是4阶累积量，一般用cum(x1,x2,x3,x4)表示v联合累积量性质：™当x各分量相互独立时，其互累计量恒等于0™比例性：cum(w1x1,w2x2,w3x3,w4x4)=w1w2w3w4cum(x1,x2,x3,x4)™当两个随机变量x, y独立时，有k4(x+y)=k4(x)+k4(y)v随机向量x=(x1,x2,…,xn)各个分量独立当且仅当其联合概率分布(密度)等于各个分量概率分布(密度)的乘积ICA计算数据的预处理解混系统解混系统B白化白化W正交系统正交系统UZ(t)混合系统混合系统Ax(t)y(t)s(t)系统简图系统简图 一般在一般在ICA求解之前，先对数据进行白化处理，这样可以求解之前，先对数据进行白化处理，这样可以使得各个分量之间是不相关的，便于进行解混计算，如下图使得各个分量之间是不相关的，便于进行解混计算，如下图所示所示样本数据的白化v对原始信号数据v假设数据已经零均值化，先计算数据的协方差矩阵v计算其特征值和特征向量，得到对角化分解：Rx=UTΛUv计算白化矩阵W： ，计算白化数据：Z=WX， 此时z的协方差矩阵是单位阵，即：Rz=I。

fastICA-基于非高斯极大化的算法v假设数据经过白化处理，即m维随机向量z的均值向量为零向量，协方差矩阵为单位阵；v采用一次提取一个分量的方法，以极大化非高斯性作为优化的目标函数，可以考虑以下两种度量™极大化峭度(Kurtosis)的绝度值™极小化负熵(Negentropy)v由于峭度：Kurt(y)=E[y4]-E[y2]2 的鲁棒性不好，一般 fastICA 算法采用负熵作为目标函数vfastICA算法也称为固定点算法或投影寻踪算法牛顿迭代法v牛顿法迭代法是用于求解方程 f(x)=0 的解v设在某一点Pk，对应x坐标为xk，其切线方程为v与x轴的交点为xk+1，则有迭代公式例.用牛顿迭代法求方程的根:解:由牛顿迭代法x0 =0.5;x1 =0.3333333333x2 =0.3472222222x3 =0.3472963532x4 =0.3472963553迭代四次精度达10-8 负熵负熵v信息熵： 当随机变量的均值和方差一定时，其 pdf 为高斯分布时其信息熵最大v负熵：任意 pdf 的 p(x) 和具有相同协方差矩阵的高斯分布pG(x)的差作为该 pdf 非高斯程度的度量，定义为负熵 J(x)≥0，p(x)=pG(x)：当且仅当J[p(x)]=0代表高斯分布。

v 由以上公式可知，计算负熵需要随机变量的概率密度，在实际计算中是不方便的，因此需要给出其计算的近似公式负熵的近似负熵的近似v高阶统计量形式：设x零均值，方差为1(白化数据)™Edgeworth展开™Gram-Charlier展开缺点：缺点：大值野点会引大值野点会引起较大误差起较大误差负熵的近似负熵的近似v非多项式函数形式： 其中G(1)(y)为奇函数，G(2)(y)为偶函数，v为与y有相同均值与方差的高斯随机变量可以考虑如下更简洁的形式vfastICA算法就以上式为基础设y=wTz，其关于w的梯度为v所以其极值点使 gradw(E[G(y)])=E[G’(wTz)z]=0 基于负熵的牛顿迭代公式基于负熵的牛顿迭代公式v考虑到约束wTw=1，由拉格朗日乘子法得极值点满足： v考虑映射： v由线性近似：v得到迭代公式：基于负熵的牛顿迭代公式基于负熵的牛顿迭代公式v于是v考虑到单位化，可以简化为v由于G(y)为偶函数，其导数g(y)为奇函数，一般可选 基于负熵的基于负熵的fastICA基于负熵提取一个分量的固定点算法（基于负熵提取一个分量的固定点算法（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）若）若p

互信息与随机变量的独立性互信息与随机变量的独立性vKL散度(kullback-leibler): 描述的是两个概率密度函数间相似程度 KL≥0；当 p(x)=q(x) 时, KL=0v互信息：当p(x)为多变量x=[x1,x2,…xm]T的联合 pdf，p(xi)为各分量的边际 pdf，则互信息(mutual information)定义为v互信息性质：™I(x) ≥0，I(x)=0当且仅当x 中各分量相互独立；™ 互信息与随机变量的独立性互信息与随机变量的独立性v假设数据经过白化处理，即m维随机向量x的均值向量为零向量，协方差矩阵为单位阵v设y=Bx，假设x的 pdf 为 px(x)，y的 pdf 为 py(y)则有 py(y)=px(B-1y)|det-1(B)|v于是有互信息与随机变量的独立性互信息与随机变量的独立性考虑考虑y的互信息的互信息输入信号的熵与B无关等效于等效于因此因此(y1,y2,…,ym)的独立性可以通过极大化下式得到的独立性可以通过极大化下式得到基于互信息的迭代公式基于互信息的迭代公式v考虑梯度 其中 即： 自然梯度： 迭代公式：基于互信息的迭代公式基于互信息的迭代公式v很明显，函数 的选取与源信号的pdf有关，一般只要区分超高斯或亚高斯即可，即 超高斯：g+(y)=-2tanh(y) 亚高斯：g-(y)=tanh(y)-yv在算法中要动态地对源信号的概率密度进行估计，可以采用如下判据： E[-tanh(yi)yi+(1-tanh(yi)2)] 当其大于零时为超高斯，小于零时为亚高斯。

Informax算法算法（1）对数据进行中心化使其均值为零；（2）随机生成分离矩阵B和指标γi，设置学习率μ、（3）计算y=Bz（4）如果非线性形式没有事先确定（a）γi=(1-μγ)γi+μγ(E[-tanh(yi)yi+(1-tanh(yi)2)])（b）如果γi>0，取gi(yi)=-2tanh(yi)，否则g(yi)=tanh(yi)-yi（5）更新分离矩阵B←B+γ[I+g(y)yT]B（6）如果尚未收敛，返回步骤（3）独立分量分析 v主成分分析的局限性：在主成分分解 Y=UX 中，当X不服从正态分布时，Y的各个分量是不相关的，但不能保证是独立的当 X 是独立信号的混合时，即 X = A S，主成分分析得不到 S v独立分量分析的目的是：当 X = A S 时，求矩阵 W，使得 Y=WX 的各个分量独立，此时W可能不是A的逆，但是 WA 是置换矩阵v由于生物信号一般具有非平稳、非正态等性质，因此ICA比PCA更有优势v独立分量分解的局限性：求解ICA的计算复杂度比PCA高，理论深奥算法复杂，各个分量需要解释判读胎儿心电提取胎儿心电提取8通道原始波形通道原始波形>> load foetal_ecg.dat –ascii>> X=foetal_ecg(:,2:9);>> plot(X(1:1000,1);…ICA分解分解>>Y=fastica (X’)’;>>plot(Y(1:1000,1));…8个分量波形个分量波形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;Common Spatial Pattern (CSP)v设 和 为代表两个类的两个p 维随机向量，E(X)=0， E(Y)=0，其协方差矩阵分别为 若这两个矩阵都是正定矩阵，则存在矩阵 Q 满足：实际上，Q的列向量和 为广义特征值问题 的特征向量和特征值。

实际上，存在矩阵 G 满足： 于是 为实对称正定矩阵，因此存在正交矩阵 P 和非负特征值 使得：令 Q = G-TP , 则有另一方面因此有定义变换：z =QT x，则向量 z 的分量方差可以更好地用于求解两分类问题 如果采用C语言编程，可以按照以上过程，结合求解LDLT分解和特征值特征向量问题求得Q求解 Q 的matlab程序如下：。