CN104809434A - 一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法 - Google Patents

Abstract

一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，本发明涉及眼电伪迹去除的睡眠分期方法。传统方法缺乏参考眼电信号，去除方法较为困难、知识的抽取和表达比较困难，并且其收敛性有时无法保证，无法对知识直接进行学习以及不能充分挖掘脑电信号中蕴含的睡眠信息的问题，而提出的一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法。该方法是通过1、得到M个小波系数；2、不包含眼电伪迹的P个小波系数，作为纯净的脑电信号；3、得到IMF分量的个数与残差之和；4、得到脑电分量以及眼电分量；5、重构纯净的脑电信号；6、得到X(n)；7、提取7个特征参数；8、得到睡眠状态分期指数等步骤实现的。本发明应用于睡眠分期领域。

Description

一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法

技术领域

本发明涉及眼电伪迹去除的睡眠分期方法，特别涉及一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法。

背景技术

睡眠作为人体休息的自然状态，约占人类寿命三分之一的时间，因此睡眠分期对于研究睡眠对神经功能的影响、睡眠障碍、睡眠和其他医学状态紊乱的关系具有十分重要的意义。脑电信号(EEG)是大脑细胞轴突与树突之间以及细胞之间的电压变化在大脑头皮上的反映，包含了丰富的脑功能状态信息。因此，对EEG的分析具有非常重要的实际意义。在睡眠相关脑电信号分析过程中，眼电信号伪迹去除，特征参数提取与分类是其中最重要的几个技术环节。

经大脑头皮采集到的脑电信号包含各种电路噪声，工频干扰，以及生理伪迹干扰信号。而这些干扰信号往往幅度较大使得真实的脑电信号被湮没在各种噪声之中。其中，眼电伪迹是影响脑电信号的重要生理伪迹之一。多通道脑电信号的眼电伪迹去除方法相对较为成熟和完善，而单通道脑电信号，由于信息量少，而且缺乏参考眼电信号，去除方法较为困难。

目前，针对睡眠状态分期，常用的脑电信号特征提取技术多集中在脑电信号的频域分析。但是最新的研究表明，脑电信号具有严重的非线性、非平稳等特征，仅仅依靠频域的特征参数，势必不能充分提取脑电信号所包含的睡眠信息。常用的非线性特征参数分类方法主要神经网络和模糊逻辑等。神经网络分类法具有较强的非线性映射能力，虽然可以对训练样本进行有效的学习，但对知识的抽取和表达比较困难，并且其收敛性有时无法保证。模糊逻辑分类法：采用规则库的方法很好地解决了知识的表达和抽取这一问题，但无法对知识直接进行学习。

综上所述，针对睡眠状态分期，现有的脑电特征提取方法只注重其某一特征的分析，不能充分挖掘脑电信号中蕴含的睡眠信息；在此基础上的脑电特征分类方法不能完全解决知识表达、抽取和知识学习的问题，分类效果不佳。

发明内容

本发明的目的是为了解决传统方法缺乏参考眼电信号，去除方法较为困难、知识的抽取和表达比较困难，并且其收敛性有时无法保证，无法对知识直接进行学习以及不能充分挖掘脑电信号中蕴含的睡眠信息的问题，而提出的一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法。

上述的发明目的是通过以下技术方案实现的：

步骤一、对采集到的单通道脑电信号X(n)进行小波变换，得到M个小波系数；其中，M个小波系数分为两类：不包含眼电伪迹的小波系数为P个以及含眼电伪迹的小波系数为M-P个；

步骤二、对于不包含眼电伪迹的P个小波系数，直接作为纯净的脑电信号；

步骤三、含有眼电伪迹的M-P个小波系数W(i)经过经验模态分解后表示为N个内蕴模态函数即IMF分量的个数与残差之和：

       $W\left(i\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{C}_j+r---\left(1\right)$

其中：W(i)为第i个含有眼电伪迹的小波系数，C_j为该小波系数的第j个内蕴模态函数，r为残差，i＝1,2,3，…，M-P，j＝1,2,3，…，N；

步骤四、将步骤三得到的每个含有眼电伪迹的小波系数的所有内蕴模态函数C_j以及残差r输入到独立分量分析系统ICA中，得到每个含有眼电伪迹的小波系数的脑电分量以及眼电分量；

步骤五、将步骤四得到的每个小波系数的眼电分量置零，并采用ICA逆运算重构纯净的脑电信号；

步骤六、将步骤二和步骤五中的得到的纯净脑电信号求和，得到去除眼电伪迹后的纯净脑电信号X(n)；

步骤七、对步骤六得到的纯净脑电信号X(n)提取7个特征参数即4个频域参数、2个非线性参数以及1个时频参数；其中，4个频域参数包括δ，θ，α和β节律的频带能量比即ratio(δ)、ratio(θ)、ratio(α)和ratio(β)，2个非线性参数包括复杂度complexity简称c(n)和最大Lyapunov指数简称Lyapunov，时频参数为希尔伯特黄变换瞬时频率均值简称hilber；

步骤八、将步骤七提取的3个频域参数、2个非线性参数以及1个时频参数到的特征参数输入到自适应模糊神经推理系统中，得到睡眠状态分期指数；其中，3个频域参数为在步骤七提取的4个频域参数选择3个频域参数；即完成了一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法。

发明效果

本发明的目的是为了解决基于单通道脑电信号的睡眠状态分期过程中已有方法无法有效去除眼电伪迹的干扰，且不能充分利用脑电信号中蕴含的睡眠信息以及分类效果不佳的缺点，提出了一种结合小波变换，经验模态分解以及独立成分分析的单通道脑电信号眼电伪迹去除方法，同时还提出了一种综合脑电信号频域分析、非线性分析以及时频分析的特征参数提取方法，结合自适应模糊神经推理系统(ANFIS)分类器实现睡眠分期。从而解决不能充分挖掘脑电信号中蕴含的睡眠信息的问题。

在单通道脑电信号眼电伪迹去除的研究方面，本发明提出了基于小波变换，经验模态分解和独立分量分析的WT-EMD-ICA方法，有效的去除了脑电信号中的眼电伪迹，而且使得有用的脑电信息得以充分保持；在与睡眠状态分期相关的脑电信号特征参数提取方面，目前常用的脑电特征提取方法只注重其某一特征的分析，不能充分挖掘脑电信号中蕴含的睡眠信息，本发明综合脑电信号频域分析、非线性分析以及时频分析的特征参数提取方法，提取睡眠脑电信号中4个节律的‘频带能量比’参数，2个非线性参数以及1个时频参数，进一步提取了脑电信号中与睡眠分期相关的信息同时，本发明采用具有较好知识表达和学习能力的自适应模糊神经推理系统，得到睡眠状态分期指标具有计算速度快，准确度高等优点。如图7～9所示，相比于偏最小二乘分类器(PLS)和最小二乘支持向量机分类器(LS-SVM)，自适应模糊神经推理系统分类器(ANFIS)的睡眠分期正确率明显提高。

附图说明

图1为具体实施方式一提出的中基于单通道脑电信号的眼电伪迹去除框图；

图2为具体实施方式一提出的睡眠脑电特征参数提取与分类框图；

图3具体实施方式九提出的为特征参数分类的自适应模糊神经推理系统网络结构图；其中，complex，hilbert，lyapunov，ratio(δ)，ratio(α)，ratio(β)分别为复杂度，希尔伯特瞬时频率均值，最大lyapunov指数，δ频带能量比，α频带能量比，β频带能量比，A，B，C，D，E，F分别为其对应的隶属度函数，符号П表示对该节点的参数求乘法运算，符号N表示对该节点的参数进行归一化处理；

图4(a)为实施例提出的应用本发明所设计的ANFIS网络对MIT-BIH睡眠数据库中slp01样本进行分类学习前参数ratio(δ)的初始隶属度函数图；其中，横坐标为输入参数ratio(δ)的范围，纵坐标为0～1的无量纲隶属度；

图4(b)为实施例提出的应用本发明所设计的ANFIS网络对MIT-BIH睡眠数据库中slp01样本进行分类学习后参数ratio(δ)的最终隶属度函数图；其中，横坐标为输入参数ratio(δ)的范围，纵坐标为0～1的无量纲隶属度；

图4(c)为实施例提出的应用本发明所设计的ANFIS网络对MIT-BIH睡眠数据库中slp01样本进行分类学习前参数ratio(α)的初始隶属函数图；其中，横坐标为输入参数ratio(α)的范围，纵坐标为0～1的无量纲隶属度；

图4(d)为实施例提出的应用本发明所设计的ANFIS网络对MIT-BIH睡眠数据库中slp01样本进行分类学习后参数ratio(α)的最终隶属度函数图；其中，横坐标为输入参数ratio(α)的范围，纵坐标为0～1的无量纲隶属度；

图5为实施例提出的眼电伪迹去除实验所用单通道脑电信号示意图；

图6(a)为实施例提出的单通道脑电信号眼电伪迹去除效果对比图；

图6(b)为实施例提出的采集到的含有眼电伪迹的原始脑电信号效果图；

图7(a)为实施例提出的MIT-BIH数据库中的slp01样本参考睡眠状态示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图7(b)为实施例提出的采用ANFIS分类器对MIT-BIH睡眠数据库中slp01样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图7(c)为实施例提出的采用PLS分类器对MIT-BIH睡眠数据库中slp01样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图7(d)为实施例提出的采用LS-SVM分类器对MIT-BIH睡眠数据库中slp01样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图7(e)为实施例提出的样本为slp01的ANFIS分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图7(f)为实施例提出的样本为slp01的PLS分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图7(g)为实施例提出的样本为slp01的LS-SVM分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图8(a)为实施例提出的MIT-BIH数据库中的slp02样本参考睡眠状态示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图8(b)为实施例提出的采用ANFIS分类器对MIT-BIH睡眠数据库中slp02样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图8(c)为实施例提出的采用PLS分类器对MIT-BIH睡眠数据库中slp02样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图8(d)为实施例提出的采用LS-SVM分类器对MIT-BIH睡眠数据库中slp02样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图8(e)为实施例提出的样本为slp02的ANFIS分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图8(f)为实施例提出的样本为slp02的PLS分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图8(g)为实施例提出的样本为slp02的LS-SVM分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图9(a)为实施例提出的MIT-BIH数据库中的slp04样本参考睡眠状态示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图9(b)为实施例提出的采用ANFIS分类器对MIT-BIH睡眠数据库中slp04样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图9(c)为实施例提出的采用PLS分类器对MIT-BIH睡眠数据库中slp04样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图9(d)为实施例提出的采用LS-SVM分类器对MIT-BIH睡眠数据库中slp04样本的分类效果示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图9(e)为实施例提出的样本为slp04的ANFIS分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图9(f)为实施例提出的样本为slp04的PLS分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图9(g)为实施例提出的样本为slp04的LS-SVM分类器输出误差示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图10为具体实施方式八提出的ANFIS建模输入参数个数与分类误差关系图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图11(a)为具体实施方式九提出的MIT-BIH睡眠数据库中slp01样本参考标准睡眠状态曲线示意图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图11(b)为具体实施方式九提出的MIT-BIH睡眠数据库中slp01样本复杂度complex曲线图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图11(c)为具体实施方式九提出的MIT-BIH睡眠数据库中slp01样本希尔伯特瞬时频率均值hilbert的曲线图。其中，横坐标为样本个数，纵坐标为睡眠等级；

图11(d)为具体实施方式九提出的MIT-BIH睡眠数据库中slp01样本最大李雅普诺夫指数lyapunov的曲线图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图11(e)为具体实施方式九提出的MIT-BIH睡眠数据库中slp01样本δ频带能量比ratio(δ)的曲线图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图11(f)为具体实施方式九提出的MIT-BIH睡眠数据库中slp01样本θ频带能量比ratio(θ)的曲线图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图11(g)为具体实施方式九提出的MIT-BIH睡眠数据库中slp01样本α频带能量比ratio(α)的曲线图；其中，横坐标为样本个数，纵坐标为睡眠等级；

图11(h)为具体实施方式九提出的MIT-BIH睡眠数据库中slp01样本β频带能量比ratio(β)的曲线图。其中，横坐标为样本个数，纵坐标为睡眠等级。

具体实施方式

具体实施方式一：本实施方式的一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，具体是按照以下步骤制备的：

步骤一、对采集到的单通道脑电信号X(n)进行小波变换，得到M个小波系数；其中，M个小波系数分为两类：不包含眼电伪迹的小波系数为P个以及含眼电伪迹的小波系数为M-P个；

步骤二、对于不包含眼电伪迹的P个小波系数，直接作为纯净的脑电信号；

步骤三、含有眼电伪迹的M-P个小波系数W(i)经过经验模态分解后表示为N个内蕴模态函数即IMF(Intrinsic Mode Function)分量的个数与残差之和：

       $W\left(i\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{C}_j+r---\left(1\right)$

其中：W(i)为第i个含有眼电伪迹的小波系数，C_j为该小波系数的第j个内蕴模态函数，r为残差，i＝1,2,3，…，M-P，j＝1,2,3，…，N；

步骤四、将步骤三得到的每个含有眼电伪迹的小波系数的所有内蕴模态函数C_j以及残差r输入到独立分量分析系统ICA(Independent component analysis，)中，得到每个含有眼电伪迹的小波系数的脑电分量以及眼电分量；

步骤五、将步骤四得到的每个小波系数的眼电分量置零，并采用ICA逆运算重构纯净的脑电信号；

步骤六、将步骤二和步骤五中的得到的纯净脑电信号求和，得到去除眼电伪迹后的纯净脑电信号X(n)如图1；

步骤七、对步骤六得到的纯净脑电信号X(n)提取7个特征参数即4个频域参数、2个非线性参数以及1个时频参数；其中，4个频域参数包括δ，θ，α和β节律的频带能量比即ratio(δ)、ratio(θ)、ratio(α)和ratio(β)，2个非线性参数包括复杂度complexity简称c(n)和最大Lyapunov指数简称Lyapunov，时频参数为希尔伯特黄变换(HHT)瞬时频率均值简称hilber；

步骤八、将步骤七提取的3个频域参数、2个非线性参数以及1个时频参数到的特征参数输入到自适应模糊神经推理系统中，得到睡眠状态分期指数(如图2)；其中，3个频域参数为在步骤七提取的4个频域参数选择3个频域参数；即完成了一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法。

本实施方式效果：

本实施方式的目的是为了解决基于单通道脑电信号的睡眠状态分期过程中已有方法无法有效去除眼电伪迹的干扰，且不能充分利用脑电信号中蕴含的睡眠信息以及分类效果不佳的缺点，提出了一种结合小波变换，经验模态分解以及独立成分分析的单通道脑电信号眼电伪迹去除方法，同时还提出了一种综合脑电信号频域分析、非线性分析以及时频分析的特征参数提取方法，结合自适应模糊神经推理系统(ANFIS)分类器实现睡眠分期。从而解决不能充分挖掘脑电信号中蕴含的睡眠信息的问题。

在单通道脑电信号眼电伪迹去除的研究方面，本实施方式提出了基于小波变换，经验模态分解和独立分量分析的WT-EMD-ICA方法，有效的去除了脑电信号中的眼电伪迹，而且使得有用的脑电信息得以充分保持；在与睡眠状态分期相关的脑电信号特征参数提取方面，目前常用的脑电特征提取方法只注重其某一特征的分析，不能充分挖掘脑电信号中蕴含的睡眠信息，本实施方式综合脑电信号频域分析、非线性分析以及时频分析的特征参数提取方法，提取睡眠脑电信号中4个节律的‘频带能量比’参数，2个非线性参数以及1个时频参数，进一步提取了脑电信号中与睡眠分期相关的信息同时，本实施方式采用具有较好知识表达和学习能力的自适应模糊神经推理系统，得到睡眠状态分期指标具有计算速度快，准确度高等优点。如图7～9所示，相比于偏最小二乘分类器(PLS)和最小二乘支持向量机分类器(LS-SVM)，自适应模糊神经推理系统分类器(ANFIS)的睡眠分期正确率明显提高。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤七中提取的频域参数ratio(δ)、ratio(θ)、ratio(α)和ratio(β)具体为：

按式(2)，(3)，(4)和(5)计算纯净脑电信号X(n)中各频带信号的频带能量比即4个频域参数的特征参数分别为：

ratio(δ)＝E(δ)/E_all  (2)

其中，ratio(δ)为纯净脑电信号X(n)中δ频带的频带能量比特征参数；

ratio(θ)＝E(θ)/E_all  (3)

其中，ratio(θ)为纯净脑电信号X(n)中θ频带的频带能量比特征参数；

ratio(α)＝E(α)/E_all  (4)

其中，ratio(α)为纯净脑电信号X(n)中α频带的频带能量比特征参数；

ratio(β)＝E(β)/E_all  (5)

其中，ratio(β)为纯净脑电信号X(n)中β频带的频带能量比特征参数。其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：E(δ)、E(θ)、E(α)和E(β)具体为：

对去除眼电伪迹的纯净脑电信号进行频域分析，提取4个节律脑电信号的‘频带能量比’特征参数；

(1)、对纯净脑电信号X(n)按式(6)进行离散傅里叶变换，得到信号X(n)的功率频谱P(k)：

       $P\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}X\left(n\right)e^{-I\frac{2\mathrm{\π}}{N}\mathrm{kn}},k=\mathrm{0,1},...,N-1---\left(6\right)$

其中，N为纯净脑电信号X(n)的点数；P(k)为纯净脑电信号X(n)的频谱；n＝0,1,...,N-1；

(2)、按式(7)计算脑电信号X(n)中频带的信号能量E(δ)：

       $E\left(\mathrm{\δ}\right)=\underset{1}{\overset{4}{\∫}}P{\left(k\right)}^2\mathrm{dk}---\left(7\right)$

其中，E(δ)为纯净脑电信号X(n)中δ频带的能量，δ的频率范围为1～4Hz；

按式(8)计算脑电信号X(n)中频带的信号能量E(θ)：

       $E\left(\mathrm{\θ}\right)=\underset{4}{\overset{8}{\∫}}P{\left(k\right)}^2\mathrm{dk}---\left(8\right)$

其中，E(θ)为纯净脑电信号X(n)中θ频带的能量，θ的频率范围为4～8Hz；

按式(9)计算脑电信号X(n)中频带的信号能量E(α)：

       $E\left(\mathrm{\α}\right)=\underset{8}{\overset{13}{\∫}}P{\left(k\right)}^2\mathrm{dk}---\left(9\right)$

其中，E(α)为纯净脑电信号X(n)中α频带的能量，α的频率范围为8～13Hz；

按式(10)算脑电信号X(n)中频带的信号能量E(β)：

       $E\left(\mathrm{\β}\right)=\underset{13}{\overset{30}{\∫}}P{\left(k\right)}^2\mathrm{dk}---\left(10\right)$

其中，E(β)为纯净脑电信号X(n)中β频带的能量，β的频率范围为13～30Hz。其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：E_all具体为：

按式(11)算纯净脑电信号X(n)中各频带的信号能量和E_all：

E_all＝E(δ)+E(θ)+E(α)+E(β)  (11)

其中，E_all为纯净脑电信号X(n)中各频带的信号能量和。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤七中提取2个非线性参数即复杂度complexity简称c(n)和最大Lyapunov指数简称Lyapunov具体为：

对去除眼电伪迹的纯净脑电信号X(n)进行非线性分析，提取复杂度和最大Lyapunov指数等2个非线性脑电特征参数；

(1)、计算纯净脑电信号X(n)的均值m，并利用均值m重构二进制序列S_i＝S₁,S₂,…,S_N，当X(n)>m时，S_i＝1，当X(n)<m时，S_i＝0；其中，n＝1,2,3，....N；

(2)、对二进制序列S_i进行子串划分；

(3)、统计符号*将字符串S分成的份数，即为复杂度c(n)的值；例如，对于(2)的序列“10110100”，通过复杂度计算，“*”将其分成了5个子序列，因此该序列的复杂度为5；

(4)、对c(n)都趋向于一定值，按照下式对复杂度c(n)进行归一化处理

c(n)＝c(n)/b(n)  (12)

其中，b(n)＝n/lgn；

(5)、采用C-C法计算纯净脑电信号X(n)的最优时间延迟τ和嵌入维数m，根据τ和m进行相空间重构，Y_n是X(n)的重构向量即Y_n＝(X(n),X(n+τ),....,X(n+(m-1)τ))∈R^m，其中，n＝1,2,…,N-(m-1)τ；

(6)、对纯净脑电信号X(n)进行离散傅里叶变换(FFT变换)，计算X(n)功率频谱 $\mathrm{Power}\left(k\right)={\left|\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}X\left(n\right)e^{-j\frac{2\mathrm{\&pi;}}{N}\mathrm{kn}}\right|}^2,k=\mathrm{0,1},...N-1,$ 其中字母与前面公式(6)完全一致，并记录Power中最大值所对应的序号num；令period(n)＝N/n,则平均周期P＝period(num)；

(7)、找Y_n相空间中每个点Y_j的最近邻点并限制短暂分离，即：

       $d_j\left(0\right)=\underset{\hat{j}}{\min }\left|\right|Y_j-Y_{\hat{j}}\left|\right|,$ 且， $|j-\hat{j}|>P---\left(13\right)$

其中d_j(0)为点Y_j的最近邻点的初始距离；

(8)、对相空间中每个点Y_j，计算出对经过i个离散时间步后的距离d_j(i)；

       $d_j\left(i\right)=|Y_{j+i}-Y_{\hat{j}+i}|,$ 其中， $i=\mathrm{1,2},...\min (M-j,M-\hat{j})---\left(14\right)$

其中，M＝N-(m-1)τ；

(9)、对第i个离散时间，求出所有j的ln d_j(i)平均y(i)，即

       $y\left(i\right)=\frac{1}{\mathrm{q\&Delta;t}}\underset{j=1}{\overset{q}{\mathrm{\&Sigma;}}}\ln d_j\left(i\right)---\left(15\right)$

其中，q是非零d_j(i)的数目，并用最小二乘法作出回归直线，回归直线的斜率就是最大Lyapunov指数；△t为相邻两次离散时间的间隔。其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：利用L-Z复杂度方法对二进制序列S_i，i＝1,2,…,N进行复杂度求解方法具体为：

(1)、令S＝S₁，S₂，…，S_i，Q＝S_i+1；将S和Q按照前后顺序组合到一起，得到一个新的字符串SQ；

(2)、定义SQv是SQ减去SQ中最后一个字符所得的字符串；

(3)、判断Q是否是SQv中的一个子串，如果Q是SQv的一个子串，令S不变，继续增长Q，令Q为S_i+1S_i+2；如果Q不是SQv的一个子串，则用*把SQ中S和Q前后分开，令S＝S₁，S₂，…，S_i，S_i+1，Q＝S_i+2重复步骤(1)～(3)么；直到Q＝S_N为止；

(4)、序列S_i，i＝1,2,…,N的复杂性定义为由符号*界定的S的子串数目；例如序列10110100

      

其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：步骤七中提取希尔伯特瞬时频率均值hilber具体过程：

对去除眼电伪迹的纯净脑电信号进行时频分析，提取希尔伯特瞬时频率均值特征参数；

(1)、对纯净脑电信号X(n)进行经验模态分解，分解后的序列利用公式(1)表示为N个内蕴模态函数与残差的和；

(2)、计算每个IMF分量的能量

       $E_{C_j}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_j^2\left(n\right)---\left(16\right)$

式中，C_j为第j个内蕴模态函数，n为纯净脑电信号X(n)的点数；

(3)、将IMF分量的能量按从大到小的顺序排列，选择能量累计贡献率大于80％的累计贡献率从大到小前f个IMF分量；

       $n_f=\underset{j=1}{\overset{f}{\mathrm{\&Sigma;}}}{E}_{C_j}/\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{E}_{C_j}---\left(17\right)$

(4)、对前f个IMF分量中每个选中的IMF分量按式(18)进行希尔伯特变换：

       $Y_j\left(n\right)=H\left\{C_j\left(n\right)\right\}=\frac{1}{\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{\&infin;}{C}_j\left(s\right)\frac{1}{t-s}\mathrm{ds}---\left(18\right)$

C_j(n)为第j个内蕴模态函数；

C_j(s)为第j个内蕴模态函数的积分形式，s为积分算子；

Y_j(n)为C_j(n)的Hilbert变换结果；

H{}表示Hilbert变换；

(5)、构造解析信号Z_j(n)；

       $Z_j\left(n\right)=C_j\left(n\right)+{\mathrm{IY}}_j\left(n\right)=a_j\left(n\right)e^{I{\mathrm{\&theta;}}_j\left(n\right)}---(19$

其中， $a_j\left(n\right)=\sqrt{C_j^2\left(n\right)+Y_j^2\left(n\right)},{\mathrm{\&theta;}}_j\left(n\right)=\arctan \frac{Y_j\left(n\right)}{C_j\left(n\right)};$

I为虚数，I²＝-1

(6)、计算C_j(n)的瞬时频率均值w_j(n)；

       $\mathrm{mean}\left[w_j\left(n\right)\right]=\mathrm{mean}\left[\frac{d{\mathrm{\&theta;}}_j\left(n\right)}{\mathrm{dt}}\right],n=\mathrm{1,2},...,N-1---\left(20\right)$

(7)、按能量权值计算所有选中的IMF分量总体的瞬时频率均值：

       $\mathrm{hilbert}=\underset{j\&Element;J}{\mathrm{\&Sigma;}}\mathrm{mean}\left[w_j\left(n\right)\right]\&times;(E_{C_j}/\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{E}_{C_k})---\left(21\right)$

其中，J为所有IMF分量集，参数hilbert即为序列X(n)的希尔伯特瞬时频率均值；

      为IMF分量集J中的j个IMF函数能量和；

      为所有N个IMF函数的能量和。其它步骤及参数与具体实施方式一至六之一相同。

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：步骤八中将步骤七提取3个频域参数具体为：由于δ，θ，α和β节律的频带能量比之和始终为1即ratio(δ)+ratio(θ)+ratio(α)+ratio(β)＝1，因此对δ，θ，α和β等4节律的频带能量比参数进行简化，采用δ，α和β等3个节律(实验验证的结果由图10得出结论得出采用δ，α和β等3个节律)的频带能量比即ratio(δ)、ratio(α)和ratio(β)表示。其它步骤及参数与具体实施方式一至七之一相同。

具体实施方式九：本实施方式与具体实施方式一至八之一不同的是：步骤八中将步骤七提取的3个频域参数、2个非线性参数以及1个时频参数到的特征参数输入到自适应模糊神经推理系统中，得到睡眠状态分期指数具体过程为如图3：

(1)、由于ratio(δ)、ratio(θ)、ratio(α)和ratio(β)的和始终为1，因此可以对输入参数的个数以及ANFIS网络结构进行精简；本发明将步骤二提取出的7个特征参数中的6个参数输入到自适应模糊神经推理系统网络的第一层；将6个输入参数分别通过式(22)的高斯型隶属度函数进行模糊化处理：

       $\mathrm{\&mu;}\left(x\right)=\exp (-\frac{{(x-c)}^2}{2{\mathrm{\&sigma;}}^2})---\left(22\right)$

其中，c为隶属度函数的中心，σ为隶属度函数的宽度；x为隶属度函数的输入参数，μ(x)为隶属度函数；

(2)、6个输入参数分别为：3个脑电频带能量比特征参数ratio(δ)、ratio(α)和ratio(β)，2个非线性参数即复杂度complex和最大李雅普诺夫指数lyapunov和1个时频参数，希尔伯特黄变换(HHT)瞬时频率均值hilber；

ratio(δ)参数与参考睡眠状态的相关系数在所有特征参数中最大(这是公知常识么)，因此设定ratio(δ)参数有4个隶属度函数，复杂度、ratio(α)和ratio(β)等3个参数与睡眠状态也具有高度相关性，因此设定复杂度complex、ratio(α)和ratio(β)的分别有3个隶属度函数，设定将希尔伯特黄变换(HHT)瞬时频率均值hilbert及最大李雅普诺夫指数lyapunov分别有1个隶属度函数，(将希尔伯特黄变换(HHT)瞬时频率均值hilbert及最大李雅普诺夫指数lyapunov与睡眠状态指数的相关系数较低)；图11(a)～(h)证明了此结论；因此自适应模糊神经推理系统网络第一层具有15个节点即15个隶属度函数；

自适应模糊神经推理系统网络的第一层含有15个隶属度函数，每个隶属度函数有2个参数，c和σ即共有15×2＝30个参数；

该网络的规则库有3×1×1×4×3×3＝108个规则组成；规则的形式为：

(1)如果complex为A₁且hilbert为B且lyapunov为C且ratio(δ)为D₁且ratio(α)为E₁且ratio(β)为F₁

则输出为：p₁ ¹×complex+p₁ ²×hilbert+p₁ ³×lyapunov+p₁ ⁴×ratio(δ)+p₁ ⁵×ratio(α)+p₁ ⁶×ratio(β)+r₁；

(2)如果complex为A₁且hilbert为B且lyapunov为C且ratio(δ)为D₁且ratio(α)为E₁且ratio(β)为F₂

则输出为：p₂ ¹×complex+p₂ ²×hilbert+p₂ ³×lyapunov+p₂ ⁴×ratio(δ)+p₂ ⁵×ratio(α)+p₂ ⁶×ratio(β)+r₂；

…

(108)如果complex为A₃且hilbert为B且lyapunov为C且ratio(δ)为D₄且ratio(α)为E₃且ratio(β)为F₃

则输出为：p₁₀₈ ¹×complex+p₁₀₈ ²×hilbert+p₁₀₈ ³×lyapunov+p₁₀₈ ⁴×ratio(δ)+p₁₀₈ ⁵×ratio(α)+p₁₀₈ ⁶×ratio(β)+r₁₀₈；

每条规则有7个参数，即规则参数共有108×7＝756个，整个网络共有786个需要调整的参数，需要调整的参数即为自适应模糊神经推理系统网络权值；

为便于表示采用O_t ^T表示自适应模糊神经推理系统网络第T层的第t个节点的输出；自适应模糊神经推理系统网络简称网络第一层为隶属度函数层，隶属度函数层中的每个节点用参数ratio(δ)的隶属度函数表示为：

       $O_{D_d}^1\left(\mathrm{ratio}\left(\mathrm{\&delta;}\right)\right)={\mathrm{\&mu;}}_{D_d}\left(\mathrm{ratio}\left(\mathrm{\&delta;}\right)\right),d=\mathrm{1,2,3,4}---\left(23\right)$

同理得其他五个参数的隶属度函数complexity，hilbert，Lyapunov，ratio(α)，ratio(β)隶属度函数分别为：

       $O_{A_{\mathrm{\&alpha;}}}^1\left(\mathrm{complex}\right)={\mathrm{\&mu;}}_{A_{\mathrm{\&alpha;}}}\left(\mathrm{complex}\right),a=\mathrm{1,2,3}$

       $O_{B_b}^1\left(\mathrm{hilbert}\right)={\mathrm{\&mu;}}_{B_b}\left(\mathrm{hilbert}\right),b=1$

       $O_{C_c}^1\left(\mathrm{lyapunov}\right)={\mathrm{\&mu;}}_{C_c}\left(\mathrm{lyapunov}\right),c=\mathrm{1,2,3}$

       $O_{E_e}^1\left(\mathrm{ratio}\left(\mathrm{\&alpha;}\right)\right)={\mathrm{\&mu;}}_{E_e}\left(\mathrm{ratio}\left(\mathrm{\&alpha;}\right)\right),e=\mathrm{1,2,3}$

       $O_{F_f}^1\left(\mathrm{ratio}\left(\mathrm{\&beta;}\right)\right)={\mathrm{\&mu;}}_{F_f}\left(\mathrm{ratio}\left(\mathrm{\&beta;}\right)\right),f=\mathrm{1,2,3}$

(2)、按式(24)计算网络第二层强度释放层的输出O_h：

       $O_h^2=w_h={\mathrm{\&mu;}}_{A_a}{\mathrm{\&mu;}}_{B_b}{\mathrm{\&mu;}}_{C_c}{\mathrm{\&mu;}}_{D_d}{\mathrm{\&mu;}}_{E_e}{\mathrm{\&mu;}}_{F_f}---\left(24\right)$

其中，h为自适应模糊神经推理系统网络第2层的第h个输出，且h＝(a-1)×4×3²+(d-1)×3²+(e-1)×3+f；

a为隶属度函数A的序号，a＝1，2，3，4

d为隶属度函数B的序号，d＝1，2，3

e为隶属度函数E的序号，e＝1，2，3

f为隶属度函数F的序号，f＝1，2，3

      为第A_a个隶属度函数的输出，为第B_b个隶属度函数的输出，

      为第C_c个隶属度函数的输出，为第D_d个隶属度函数的输出,

      为第E_e个隶属度函数的输出，为第F_f个隶属度函数的输出；

(3)、按式(25)计算网络第三层强度归一化层的输出；

       $O_i^3={\stackrel{\&OverBar;}{w}}_i=\frac{w_i}{\mathrm{\&Sigma;}{w}_g},g=\mathrm{1,2},3...108---\left(25\right)$

(4)、按式(26)计算网络第四层模糊规则层的输出

       $\begin{array}{c}{O}_i^4={\stackrel{\&OverBar;}{w}}_i(p_i^1\&times;\mathrm{complex}+p_i^2\&times;\mathrm{hilbert}+p_i^3\&times;\mathrm{lyapunov}+\\ p_i^4\&times;\mathrm{ratio}\left(\mathrm{\&delta;}\right)+p_i^5\&times;\mathrm{ratio}\left(\mathrm{\&alpha;}\right)+p_i^6\&times;\mathrm{ratio}\left(\mathrm{\&beta;}\right)+r_i)\end{array}---\left(26\right)$

(5)、按式(27)计算自适应模糊神经推理系统网络第五层的输出，因自适应模糊神经推理系统网络的第五层只有1个输出；即为睡眠状态分期结果

       $O_1^5=\underset{c=1}{\overset{108}{\mathrm{\&Sigma;}}}{O}_c^4---\left(27\right)$

其它步骤及参数与具体实施方式一至八之一相同。

采用以下实施例验证本发明的有益效果：

实施例一：

本实施例一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，具体是按照以下步骤制备的：

1.对本发明步骤一中所述WT-EMD-ICA单通道脑电信号眼电伪迹去除方法的实验验证，图5和图6(a)和(b)为相关的实验效果图。

本实验所用脑电采集设备为哈尔滨工业大学智能测试集信息处理技术研究所研制的便携式三电极脑电采集设备。采集部位为实验对象的前额。在三个电极中，其中一个作为参考电极，另外两个电极的电势差为前额脑电的采集数据。脑电采样率为512Hz。实验过程中，实验对象首先保持安静和放松状态，眼睛保持闭眼状态，一段时间后，刻意快速眨眼若干次，然后恢复初始状态。采集整个过程中的前额脑电数据。图5为其中一段脑电信号。图5中，(a)处为由眨眼造成的眼电伪迹，(b)处为未受眼电伪迹干扰的脑电信号。

采用‘sym7’母小波对该段脑电信号进行7层小波分解。应用WT-EMD-ICA算法对其进行眼电伪迹去除。图5中眼电伪迹去除前后的效果对比如图6所示。图6中，眼电伪迹被明显去除。同时，本实验选择了图5中信号P1：1.37～1.95s，P2：4.49～4.98s，P3：5.47～6.05s以及P4：6.64～7.42s不含眼电伪迹的信号部分，P1～P4部分信号经过WT-EMD-ICA算法处理前后的相关系数分别为：0.88，0.93，0.96和0.97。由此证明，经过该算法，眼电伪迹信号被明显消除，而且脑电信号中包含的有用信息得以较好的保持，为下一步睡眠脑电信号特征提取与分类提供最优的原始数据。

2.对本发明步骤二和步骤三中所述脑电信号特征参数提取与分类方法的实验验证。

对来自MIT-BIH多导睡眠数据库的睡眠脑电进行睡眠状态分期，结果证明本发明方法有效可行，大大提高了睡眠状态分期精度，为睡眠障碍疾病的诊断提供了有效的数据支持。

MIT-BIH数据库包含18组睡眠脑电数据。脑电信号采样率为250Hz。本专利采用slp01、slp02及slp04等3个样本作为分析对象。其中，slp01和slp02样本分别包括slp01a和slp01b，slp02a和slp02b两部分。slp01样本睡眠脑电采集时间为5小时，数据长度为4.5×10⁶。该数据库每30s给出一个睡眠状态指数，即每7500个点脑电数据对应一个睡眠状态，共600组数据。slp02样本采集时间为5小时15分钟，数据长度为4.725×10⁶，共630组数据。slp04样本采集时间为6小时，脑电数据长度为5.4×10⁶，共720组数据。数据库参考睡眠状态包括：清醒期(W)、睡眠1期(S1)、快速眼动期(R)、睡眠2期(S2)、睡眠3期(S3)以及睡眠4期(S4)等6个状态，分别用0～5这6个数字表示，数字越大表示睡眠程度越深。其中，slp02号样本没有S3和S4这两个睡眠分期的数据，slp04号样本没有S4睡眠分期的数据。

附图4(a)～(d)为应用本发明所设计的ANFIS网络对slp01样本进行分类学习前后参数ratio(δ)和ratio(α)隶属度变化对比图。经过学习后，各参数的隶属度函数的形状发生了较大的改变。学习后的隶属度函数是对输入参数最佳的模糊化处理函数。

采用PLS以及LS-SVM分类器与本发明所提方法进行比较，附图7(a)～(g)、附图8(a)～(g)和附图9(a)～(g)分别为采用不同方法对slp01、slp02以及slp04号样本的分类效果对比图。

采用分类标准偏差(standard error of classification，SEC)以及分类正确率作为分类效果的评价标准。

       $\mathrm{SEC}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(y_{i,\mathrm{actual}}-y_{i,\mathrm{predicted}})}^2}{n-1}}---\left(28\right)$

其中，y_i,actual为第i个样本的参考标准睡眠状态，y_i,predicted为所用分类器对第i个样本的输出，n为样本数目。

分类正确率分为总体分类正确率和各睡眠状态分类正确率。总体分类正确率为个样本的正确分类数目与样本总数的比值。各睡眠状态分类正确率为各睡眠状态的正确分类数目与各睡眠状态数目的比值。表1～3为各样本分类效果数据表。表中NaN代表该样本无此状态数据。

表1 slp01样本特征分类对比

      

表2 slp02样本特征分类对比

      

表3 slp04样本特征分类对比

      

从附图7中图a),c),e)中可以看出，三种分类器都可以实现通过脑电信号对睡眠状态进行分期处理，而且基本都能够反映各个样本的睡眠状态变化趋势。但从图b),d),f)中可以明显看出，三者的分类效果差别还是较大的。PLS分类器的误差相比其他两种分类器明显较高，而且波动较大。究其原因，PLS分类器是基于线性回归的原理。而脑电信号具有严重的非线性特性，因此采用线性原理方法进行分类，分类误差一定会偏大。相比之下，ANFIS和LS-SVM采用非线性原理进行分类，三个样本的分类标准偏差SEC都小于0.60。

ANFIS方法融合了神经网络和模糊推理系统的优点，相比LS-SVM分类器具有更高的学习拟合能力。由表1～3可知，采用ANFIS分类器的三个样本分类标准偏差SEC都在0.30左右，而LS-SVM分类器的分类标准偏差在0.50左右。对于slp02样本，采用ANFIS分类器的总体分类正确率为96.03％，而采用LS-SVM分类器的总体分类正确率仅为77.46％。三个样本中，LS-SVM分类器的最高总体分类正确率为slp01样本的82.83％。

综上所述，本发明通过引入脑电信号的非线性参数结合频域参数进行特征提取，最后采用ANFIS分类器，相比已有方法，大大提高了睡眠脑电分期的分类精度，对于研究睡眠对神经功能的影响、睡眠障碍、睡眠和其他医学状态紊乱的关系具有十分重要的意义。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (9)

1.一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法具体是按照以下步骤进行的：

步骤一、对采集到的单通道脑电信号X(n)进行小波变换，得到M个小波系数；其中，M个小波系数分为两类：不包含眼电伪迹的小波系数为P个以及含眼电伪迹的小波系数为M-P个；

步骤二、对于不包含眼电伪迹的P个小波系数，直接作为纯净的脑电信号；

步骤三、含有眼电伪迹的M-P个小波系数W(i)经过经验模态分解后表示为N个内蕴模态函数即IMF分量的个数与残差之和：

其中：W(i)为第i个含有眼电伪迹的小波系数，C_j为该小波系数的第j个内蕴模态函数，r为残差，i＝1,2,3，…，M-P，j＝1,2,3，…，N；

步骤四、将步骤三得到的每个含有眼电伪迹的小波系数的所有内蕴模态函数C_j以及残差r输入到独立分量分析系统ICA中，得到每个含有眼电伪迹的小波系数的脑电分量以及眼电分量；

步骤五、将步骤四得到的每个小波系数的眼电分量置零，并采用ICA逆运算重构纯净的脑电信号；

步骤六、将步骤二和步骤五中的得到的纯净脑电信号求和，得到去除眼电伪迹后的纯净脑电信号X(n)；

步骤七、对步骤六得到的纯净脑电信号X(n)提取7个特征参数即4个频域参数、2个非线性参数以及1个时频参数；其中，4个频域参数包括δ，θ，α和β节律的频带能量比即ratio(δ)、ratio(θ)、ratio(α)和ratio(β)，2个非线性参数包括复杂度complexity简称c(n)和最大Lyapunov指数简称Lyapunov，时频参数为希尔伯特黄变换瞬时频率均值简称hilber；

步骤八、将步骤七提取的3个频域参数、2个非线性参数以及1个时频参数到的特征参数输入到自适应模糊神经推理系统中，得到睡眠状态分期指数；其中，3个频域参数为在步骤七提取的4个频域参数选择3个频域参数；即完成了一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法。

2.根据权利要求1所述一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于：步骤七中提取的频域参数ratio(δ)、ratio(θ)、ratio(α)和ratio(β)具体为：

按式(2)，(3)，(4)和(5)计算纯净脑电信号X(n)中各频带信号的频带能量比即4个频域参数的特征参数分别为：

ratio(δ)＝E(δ)/E_all            (2) 

其中，ratio(δ)为纯净脑电信号X(n)中δ频带的频带能量比特征参数；

ratio(θ)＝E(θ)/E_all      (3) 

其中，ratio(θ)为纯净脑电信号X(n)中θ频带的频带能量比特征参数；

ratio(α)＝E(α)/E_all          (4) 

其中，ratio(α)为纯净脑电信号X(n)中α频带的频带能量比特征参数；

ratio(β)＝E(β)/E_all          (5) 

其中，ratio(β)为纯净脑电信号X(n)中β频带的频带能量比特征参数。

3.根据权利要求2所述一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于：E(δ)、E(θ)、E(α)和E(β)具体为：

(1)、对纯净脑电信号X(n)按式(6)进行离散傅里叶变换，得到信号X(n)的功率频谱P(k)：

其中，N为纯净脑电信号X(n)的点数；P(k)为纯净脑电信号X(n)的频谱；n＝0,1,...,N-1；

(2)、按式(7)计算脑电信号X(n)中频带的信号能量E(δ)：

其中，E(δ)为纯净脑电信号X(n)中δ频带的能量，δ的频率范围为1～4Hz；

按式(8)计算脑电信号X(n)中频带的信号能量E(θ)：

其中，E(θ)为纯净脑电信号X(n)中θ频带的能量，θ的频率范围为4～8Hz；

按式(9)计算脑电信号X(n)中频带的信号能量E(α)：

其中，E(α)为纯净脑电信号X(n)中α频带的能量，α的频率范围为8～13Hz；

按式(10)计算脑电信号X(n)中频带的信号能量E(β)：

其中，E(β)为纯净脑电信号X(n)中β频带的能量，β的频率范围为13～30Hz。

4.根据权利要求2所述一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于：E_all具体为：

按式(11)计算纯净脑电信号X(n)中各频带的信号能量和E_all：

E_all＝E(δ)+E(θ)+E(α)+E(β)       (11) 

其中，E_all为纯净脑电信号X(n)中各频带的信号能量和。

5.根据权利要求1所述一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于：步骤七中提取2个非线性参数即复杂度complexity简称c(n)和最大Lyapunov指数简称Lyapunov具体为：

(1)、计算纯净脑电信号X(n)的均值m，并利用均值m重构二进制序列S_i＝S₁,S₂,…,S_N，当X(n)>m时，S_i＝1，当X(n)<m时，S_i＝0；其中，n＝1,2,3，….N；

(2)、对二进制序列S_i进行子串划分；

(3)、统计符号*将字符串S分成的份数，即为复杂度c(n)的值；

(4)、按照下式对复杂度c(n)进行归一化处理

c(n)＝c(n)/b(n)         (12) 

其中，b(n)＝n/lgn；

(5)、采用C-C法计算纯净脑电信号X(n)的最优时间延迟τ和嵌入维数m，根据τ和m进行相空间重构，Y_n是X(n)的重构向量即Y_n＝(X(n),X(n+τ),....,X(n+(m-1)τ))∈R^m，其中，n＝1,2,…,N-(m-1)τ；

(6)、对纯净脑电信号X(n)进行离散傅里叶变换，计算X(n)功率频谱 k＝0,1,...N-1，其中字母与前面公式(6)完全一致，并记录Power中最大值所对应的序号num；令period(n)＝N/n,则平均周期P＝period(num)；

(7)、找Y_n相空间中每个点Y_j的最近邻点并限制短暂分离，即：

且，

其中d_j(0)为点Y_j的最近邻点的初始距离；

(8)、对相空间中每个点Y_j，计算出对经过i个离散时间步后的距离d_j(i)；

其中，

其中，M＝N-(m-1)τ；

(9)、对第i个离散时间，求出所有j的ln d_j(i)平均y(i)，即

其中，q是非零d_j(i)的数目，并用最小二乘法作出回归直线，回归直线的斜率就是最大Lyapunov指数_；△t为相邻两次离散时间的间隔。

6.根据权利要求5所述一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于：利用L-Z复杂度方法对二进制序列S_i，i＝1,2,…,N进行复杂度求解方法具体为：

(1)、令S＝S₁，S₂，…，S_i，Q＝S_i+1；将S和Q按照前后顺序组合到一起，得到一个新的字符串SQ；

(2)、定义SQv是SQ减去SQ中最后一个字符所得的字符串；

(3)、判断Q是否是SQv中的一个子串，如果Q是SQv的一个子串，令S不变，继续增长Q，令Q为S_i+1S_i+2；如果Q不是SQv的一个子串，则用*把SQ中S和Q前后分开，令S＝S₁，S₂，…，S_i，S_i+1，Q＝S_i+2重复步骤(1)～(3)么；直到Q＝S_N为止；

(4)、序列S_i，i＝1,2,…,N的复杂性定义为由符号*界定的S的子串数目。

7.根据权利要求1所述一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于：步骤七中提取希尔伯特瞬时频率均值hilber具体过程：

(1)、对纯净脑电信号X(n)进行经验模态分解，分解后的序列利用公式(1)表示为N个内蕴模态函数与残差的和；

(2)、计算每个IMF分量的能量

式中，C_j为第j个内蕴模态函数，n为纯净脑电信号X(n)的点数；

(3)、将IMF分量的能量按从大到小的顺序排列，选择能量累计贡献率大于80％的累计贡献率从大到小前f个IMF分量；

(4)、对前f个IMF分量中每个选中的IMF分量按式(18)进行希尔伯特变换：

C_j(n)为第j个内蕴模态函数；

C_j(s)为第j个内蕴模态函数的积分形式，s为积分算子；

Y_j(n)为C_j(n)的Hilbert变换结果；

H{}表示Hilbert变换；

(5)、构造解析信号Z_j(n)；

其中，

I为虚数，I²＝-1

(6)、计算C_j(n)的瞬时频率均值w_j(n)；

(7)、按能量权值计算所有选中的IMF分量总体的瞬时频率均值：

其中，J为所有IMF分量集，参数hilbert即为序列X(n)的希尔伯特瞬时频率均值；

为IMF分量集J中的j个IMF函数能量和；

为所有N个IMF函数的能量和。

8.根据权利要求1所述一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于：步骤八中将步骤七提取3个频域参数具体为：由于δ，θ，α和β节律的频带能量比之和始终为1即ratio(δ)+ratio(θ)+ratio(α)+ratio(β)＝1，因此对δ，θ，α和β等4节律的频带能量比参数进行简化，采用δ，α和β等3个节律的频带能量比即ratio(δ)、ratio(α)和ratio(β)表示。

9.根据权利要求1所述一种基于单通道脑电信号眼电伪迹去除的睡眠分期方法，其特征在于：步骤八中将步骤七提取的3个频域参数、2个非线性参数以及1个时频参数到的特征参数输入到自适应模糊神经推理系统中，得到睡眠状态分期指数具体过程为：

(1)、将6个输入参数分别通过式(22)的高斯型隶属度函数进行模糊化处理：

其中，c为隶属度函数的中心，σ为隶属度函数的宽度；x为隶属度函数的输入参数，μ(x)为隶属度函数；

(2)、6个输入参数分别为：3个脑电频带能量比特征参数ratio(δ)、ratio(α)和ratio(β)，2个非线性参数即复杂度complex和最大李雅普诺夫指数lyapunov和1个时频参数，希尔伯特黄变换(HHT)瞬时频率均值hilber；

自适应模糊神经推理系统网络的第一层含有15个隶属度函数，每个隶属度函数有2个参数，c和σ即共有15×2＝30个参数；

表示采用O_t ^T表示自适应模糊神经推理系统网络第T层的第t个节点的输出；自适应模糊神经推理系统网络简称网络第一层为隶属度函数层，隶属度函数层中的每个节点用参数ratio(δ)的隶属度函数表示为：

同理得其他五个参数的隶属度函数complexity，hilbert，Lyapunov，ratio(α)，ratio(β)隶属度函数分别为：

(2)、按式(24)计算网络第二层强度释放层的输出O_h：

其中，h为自适应模糊神经推理系统网络第2层的第h个输出，且h＝(a-1)×4×3²+(d-1)×3²+(e-1)×3+f；

a为隶属度函数A的序号，a＝1，2，3，4

d为隶属度函数B的序号，d＝1，2，3

e为隶属度函数E的序号，e＝1，2，3

f为隶属度函数F的序号，f＝1，2，3

为第A_a个隶属度函数的输出，为第B_b个隶属度函数的输出，

为第C_c个隶属度函数的输出，为第D_d个隶属度函数的输出,

为第E_e个隶属度函数的输出，为第F_f个隶属度函数的输出；

(3)、按式(25)计算网络第三层强度归一化层的输出；

(4)、按式(26)计算网络第四层模糊规则层的输出

(5)、按式(27)计算自适应模糊神经推理系统网络第五层的输出，因自适应模糊神经推理系统网络的第五层只有1个输出；即为睡眠状态分期结果

。