第二讲 EEG信号处理基础

1、第二讲 EEG信号处理基础主讲人：谢宏主讲人：谢宏信息工程学院信息工程学院vEEG信号模型vEEG信号特征EEG信号的特点v随机信号v非平稳v非线性v非高斯过程EEG信号模型v基于神经元的生化物理模型Hodgkin and Huxley 模型MorrisLecar 模型v基于随机信号的动态模型线性模型：AR模型和ARMA模型非线性模型：GARCH模型EEG信号线性模型vAR模型)()()(1nxknyanypkk), 0()(2NnxvARMA模型qkpkkknxknyany01)()()(), 0()(2NnxvAIC准则确定模型阶数：即最小化如下目标函数)(2)ln(),(2,qpNqpA

2、ICqpv多维AR模型：对应多通道)()()()(111nxknyaknyanyimijjpkjjkpkiiki诱发电位的Prony方法方法v首先建立输入为脉冲函数的AR模型：)()() 1()(1nApnanyanypv假设有样本值y(1),y(2),y(N), 令)()2() 1()2()(1() 1 () 1()(pNYNyNyypypyypypyY)()2() 1(NYpYpyyv解得AR模型系数：yYaaaap21v求解特征方程：0111ppppaaa 得到其p个根pjpjjeAeAeA,2121v方程的通解为：njnppnjnnjnpeAweAweAwny212211)(v求解其中

3、kjkkeAr1) 1() 1 ()0(111211121121pyyywwwrrrrrrRppppppEEG信号非线性模型vAR模型)()()(1nxknyanypkk), 0()(2NnxvARMA模型qkpkkknxknyany01)()()(), 0()(2NnxvAIC准则确定模型阶数：即最小化如下目标函数)(2)ln(),(2,qpNqpAICqpv多维AR模型：对应多通道)()()()(111nxknyaknyanyimijjpkjjkpkiiki由模型化神经元活动产生EEG信号vAR模型vARMA模型由模型化神经元活动产生EEG信号vAR模型EEG信号常用统计特征v随机变量或向

4、量统计量均值、方差、协方差与相关系数矩阵偏度(skewness)：峭度或峰度(kurtosis)：v随机过程或随机信号统计量自相关函数与功率谱双谱等高阶谱EEG信号常用统计特征v设m个电极通道采集的脑电信号表示为Tmnxnxnxnx)()()()(21) 1() 1 ()0() 1() 1 ()0() 1() 1 ()0(222111NxxxNxxxNxxxXmmm11)(1)(NkiiikxNnxEv均值v其N点的采样数据为EEG信号常用统计特征v方差11222)(1)(NkiiiiikxNnxEv偏度113333/)(1)()(NkiiiiiiikxNnxEnxSkewv峭度(峰度)3)(

5、)(24iiiinxEnxKurtv协方差矩阵)()()()(TnxEnxnxEnxEA1)()()(1)()()(NkjjiijjiiijkxkxNnxnxEaEEG信号常用统计特征v相关系数矩阵v注意：以上统计量只涉及随机变量取值，不涉及时序计算时要考虑实时性，如均值的计算，可以考虑采用算法)(ijrR 1)()(1)()(NkjjjiiijjjiiiijkxkxNnxnxEr) 1(1) )(11(1)(1)(2010NxNkxNNNkxNnxEiNkiNkiii)(1) 1(1)(nxnnnnniii平稳EEG信号时序特征v自相关函数v功率谱：信号的功率谱是其自相关函数的傅里叶变换，即

6、10)()(1)()()(kNniiiiiknxnxNknxnxEkv多通道信号的自相关函数矩阵10)()(1)()()(kNnTTknxnxNknxnxEkRkkjiikkjiieknxnxEekS)()()()()()()()(injikkjiiXenxEeknxnxE2*| )(|)()()()(iiinnjiiXXXenxX平稳EEG信号时序特征v功率谱的计算：周期图法与线性时序模型法v周期图法计算频谱得到，简单、功率谱不光滑、误差大；v假设EEG信号的时序模型为AR(p)，即)()()(1nknxanxpkkv其功率谱为212|1 |)(pkkjkpeaS平稳EEG信号时序特征vEE

7、G信号节律波能量提取： 设信号x(n)的傅里叶变换为X()，则有对有限个采样值x(0),x(1),x(N-1), 其DFT为X(k), 对应频率kfs/N, 则有：dXnxin2022| )(|21| )(|sff2则EEG信号节律波平均能量提取算法可以考虑两种方案。102102| )(|1| )(|NkNnkXNnx平稳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第一步：对

8、脑电信号x(n)进行带通fL,fH滤波得到输出信号y(n);v第二步：HLNNkavgkXNE2| )(|1102| )(|1NnavgnyNEEEG信号滤波v信号滤波涉及：低通、高通、带通、陷波v滤波器的比较：IIR满足相同特性阶数较低，只能近似线性相位，必须浮点运算FIR满足相同特性阶数较高，可以做到严格线性相位，可以采用整数运算v滤波器一般要结合实现时的计算效率和滤波器特性等综合考虑阶数通带、阻带和过渡带特性延迟IIR滤波器vIIR滤波器的模型为：vIIR滤波器的类型：贝塞尔、巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆型v相位特性：贝塞尔巴特沃斯切比雪夫椭圆v过渡带宽度：贝塞尔巴特沃斯

9、切比雪夫0发散，a=0为常数，a3.6时进入混沌状态时进入混沌状态64正 文v分析论证过程65罗沦兹方程 其中其中x x、y y、z z为无量纲量为无量纲量，分别表征对流强度，对流中，分别表征对流强度，对流中升流与降流间的温差和竖直方升流与降流间的温差和竖直方向温度分布的非线性度。任意向温度分布的非线性度。任意给定初值，系统最终都会回到给定初值，系统最终都会回到状态空间的特定区域内，其吸状态空间的特定区域内，其吸引子具有精巧而奇特的结构，引子具有精巧而奇特的结构，如图所示，表明系统进入了混如图所示，表明系统进入了混沌状态。沌状态。( (图为图为LorenzLorenz混沌吸混沌吸引子引子) )

10、混沌识别 混沌识别主要包括定性和定量两种方法，定性方法主要通过揭示混沌信号在时域或频域中表现出的特殊空间结构或频率特性来判别，这种方法简单直观，但是过于笼统。 定量方法通过计算混沌信号奇异吸引子的特性参数来辨别混沌行为的方法。主要有两个： (1)描述邻近轨道发散率的Laypunov指数 (2)描述吸引子维数的关联维数和反映信息产生频率的Kolmogorov熵Lyapunov指数 混沌系统初值敏感性是指相空间中初始距离很近的两条轨迹会以指数速率发散，Lyapunov指数即是根据相轨迹有无扩散运动特征来判别系统的混沌特性。在相空间中，轨迹间的距离分别表现为线度、面积和体积。 对一维映射xn+1 =

11、 fxn,假设初始位置有两个相邻的状态x0和y0=x0+, 则经过n次迭代后，有 所以)()(11nnnnnxfyfxyE11111)()(nnnnnExfxyxfnnkknnnnnexfEEEEEEEE10012110)(Lyapunov指数vLyapunov指数定义：设相轨迹上两点之间的初始距离为|E0|=，用|En|表示经过n次迭代后该两点之间的距离，则称 为系统Lyapunov指数。 当 0 时，系统具有混沌特征。10110)(ln1limln1limnkknnkkknxfnEEn时间序列Lyapunov指数的计算 在实际时间序列混沌识别中，通常只估计最大Lyapunov指数，下面介绍

12、一种算法: 小数据量法。 设时间序列x1,x2,xN, 嵌入维数m，时间延迟，重构相空间后有： 其中kN-(m-1) ,选取一距离阈值，对初始点x(0)，选取与其距离最近的点不妨设为x(k0)=z0(0)，当经过t0次迭代后，z0(t0)与x(t0)的距离大于，令mTmkkkRxxxkx),()()1(|)0()0(|)()(|ln1000000 xztxtztl时间序列Lyapunov指数的计算 取与x(t0)的距离小于的点x(k1)=z1(0)，且z1(0)-x(t0)与z0(t0)-x(t0)的夹角最小，如图所示重复以上过程，假设直到M-1步结束，此时tM-1=N-(m-1)，则有 10

13、1MkklMKolmogorov熵 Kolmogorov熵K在混沌的度量中有着相当重要的应用。对于规则运动，K=0; 随机系统K=无穷，若系统表现确定性混沌 ， 则 Kolmogorov熵是大于0的常数。Kolmogorov熵越大，那么信息的损失速度越大，系统的混沌程度越大。混沌相空间重构理论v相空间重构是对一维时间序列，按照延迟时间和嵌入维数重构一个与原动力系统等价的相空间;v相空间重构的理论基础是Takens相空间嵌入定理，即通过对一个m维流形上的连续流的观测值h(t)的延迟(h(t),h(t+ ),h(t+2m)，可以将该连续流光滑嵌入到2m+1维空间(相空间)中；v相空间重构过程中有两

14、个参数选取特别重要: 延迟时间 和嵌入维数.延迟时间间隔延迟时间间隔 的计算的计算主要方法v 线性自相关函数法v平均互信息法延迟时间间隔延迟时间间隔 的选取的选取线性自相关函数法NnnNnnNnnnxNxxxNxxxxNC11211,)(1)(1)(其中定义自相关函数为 选择使得自相关函数C()第一次为零时的 的值。缺点：对嵌入维数大于2时不是最优的。延迟时间间隔延迟时间间隔 的选取的选取平均互信息法为概率，其中定义平均互信息为对于时间序列),()(,)()(),(log),()(,21nnnnnnnnnNnnxxPxPxPxPxxPxxPIx选择使I( ) 为第一个局部极小的为延迟时间间隔。

15、缺点：要对概率密度进行估计，需要的样本点多，计算复杂，且对嵌入维数大于2维不是最优的。嵌入维数嵌入维数m的计算的计算主要方法vLyapunov指数法v关联积分法嵌入维数嵌入维数m的计算的计算Lyapunov指数法（1）首先确定延迟 ；（2）依次递增取嵌入维数m，计算嵌入相空间(x(n),x(n- ),x(n-(m-1)的Lyapunov指数，选取指数趋于常数的m作为嵌入相空间维数。缺点：当样本点数不多时可信度较低，当嵌入维数较高时面临“维数灾”。嵌入维数嵌入维数m的的计算v关联积分: 定义：设 为时间序列的延时，m是嵌入维数，N是样本数据集的大小，数据点个数M=N-(m-1) ，则重构相空间中

16、嵌入时间序列Y(i)的关联积分定义为： 关联积分是一个累积分布函数，表示相空间中任意两点之间距离小于半径r的概率，这里点与点之间的距离用矢量之差的无穷范数表示。0,01,00,| , ( )ZijijZrdYYZMjijidrMMrC1)() 1(2)(嵌入维数m的计算关联积分法的主要步骤(1)利用时间序列X1,X2,XN,先给定一个较小的嵌入维数m0,重构相空间，得到新的序列Yi(2)计算关联积分C(r)(3)对于r的某个取值范围，吸引子的维数d与累积分布函数C(r)应满足对数线性关系，即d(m)=LnC(r)/Lnr,从而可用最小二乘拟合得到对应于m0的关联维数估计d(m0)(4)增加嵌入

17、维数m0,重新计算步骤(2)和(3),直到相应的维数估计值d(m)不再随着m的增加而在一定误差范围内不变为止。缺点：当样本点数不多时可信度较低，当嵌入维数较高时面临“维数灾”。总 结 关于混沌判定，一般应用最大Lyapunov指数，或者Kolmogorov熵，或者结合两者判定。在计算延迟时间方法上，常用的有自相关函数法和互信息法，计算嵌入维常用的方法有G-P算法，这种方法是先计算出混沌时间序列的关联维，然后再计算出嵌入维数。自适应滤波噪声消除v自适应滤波是一种采用参数可在线调整的有限脉冲响应滤波器进行滤波；v相对于固定参数滤波器，自适应滤波器有可能达到更好的去除噪声的效果；v相对于主成分和独立

18、分量分析，自适应滤波算法简单；v自适应滤波需要有欲滤除噪声的参考信号。)()()( )()(nnxnvnxneTrw)()()()(2)(2wrrwwrnnnnxnxETTT)()()(22nnxEneETrwwrrwwr)()()()(2)(2nnEnnxEnxETTT) 1() 1()()(Mnrnrnrnr)()()()(1nnxEnnETrrrw自适应滤波算法一v可以考虑由每一次的采样值修正参数w，由22)()()(nnxneTrwv为减小该式误差，可以按照负梯度方向修正w，即)()(2) 1()() 1()(2nnendndennrwwww)()() 1()(2) 1(nnnnxnw

19、rrwv考虑到稳定性，一般应比较小，满足)()(max10nnETrr自适应滤波算法二v假设有n个样本，考虑遗忘因子加权的误差niinie02)()(wv令其关于w的梯度为零得：v得到niiniiedd0)()(2)(rwwniTiniiix0)()()(2rrw0wrrrniniTininiiiix11)()( 2)()( 2niinniTiniixii111)()( )()(rrrw自适应滤波算法二v令v有)()() 1()()()(1nnniinTniTinrrRrrR)()()() 1()() 1()(1nnnnnxnnTrRrwww11)()() 1()(nnnnTrrRR)() 1

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

21、统计特性得到v该正交信号系张成的空间称为信号子空间，其正交补空间称为噪声子空间v主成分分析或K-L变换的应用用于分离信号与噪声对信号的维数进行压缩降维主成分的定义及导出v设 x(n)=(x1(n),x2(n),xp(n)T 为一个 p 维平稳随机信号，E(x(n)=0，其协方差矩阵为 该矩阵为实对称矩阵，其非负特征值设为 ，则存在正交矩阵U，使得)()()()()()()()()()()()()()()()()(2212221212121nxEnxnxEnxnxEnxnxEnxEnxnxEnxnxEnxnxEnxEnxnxERpppppT)()(21TpTnUxnUxEURU021p令：)()

22、()()()()()(2121222211121121nxnxnxuuuuuuuuunUxnyynynypppppppp则有： ， ，因此 y(n) 的任意两个分量不相关。y(n)的分量称为X的主分量。由于iinyVar)(0)(),(cov(nynyjipiiTpiipiixVartracePPtraceyVar111)()()()(总方差(平均能量)中属于第 i 主成分 yi 的比例为 称为主成分 yi 的贡献率。1piii主成分的计算v对原始信号数据)()2() 1 ()()2() 1 ()()2() 1 (222111NxxxNxxxNxxxXpppv由其奇异值分解，知存在 p 阶正交

23、矩阵U和 N 阶正交矩阵V，使得 其中VUXH0S p1Sv样本协方差矩阵v因此，矩阵R的特征值为UUNUVVUNXXNRTTTHT2111SNii2v主分量为pSVUXYvV的前p行为正交信号，由其张成的空间为信号子空间)()2() 1 ()()2() 1 ()()2() 1 (222111NvvvNvvvNvvvVpppp主成分分析在降维中的应用v首先对样本协方差矩阵R计算正交变换矩阵U和特征值v计算前 m 个主成分的贡献率之和v对给定的阈值(一般在0.80.95), 选取累计贡献达到最小的m。v由 y(n)=Umx(n) (Um为矩阵U的前m行子矩阵)计算得到m维向量，可用来代替p维向量

24、x(n)，从而达到降维的目的，而信息的损失却不多。11pmiiii主成分分析在降噪中的应用v首先对样本数据矩阵X计算正交矩阵U和Vv对实际采样信号xi(n)计算分解)()()(1nnvwnxipkkikiNnkikiiknvnxnvnxnw1)()()(),()(v也可以采用在线最小二乘算法(LMS)计算分解，即由212)()()()()(iTipkkikiiwnvnxnvwnxniiidwndnwnw)() 1()(2)() 1()()(2) 1(nvnwnvnxnwTiTii主成分分析总结v要求信号服从高斯分布v各个主分量之间不相关，但不保证独立独立分量分析 v独立分量分析的目的是：当 X

25、 = A S 时，求矩阵 W，使得 Y=WX 的各个分量独立，此时W可能不是A的逆，但是 WA 是置换矩阵。v典型例子： “鸡尾酒会”问题，从酒会嘈杂人声中提取所关心对象的语音, 人的大脑可以很快辨出或集中听某种需要关注声音。)()()()()()()()()()()()(333232131332322212123132121111tsatsatsatxtsatsatsatxtsatsatsatx麦克风1麦克风2麦克风3)(1tx) (2tx)(3tx11a12a13a21a22a)(1ts)(2ts)(3ts23a31a32a33a独立分量分析主要研究结构:(1) 美国加州大学生物系，计算神

26、经生物实验室，提出信息极大化( infomax )算法。 .(2) 日本Riken的数量神经科学实验室，互信息极小化（minimization of mutual information MMI)采用人工神经网络优化。 http:/www.brain.riken.jp/lab/mns/amari(3) 芬兰赫尔辛基工业大学神经网络研究中心，提出了立足于逐次提取独立分量的固定点算 法（fixed point algorithm) ：fastICA。 www.cis.hut.fi/oja www.cs.helsinki.fi/aapo.hyvarinen(4) 法国

27、学者：J.F.Cardoso，提出了JADE算法、批数据处理算法、近年来引人注意的稀疏分量分析。 http:/tsi.enst.fr/cardoso. ICA相关的基本概念vn阶矩(moment)：mn=E(xn) v特征函数v第二特征函数vn阶累计量(cumulant) 0)()()(nnnnddjdxxpx0)()(nnnnddjk)(log)(11mk 2122mmkdxexpxj)()(31123323mmmmk4121213224461243mmmmmmmkv对高斯型信号，二阶以上的累计量都为0，因此可由一、二阶统计特征来完整描述。vk40 超高斯, k40 亚高斯, 用k4大小作为

28、衡量信号距离高斯型程度的度量。均值 方差 峭度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)混合系

29、统混合系统Ax(t)y(t)s(t)系统简图系统简图 一般在一般在ICA求解之前，先对数据进行白化处理，这样可以求解之前，先对数据进行白化处理，这样可以使得各个分量之间是不相关的，便于进行解混计算，如下图使得各个分量之间是不相关的，便于进行解混计算，如下图所示所示样本数据的白化v对原始信号数据)()2() 1 ()()2() 1 ()()2() 1 (222111NxxxNxxxNxxxXmmmv假设数据已经零均值化，先计算数据的协方差矩阵TxXXNR1v计算其特征值和特征向量，得到对角化分解：Rx=UTUv计算白化矩阵W： ，计算白化数据：Z=WX， 此时z的协方差矩阵是单位阵，即：Rz=I

30、。UW21fastICA-基于非高斯极大化的算法v假设数据经过白化处理，即m维随机向量z的均值向量为零向量，协方差矩阵为单位阵；v采用一次提取一个分量的方法，以极大化非高斯性作为优化的目标函数，可以考虑以下两种度量极大化峭度(Kurtosis)的绝度值极小化负熵(Negentropy)v由于峭度：Kurt(y)=Ey4-Ey22 的鲁棒性不好，一般 fastICA 算法采用负熵作为目标函数。vfastICA算法也称为固定点算法或投影寻踪算法。牛顿迭代法v牛顿法迭代法是用于求解方程 f(x)=0 的解。v设在某一点Pk，对应x坐标为xk，其切线方程为v与x轴的交点为xk+1，则有迭代公式)()(

31、kkkxfxxxfy)()(1kkkkxfxfxx例.用牛顿迭代法求方程的根:0133xx解:13)(3xxxf设33)(2xxf由牛顿迭代法)()(1kkkkxfxfxx得取初值,5 .00 xx0 =0.5;x1 =0.3333333333x2 =0.3472222222x3 =0.3472963532x4 =0.3472963553331323kkkkxxxx迭代四次精度达10-8 1kx*x)(xfy kx负熵负熵v信息熵： 当随机变量的均值和方差一定时，其 pdf 为高斯分布时其信息熵最大。v负熵：任意 pdf 的 p(x) 和具有相同协方差矩阵的高斯分布pG(x)的差作为该 pdf

32、 非高斯程度的度量，定义为负熵dxxpxpXH)(log)()()()()(xHxHxJG J(x)0，p(x)=pG(x)：当且仅当Jp(x)=0代表高斯分布。v 由以上公式可知，计算负熵需要随机变量的概率密度，在实际计算中是不方便的，因此需要给出其计算的近似公式。负熵的近似负熵的近似v高阶统计量形式：设x零均值，方差为1(白化数据)Edgeworth展开Gram-Charlier展开3434( )1( )( )( )3!4!Gp xkkHxHxpx ( ),nH x Hermite多项式2232344341 ( )47648J p xkkkk k2242343341 ( )431848J

33、p xkkkk k缺点：缺点：大值野点会引大值野点会引起较大误差起较大误差负熵的近似负熵的近似v非多项式函数形式： 其中G(1)(y)为奇函数，G(2)(y)为偶函数，v为与y有相同均值与方差的高斯随机变量。可以考虑如下更简洁的形式vfastICA算法就以上式为基础。设y=wTz，其关于w的梯度为v所以其极值点使 gradw(EG(y)=EG(wTz)z=0 2)2()2(22)1(1)()()()(vGEyGEkyGEkyJ2)()()(vGEyGExJ)(grad)()(2)(gradyGEvGEyGEyJww基于负熵的牛顿迭代公式基于负熵的牛顿迭代公式v考虑到约束wTw=1，由拉格朗日乘

34、子法得极值点满足： v考虑映射： v由线性近似：v得到迭代公式：wzzwgEwFT)()()()(yGyg)()()(kkwTwFwwdwwdFFk)()(11kwTkkwFdwwdFwwkIzzzwgEdwwdFTTT)()(IzwgT)(0 wzzwgET)(IzzEzwgTT)(基于负熵的牛顿迭代公式基于负熵的牛顿迭代公式v于是v考虑到单位化，可以简化为v由于G(y)为偶函数，其导数g(y)为奇函数，一般可选 )()(11kTkTkkkwzzwgEzwgEww)()()(1zzwgEwzwgEzwgETkkTkTkkTkTkkwzwgEzzwgEw)()(1332211)()2/exp(

35、)()tanh()(yygyyygyayg2322212113)()2/exp()1 ()()(tanh1 ()(yygyyygyaayg基于负熵的基于负熵的fastICA基于负熵提取一个分量的固定点算法（基于负熵提取一个分量的固定点算法（2001）（1）X去均值，进行白化得去均值，进行白化得Z；（2）任选）任选 w1 的初值，使的初值，使 w1的的2-范数为范数为1，k=1；（3）令）令 ，均值用样本，均值用样本均值；均值；（4）归一化：）归一化：wk+1=wk+1/|wk+1|;（5）若）若 不接近不接近1，返回步骤（，返回步骤（3），否则迭代结束），否则迭代结束（6）提取分量）提取分量该方法收敛具有三阶收敛速度该方法收敛具有三阶收敛速度kTkTkkwzwgEzzwgEw)()(1|1kTkwwZwyT基于负熵的基于负熵的fastICA基于负熵逐个提取多个分量的固定点算法基于负熵逐个提取多个分量的固定点算法（1）X去均值，进行白化得去均值，进行白化得Z；（2）m为待提取独立分量数目，为待提取独立分量数目，p=1；（3）任选）任选 wp,1 的初值，使的初值，使 wp,1 的的2-范数为范数为1，k=1；（4）令）令 ，均值用样本均值；，均值用样本均值；（5）正交化：）正交