本发明涉及emd和csp融合最优波长空间滤波脑电特征提取方法，属于智能信息处理技术领域。
背景技术：
传统的运动通道是由大脑神经和肌肉组成，神经传导冲动，肌肉配合完成相应的动作，而脑-机接口(brain-computerinterface,bci)则提供了另一种运动通道，不依赖传统的运动通道，由大脑意识直接与外部设备相连，建立运动通道，用人的大脑意识控制外部设备，不需要经过神经传导和肌肉运动，为神经损伤或肌肉损伤的病人提供了一种新的运动方式，不再需要依靠别人的照顾，自己即可完成运动。脑机接口技术领域的发展不仅可以帮助瘫痪病人使用计算机、神经假体、机械臂等电子设备,也实现了包括:运动恢复、通信、环境控制甚至娱乐等其他功能。
脑机接口技术主要包含信号采集、预处理、特征提取、特征分类和接口设备控制等五个步骤。其中，由特征提取得到的特征信号可以标识不同想象运动脑电信号的判别信息，对后续的分类识别有很大的影响，因此特征提取在bci研究界受到广泛的关注。
有效的特征提取方法是改善识别精度的关键，目前，时间-频率方法被广泛应用于脑电信号的研究。传统的时间-频率方法包括：短时傅里叶变换(short-timefouriertransform,stft)，小波变换(wavelettransform,wt)等，但是这些方法的本质都是基于傅里叶变换，根据海森堡不确定性原理，该方法不可能同时得到时间-频率的良好分辨率。近年，希尔伯特黄变换(hilbert-huangtransform,hht)作为另一种时间-频率分析法已经变得越来越流行，同时它也十分适合分析非线性和非平稳信号。原信号经过经验模式分解(empiricalmodedecomposition,emd)被分解为一系列固有模态函数(intrinsicmodefunctions,imfs)，随后对每个固有模态函数进行希尔伯特黄变换，求其相应的能量谱和边际谱，作为特征进行分类。hht不涉及海森堡不确定性原理可以获得时域和频域的高分辨率。目前被广泛应用于许多信号处理领域，如雷达探测、地震信号和生物医学信号等。
再者，由于eeg信号低空间分辨率，eeg信号构成的bci系统需要进行有效的空间滤波，从而确保从受试的相关脑域中提取特征信息。在这一方面，常用的算法有：共空间模式(commonspatialpattern,csp)、独立主成分分析(independentcomponentanalysis,ica)和共域空间谱模式(commonspatialspectralpattern，cssp)、滤波器csp(filterbankcommonspatialpattern,fbcsp)、判别滤波csp(discriminantfilteringcommonspatialpattern,dfbcsp)等多种csp改进算法。
然而，传统的csp需要大量的输入通道，同时缺乏频率信息。csp算法主要是通过对脑电源的空间信息的判定来实现不同分类间的差异极大化，算法本质是利用代数上矩阵同时对角化的理论，寻找一组空间滤波器，脑电信号通过这组滤波器的投影获得较明显的特征向量。由于脑电信号是极复杂的非线性，不平稳信号，采集到的信号中，不止有脑电信号，还有其他的杂质信号，即使在脑电信号的频带范围内依旧存在其他的混叠信号，这些杂质与混叠信号将会影响特征提取的有效性，而单纯的csp滤波，只是通过空间信息来判定脑电信号的有效性，存在一定的不足，可能不会将杂质与混叠信号从频段中去除干净，如果许多与运动想象无关的频率信号混在其中，严重影响了特征向量的有效性。而本发明能够很好的解决上面的问题。
技术实现要素：
本发明针对上述现有技术的不足，提出了一种emd和csp融合最优波长空间滤波脑电特征提取方法，大大提供了模式识别的准确率。该方法是基于如下内容：1)基于经验模式分解的eeg信号处理。2)根据频谱分析筛选固有模态函数，组成新的信号矩阵。3)对信号矩阵进行最优波长计算，再进行公共空间模式分解，解决csp多输入、缺频域信息的问题。4)支持向量机分类。
本发明解决其技术问题所采取的技术方案是：本发明首先对预处理后信号进行经验模式分解，根据各阶固有模态函数的频谱分析筛选出集中在想象运动频段的前2阶固有模态函数，对其进行最优波长计算，再进一步进行csp空间滤波，通过支持向量机进行特征选择得到最终分类结果。
本发明对通道信号的提取进行了改进包括：
对于左手想象运动，提取c3和c4两个通道的数据作为原始信号；对于右手想象运动，只提取c3通道的数据作为原始信号。将右手想象运动的c3通道分别与左手想象运动的c3，c4通道进行对比，提取特征向量，进行分类识别。
本发明对传统的最优波长计算进行了改进，包括：
对信号进行最优波长计算，如下
2：n1：n-1
δxi＝xi-xi
其中n为信号的采样点数，xi:j代表信号矩阵从第i列到第j列，同时保留第一列的原始数据；传统的是以矩阵行进行变换。
方法流程：
步骤1：选取9位受试的脑电信号作为训练集和测试集，分别对单个受试c3、c4两个通道中的信号进行预处理；
步骤2：对预处理后的eeg信号x(t)进行经验模式分解；得到一系列固有模态函数imfi(i为固有模态函数的阶数)并绘制所有固有模态函数能谱图；
eeg信号进行经验模式分解的具体步骤如下：
(1)判断每个x(t)的局部极值，用三次样条曲线进行曲线拟合，局部极大值形成上包络emax(t),局部极小值形成下包络emin(t)。
(2)求emax(t)和emin(t)的均值：
(3)计算输入信号x(t)和m(t)的差值:
c(t)＝x(t)-m(t)(2)
如果c(t)不能满足imf的定义截至条件，重复上述过程(1)-(3)，否则，提取c(t)作为固有模态函数，剩余量r(t)计算如下：
r(t)＝x(t)-c(t)(3)
(4)剩余量作为一个新的数据经过相同的筛选过程以获得下一个更低频率的固有模态函数。直到剩余函数r(t)为一个单调函数或者仅有一个极致时，分解过程停止。假设原始信号x(t)被分解为n个固有模态函数和一个剩余函数量r(t)，可以得重构信号：
步骤3：将单次试验的c3，c4通道前2阶imf分量合并，构成一个4*2000的矩阵xi(i＝l表示想象左手运动，i＝r表示想象右手运动)，其中，4为imf个数，可看作通道数，2000为一次试验的采样点个数，即窗口长度。
步骤4：对上述imf矩阵进行最优波长计算：
2：n1：n-1
δxi＝xi-xi(5)
其中n为信号的采样点数，xi:j代表信号矩阵从第i列到第j列；对变换后的矩阵进行公共空间模式分解；
公共空间模式算法具体过程如下：
对两类想象运动a和b分别进行ta,tb次实验，ta,tb为正整数；
(1)计算混合空间的协方差；
首先，计算两类运动想象信号每次实验的协方差，公式如下：
其中，trace(xxt)为矩阵xxt的迹，即矩阵xxt的对角线元素之和；
然后，分别计算两类运动想象的平均协方差：
其中，ca,i、cb,i分别表示运动想象a和b的第i次实验的协方差；
进而求得混合空间的协方差：
cm＝ca+cb(8)
(2)对混合空间协方差进行特征值分解，公式如下：
其中，um为特征向量矩阵，λm为特征值对角矩阵；
(3)进行白化处理；
对λm进行降序排序得到λmd，并对um做同样的行列变换得到umd；令
sa＝pcapt
sb＝pcbpt(10)
利用sa、sb具有相同的特征向量的特点，经特征值分解后可得：
sa＝bλabt
sb＝bλbbt(11)
其中，b为sa与sb的共同特征向量，λa、λb分别为sa和sb的特征对角矩阵，且λa+λb＝i，i为单位矩阵；
因此求得空间滤波器矩阵为：
w＝btp(12)
对x进行w滤波得：
z0＝wx(13)
(4)求特征向量f；
提取z0的前m行和后m行，构成z＝[z1,z2,…,z2m]t,然后进行特征提取，计算公式如下：
其中，var(·)表示方差，i＝1,2,…,2m，则特征向量f＝[f1,f2,…,f2m]t；步骤5：使用支持向量机进行特征分类，其具体步骤如下：
(1)对于非线性eeg问题，将特征集通过非线性变换转化为另一个空间中的线性问题，构造最优分类面。相应的最优决策函数为：
其中n为支持向量个数，αi为lagrangue乘子。从而目标函数变为使下式最小化：
核参数γ和误差惩罚因子c是影响svm性能的主要参数。γ的取值影响空间变换后的数据分布，惩罚因子c则决定了支持向量机的收敛速度及推广能力；因此，对γ和c的选择很大程度上影响了脑电信号的识别率。
(2)采用网格化交叉验证方法进行γ和c最优参数的选择，将训练集作为原始的数据集，在一定范围内改变核函数和惩罚因子的值，运用交叉验证方法进行分类，选择分类准确率最高的γ和c作为最佳参数。
(3)确定γ和c后，将测试集输入训练好的支持向量机进行特征分类。
有益效果：
1、本发明利用三个通道emd分解后的固有模态函数进行csp滤波，在csp的基础上加入emd的频域信息，很好地解决csp缺乏频域信息的问题。
2、本发明将经验模式分解后的多阶固有模态函数看作多输入信号进行公共空间模式分解，在仅使用c3，c4，cz三个通道的情况下，获得较好的特征分类结果，解决一般csp算法需要大量输入通道问题。
3、针对最优波长算法做出改进：由于矩阵行数较少，本发明将行计算改为列计算，同时保留第一列的原始数据，为后续算法做准备。
附图说明
图1为本发明方法流程图。
具体实施方式
下面结合说明书附图对本发明作进一步的详细说明。
如图1所示，本发明所述方法包括如下步骤：
步骤1：采集各位受试的脑电信号。选取9位受试的脑电信号作为训练集和测试集，分别对单个受试c3、c4两个通道中的信号进行预处理；
步骤2：对预处理后的eeg信号x(t)进行经验模式分解；得到一系列固有模态函数imfi(i为固有模态函数的阶数)并绘制所有固有模态函数能谱图；
eeg信号进行经验模式分解的具体步骤如下：
(1)判断每个x(t)的局部极值，用三次样条曲线进行曲线拟合，局部极大值形成上包络emax(t),局部极小值形成下包络emin(t)。
(2)求emax(t)和emin(t)的均值：
(3)计算输入信号x(t)和m(t)的差值:
c(t)＝x(t)-m(t)(2)
如果c(t)不能满足imf的定义截至条件，重复上述过程(1)-(3)，否则，提取c(t)作为固有模态函数，剩余量r(t)计算如下：
r(t)＝x(t)-c(t)(3)
(4)剩余量作为一个新的数据经过相同的筛选过程以获得下一个更低频率的固有模态函数。直到剩余函数r(t)为一个单调函数或者仅有一个极致时，分解过程停止。假设原始信号x(t)被分解为n个固有模态函数和一个剩余函数量r(t)，可以得重构信号：
步骤3：将次试验的c3，c4通道前2阶imf分量合并，构成一个4*2000的矩阵xi(i＝l表示想象左手运动，i＝r表示想象右手运动)，其中，4为imf个数，可看作通道数，2000为一次试验的采样点个数，即窗口长度。
步骤4：对上述imf矩阵进行最优波长计算：
2：n1：n-1
δxi＝xi-xi(5)
其中n为信号的采样点数，xi:j代表信号矩阵从第i列到第j列；对变换后的矩阵进行公共空间模式分解；
公共空间模式算法具体过程如下：
对两类想象运动a和b分别进行ta,tb次实验，ta,tb为正整数；
(1)计算混合空间的协方差；
首先，计算两类运动想象信号每次实验的协方差，公式如下：
其中，trace(xxt)为矩阵xxt的迹，即矩阵xxt的对角线元素之和；
然后，分别计算两类运动想象的平均协方差：
其中，ca,i、cb,i分别表示运动想象a和b的第i次实验的协方差；
进而求得混合空间的协方差：
cm＝ca+cb(8)
(2)对混合空间协方差进行特征值分解，公式如下：
其中，um为特征向量矩阵，λm为特征值对角矩阵；
(3)进行白化处理；
对λm进行降序排序得到λmd，并对um做同样的行列变换得到umd；令
sa＝pcapt
sb＝pcbpt(10)
利用sa、sb具有相同的特征向量的特点，经特征值分解后可得：
sa＝bλabt
sb＝bλbbt(11)
其中，b为sa与sb的共同特征向量，λa、λb分别为sa和sb的特征对角矩阵，且λa+λb＝i，i为单位矩阵；
因此求得空间滤波器矩阵为：
w＝btp(12)
对x进行w滤波得：
z0＝wx(13)
(4)求特征向量f；
提取z0的前m行和后m行，构成z＝[z1,z2,…,z2m]t,然后进行特征提取，计算公式如下：
其中，var(·)表示方差，i＝1,2,…,2m，则特征向量f＝[f1,f2,…,f2m]t；步骤5：使用支持向量机进行特征分类，其具体步骤如下：
(1)对于非线性eeg问题，将特征集通过非线性变换转化为另一个空间中的线性问题，构造最优分类面。相应的最优决策函数为：
其中n为支持向量个数，αi为lagrangue乘子。从而目标函数变为使下式最小化：
核参数γ和误差惩罚因子c是影响svm性能的主要参数。γ的取值影响空间变换后的数据分布，惩罚因子c则决定了支持向量机的收敛速度及推广能力；因此，对γ和c的选择很大程度上影响了脑电信号的识别率。
(2)采用网格化交叉验证方法进行γ和c最优参数的选择，将训练集作为原始的数据集，在一定范围内改变核函数和惩罚因子的值，运用交叉验证方法进行分类，选择分类准确率最高的γ和c作为最佳参数。
(3)确定γ和c后，将测试集输入训练好的支持向量机进行特征分类。
综上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此。在发明所披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明所揭露的技术范围之内。因此，本发明的保护范围应以权利要求书的保护范围为准。