EMD和CSP融合功率谱密度脑电特征提取方法与流程

本发明涉及emd和csp融合功率谱密度脑电特征提取方法，属于智能信息处理技术领域。

背景技术：

脑-机接口(braincomputerinterface，bci)是涉及众多学科和知识领域的控制技术。二十世纪九十年代年提出了脑-机接口的概念：它是在人或动物脑与外部设备间建立的直接传递信息的通路，通过采集与提取大脑产生的脑电信号来识别人的思想。脑机接口技术的发展可以帮助有运动障碍的病人提高自理能力和生活质量，可以通过大脑控制外部辅助设备如计算机、语音合成器、辅助应用和神经假肢等来加强他们与外界环境的交流和交互。脑机接口技术还被广泛地应用于其他一些领域，如游戏应用和导航。

经过多年的研究脑机接口技术主要包含信号采集、预处理、特征提取、特征分类和接口设备控制等五个步骤。特征提取可以识别不同想象运动的脑电信号的判别信息，因此在脑机接口的研究领域备受关注。由于eeg信号的非线性和非平稳性，传统的时间-频率方法：短时傅里叶变换(short-timefouriertransform,stft)，小波变换(wavelettransform,wt)等不可能同时得到时间-频率的良好分辨率。近年，希尔伯特黄变换(hilbert-huangtransform,hht)作为另一种时间-频率分析法已经变得越来越流行，同时它也十分适合分析非线性和非平稳信号。原信号经过经验模式分解(empiricalmodedecomposition,emd)被分解为一系列固有模态函数(intrinsicmodefunctions,imfs)，随后对每个固有模态函数进行希尔伯特黄变换，求其相应的能量谱和边际谱，作为特征进行分类。hht不涉及海森堡不确定性原理可以获得时域和频域的高分辨率。目前它被应用于多个领域，如雷达探测，地震信号和生物医学信号等。

在此基础上，由于eeg信号低空间分辨率，eeg信号构成的bci系统需要进行有效的空间滤波，从而确保从受试的相关脑域中提取特征信息。常用的算法有：共空间模式(commonspatialpattern,csp)、独立主成分分析(independentcomponentanalysis,ica)和共域空间谱模式(cssp)、滤波器csp(filterbankcommonspatialpattern,fbcsp)、判别滤波csp(discriminantfilteringcommonspatialpattern,dfbcsp)等多种csp的改进算法。

然而，传统的csp需要大量的输入通道，同时缺乏频率信息。csp算法主要是通过对脑电源的空间信息的判定来实现不同分类间的差异极大化，算法本质是利用代数上矩阵同时对角化的理论，寻找一组空间滤波器，脑电信号通过这组滤波器的投影获得较明显的特征向量。由于脑电信号是极复杂的非线性，不平稳信号，采集到的信号中，不止有脑电信号，还有其他的杂质信号，即使在脑电信号的频带范围内依旧存在其他的混叠信号，这些杂质与混叠信号将会影响特征提取的有效性，而单纯的csp滤波，只是通过空间信息来判定脑电信号的有效性，存在一定的不足，可能不会将杂质与混叠信号从频段中去除干净，如果许多与运动想象无关的频率信号混在其中，严重影响了特征向量的有效性。而本发明能够很好的解决上面的问题。

公开于该背景技术部分的信息仅仅旨在增加对本发明的总体背景的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域一般技术人员所公知的现有技术。

技术实现要素：

本发明的目的在于提供一种emd和csp融合功率谱密度脑电特征提取方法，从而克服上述现有技术中的缺陷。

为实现上述目的，本发明提供了一种emd和csp融合功率谱密度脑电特征提取方法，该方法是基于如下内容：步骤1，基于经验模式分解的eeg信号处理。步骤2，根据相关系数的计算分析筛选固有模态函数，构成新的信号矩阵。步骤3，计算信号矩阵的功率谱密度值对信号进行优化，再进行公共空间模式分解，解决csp多输入、缺频域信息的问题。步骤4，支持向量机分类。本发明首先对预处理后信号进行经验模式分解得到一系列固有模态函数，通过计算每次实验原始脑电信号与各阶imf分量之间的相关系数，并计算所有实验得出的相关系数的绝对值的平均数，选择具有较大相关系数绝对值平均数的imf,计算其功率谱密度作为特征，经csp投影映射再提取相应的特征向量，并用支持向量机进行特征选择得到分类结果。

本发明进一步限定的技术方案为：

优选地，上述技术方案中，该方法包括以下步骤：

步骤1：选取9位受试者的脑电信号作为训练集和测试集，分别对单个受试c3、c4两个通道中的信号进行预处理；

步骤2：对预处理后的eeg信号x(t)进行经验模式分解；得到一系列固有模态函数imfi，i为固有模态函数的阶数；

eeg信号进行经验模式分解的步骤如下：

2.1找出原信号x(t)的所有极大值点和极小值点，将其用三次样条函数拟合为原序列的上、下包络线，分别为vmax(t)和vmin(t)，可通过公式1得包络线的平均值m(t):

2.2公式2计算原信号与均值的差值c(t)：

c(t)＝x(t)-m(t)，

如果c(t)不能满足imf的定义截止条件，重复上述过程，否则提取作为固有模态函数，剩余量r(t)计算如公式3：

r(t)＝x(t)-c(t)；

2.3剩余量作为一个新的数据经过相同的筛选过程以获得下一个更低频率的固有模态函数；直到剩余函数r(t)为一个单调函数或者仅有一个极值时，分解过程停止；设原始信号x(t)被分解为n个固有模态函数和一个剩余函数量r(t)，通过公式4可以得重构信号：

步骤3：计算每次试验原始信号和各个imf分量之间相关系数的绝对值，再计算所有实验的平均值，选择具有较大相关系数绝对值平均数的固有模态函数作为特征信号；

相关系数是用以反映变量之间相关关系密切程度的统计指标；相关系数是按积差方法进行计算，同样以两变量与各自平均值的离差为基础，通过两个离差相乘来反映两变量之间的相关程度；相关系数的定义:假设有两个随机变量x和y，则它们的相关系数为公式5:

其中，cov(x,y)为随机变量x和y的协方差函数，σx、σy分别指x、y的标准差，e(x)、e(y)为两者的平均值；相关系数r的取值范围是[-1,1]，表示变量之间相关程度的高低，r的绝对值越大，说明两个变量之间的相关程度越高；r>0表示正相关，r<0表示负相关，r＝1称为完全正相关，r＝-1称为完全负相关，r＝0称为不相关；

将单次试验的c3，c4通道第1阶imf分量构成一个2*2000的矩阵xi，i＝l表示想象左手运动，i＝r表示想象右手运动，其中，2为imf个数，可看作通道数，2000为一次试验的采样点个数，即窗口长度；

步骤4：对上述imf矩阵进行功率谱密度计算，步骤如下：

4.1从无限长随机序列x(n)中截取长度n的有限长序列xn(n)；

4.2由n长序列xn(n)求(2m-1)点的自相关函数序列如公式6：

这里，m＝-(m-1)…，-1，0，1…，m-1，m≤n，rx(m)是双边序列，但是由自相关函数的偶对称性式，只要求出m＝0,……,m-1的傅里叶变换，另一半也就知道了；

4.3由相关函数的傅式变换求功率谱，即公式7：

优选地，上述技术方案中，计算选择的第1阶imf分量的psd值作为特征信号对应步骤2，psd计算包括如下步骤：

(1)从无限长随机序列x(n)中截取长度n的有限长序列xn(n)；

(2)由n长序列xn(n)求(2m-1)点的自相关函数序列，公式8：

这里，m＝-(m-1)…，-1，0，1…，m-1，m≤n，rx(m)是双边序列；

(3)由相关函数的傅式变换求功率谱，即公式9：

对输入的第1阶imf信号，即2*2000的矩阵信号，计算其功率谱密度值采样点为1000，并构成新的2*1000矩阵信号，再进行csp滤波。

优选地，上述技术方案中，所述方法的共空间模式算法包括如下步骤对应权步骤3：

(1)假设运动想象单次任务的脑电信号矩阵为x，x为n*t维，n表示脑电信号采集时的通道数，t表示脑电信号采集时每次任务各通道的采样点数。于是将x归一化处理后可得到协方差矩阵，公式10：

式中：xd表示类别样本为d的脑电信号，d∈{1,2}；表示xd的转置；tr(·)表示矩阵的迹，即矩阵的对角元素之和；

(2)对合成的协方差矩阵进行分解，公式11：

其中∑为特征值对角矩阵，u0为其对应的特征向量矩阵，将∑的特征值按降序排列，其对应的u0也重新排列；

(3)求白化矩阵.白化矩阵p定义如公式12：

(4)白化协方差矩阵r1和r2，公式13:

s1＝pr1p^t，s2＝pr2p^t；

(5)主成分分解，公式14:

其中∑1+∑2＝i，因为∑1越大，∑2越小，反之亦然，所以取∑1和∑2中最大的m个特征值对应的特征向量组成投影矩阵u，u＝(u'1u'2)，其中u'1、u'2为最大的m个特征值对应的特征向量，因此，白化后的脑电信号再经过u滤波处理后，能够得到最佳的分类特征值；

(6)求得csp空间滤波器w，公式15：

w＝u^tp，

运动想象脑电信号x经空间滤波后变为z，z＝wx，取特征向量个数为2m，并且该2m个特征向量中的前m维与一类的运动想象脑电信号的方差达到最大化，而与另外一类的运动想象脑电信号的方差则达到最小化，这样两类的样本间的距离达到最大化；后m维与之相反。其中第1个和最后1个特征向量包含了最大区别两类任务的信息，第2个和倒数第2个次之，以此类推。在实际应用中，应选取合适的向量个数，包含最佳的特征信息。

(7)脑电信号进行特征提取后，再将空间投影后的信号zp(p＝1,…,2m)进行取对数变化，从而使特征值差异更明显，公式16：

其中var表示求向量的方差。

优选地，上述技术方案中，所述方法使用支持向量机进行特征分类对应步骤4，包括如下步骤：

(1)对于非线性eeg问题，将特征集通过非线性变换转化为另一个空间中的线性问题，构造最优分类面，相应的最优决策函数为公式17：

其中n为支持向量个数，αi为lagrangue乘子，从而目标函数变为使公式18最小化：

核参数γ和误差惩罚因子c是影响svm性能的主要参数，γ的取值影响空间变换后的数据分布，惩罚因子c则决定了支持向量机的收敛速度及推广能力；因此，对γ和c的选择很大程度上影响了脑电信号的识别率；

(2)采用网格化交叉验证方法进行γ和c最优参数的选择，将训练集作为原始的数据集，在一定范围内改变核函数和惩罚因子的值，运用交叉验证方法进行分类，选择分类准确率最高的γ和c作为最佳参数；

(3)确定γ和c后，将测试集输入训练好的支持向量机进行特征分类。

与现有技术相比，本发明具有如下有益效果：

1、本发明利用三个通道emd分解后的固有模态函数进行csp滤波，在csp的基础上加入emd的频域信息，很好地解决csp缺乏频域信息的问题。

2、本发明将经验模式分解后的多阶固有模态函数看作多输入信号进行公共空间模式分解，在仅使用c3、c4、cz三个通道的情况下，获得较好的特征分类结果，解决一般csp算法需要大量输入通道问题。

3、在emd和csp脑电特征提取的基础上融合了psd算法，emd分解成多个固有模态函数并用相关系数法选择合适的imf，计算imf的psd值，对信号进行优化得到新的特征值，再进行csp滤波并进行分类。

附图说明：

图1为本发明方法流程图。

具体实施方式：

下面对本发明的具体实施方式进行详细描述，但应当理解本发明的保护范围并不受具体实施方式的限制。

除非另有其它明确表示，否则在整个说明书和权利要求书中，术语“包括”或其变换如“包含”或“包括有”等等将被理解为包括所陈述的元件或组成部分，而并未排除其它元件或其它组成部分。

下面结合说明书附图对本发明作进一步的详细说明。

如图1所示，本发明所述方法包括如下步骤：

步骤1：选取9位受试者的脑电信号作为训练集和测试集，分别对单个受试c3、c4两个通道中的信号进行预处理；

步骤2：对预处理后的eeg信号x(t)进行经验模式分解；得到一系列固有模态函数imfi(i为固有模态函数的阶数)；

eeg信号进行经验模式分解的步骤如下：

(1)找出原信号x(t)的所有极大值点和极小值点，将其用三次样条函数拟合为原序列的上、下包络线，分别为vmax(t)和vmin(t)，可得包络线的平均值m(t):

(2)计算原信号与均值的差值c(t)：

c(t)＝x(t)-m(t)(2)

如果c(t)不能满足imf的定义截止条件，重复上述过程，否则提取作为固有模态函数，剩余量r(t)计算如下：

r(t)＝x(t)-c(t)(3)

(3)剩余量作为一个新的数据经过相同的筛选过程以获得下一个更低频率的固有模态函数。直到剩余函数r(t)为一个单调函数或者仅有一个极值时，分解过程停止。假设原始信号x(t)被分解为n个固有模态函数和一个剩余函数量r(t)，可以得重构信号：

步骤3：计算每次试验原始信号和各个imf分量之间相关系数的绝对值，再计算所有实验的平均值，选择具有较大相关系数绝对值平均数的固有模态函数作为特征信号。假设有两个随机变量x和y，则它们的相关系数为:

其中，cov(x,y)为随机变量x和y的协方差函数，σx、σy分别指x、y的标准差，e(x)、e(y)为两者的平均值。相关系数r的取值范围是[-1,1]，表示变量之间相关程度的高低，r的绝对值越大，说明两个变量之间的相关程度越高。r>0表示正相关，r<0表示负相关，r＝1称为完全正相关，r＝-1称为完全负相关，r＝0称为不相关。

将单次试验的c3、c4通道第1阶imf分量构成一个2*2000的矩阵xi(i＝l表示想象左手运动，i＝r表示想象右手运动)，其中，2为imf个数，可看作通道数，2000为一次试验的采样点个数，即窗口长度；

步骤4：对上述imf矩阵进行功率谱密度计算：

(1)从无限长随机序列x(n)中截取长度n的有限长序列xn(n)。

(2)由n长序列xn(n)求(2m-1)点的自相关函数序列。

这里，m＝-(m-1)…，-1，0，1…，m-1，m≤n，rx(m)是双边序列，但是由自相关函数的偶对称性式，只要求出m＝0,……,m-1的傅里叶变换，另一半也就知道了。

(3)由相关函数的傅式变换求功率谱，即：

对输入的第1阶imf信号，即2*2000的矩阵信号，计算其功率谱密度值采样点为1000，并构成新的2*1000矩阵信号，再进行csp滤波。

共空间模式算法包括如下步骤：

(1)假设运动想象单次任务的脑电信号矩阵为x，x为n*t维，表示脑电信号采集时的通道数，t表示脑电信号采集时每次任务各通道的采样点数。于是将x归一化处理后可得到协方差矩阵：

式中：xd表示类别样本为d的脑电信号，d∈{1,2}；表示xd的转置；tr(·)表示矩阵的迹，即矩阵的对角元素之和。

(2)对合成的协方差矩阵进行分解：

其中σ为特征值对角矩阵，u0为其对应的特征向量矩阵。将σ的特征值按降序排列，其对应的u0也重新排列。

(3)求白化矩阵.白化矩阵p定义如下：

(4)白化协方差矩阵r1和r2:

s1＝pr1p^ts2＝pr2p^t(13)

(5)主成分分解:

其中∑1+∑2＝i。因为∑1越大，∑2越小，反之亦然，所以取∑1和∑2中最大的m个特征值对应的特征向量组成投影矩阵u，u＝(u'1u'2)，其中u'1、u'2为最大的m个特征值对应的特征向量。因此，白化后的脑电信号再经过u滤波处理后，能够得到最佳的分类特征值。

(6)求得csp空间滤波器w。

w＝u^tp(15)

运动想象脑电信号x经空间滤波后变为z，z＝wx，取特征向量个数为2m，并且该2m个特征向量中的前m维与一类的运动想象脑电信号的方差达到最大化，而与另外一类的运动想象脑电信号的方差则达到最小化，这样两类的样本间的距离达到最大化；后m维与之相反。其中第1个和最后1个特征向量包含了最大区别两类任务的信息，第2个和倒数第2个次之，以此类推。在实际应用中，应选取合适的向量个数，包含最佳的特征信息。

(7)脑电信号进行特征提取后，再将空间投影后的信号zp(p＝1,…,2m)进行取对数变化，从而使特征值差异更明显。

其中var表示求向量的方差。

步骤5：使用支持向量机进行特征分类，包括如下步骤：

(1)对于非线性eeg问题，将特征集通过非线性变换转化为另一个空间中的线性问题，构造最优分类面，相应的最优决策函数为：

其中n为支持向量个数，αi为lagrangue乘子，从而目标函数变为使下式最小化：

核参数γ和误差惩罚因子c是影响svm性能的主要参数，γ的取值影响空间变换后的数据分布，惩罚因子c则决定了支持向量机的收敛速度及推广能力；因此，对γ和c的选择很大程度上影响了脑电信号的识别率；

(2)采用网格化交叉验证方法进行γ和c最优参数的选择，将训练集作为原始的数据集，在一定范围内改变核函数和惩罚因子的值，运用交叉验证方法进行分类，选择分类准确率最高的γ和c作为最佳参数；

(3)确定γ和c后，将测试集输入训练好的支持向量机进行特征分类。

综上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此。在发明所披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明所揭露的技术范围之内。因此，本发明的保护范围应以权利要求书的保护范围为准。

前述对本发明的具体示例性实施方案的描述是为了说明和例证的目的。这些描述并非想将本发明限定为所公开的精确形式，并且很显然，根据上述教导，可以进行很多改变和变化。对示例性实施例进行选择和描述的目的在于解释本发明的特定原理及其实际应用，从而使得本领域的技术人员能够实现并利用本发明的各种不同的示例性实施方案以及各种不同的选择和改变。本发明的范围意在由权利要求书及其等同形式所限定。