CN113499087A - 一种基于itd模态参数识别的脑电信号特征提取方法 - Google Patents

Abstract

本发明公开了一种基于ITD模态参数识别的脑电信号特征提取方法，其根据脑电信号中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号各自的频率范围，设计各自对应的有限脉冲响应带通滤波器，再利用各自对应的有限脉冲响应带通滤波器对脑电信号进行滤波处理，提取得到δ波节律信号、θ波节律信号、α波节律信号、β波节律信号；对四个节律信号分别进行ITD模态参数识别，提取得到各自的四个ITD模态参数，进而得到四个节律信号各自的特征向量；将四个节律信号的特征向量构成脑电信号的特征向量；优点是其能够从运动想象的脑电信号中提取出有效的特征，从而能够极大地提高后续特征分类的正确率。

Description

一种基于ITD模态参数识别的脑电信号特征提取方法

技术领域

本发明涉及一种信号特征提取技术，尤其是涉及一种基于ITD(Ibrahim TimeDomain)模态参数识别的脑电信号特征提取方法。

背景技术

脑-机接口(Brain-Computer Interface，BCI)是一种能够实现人类大脑与外部环境直接进行交流的通讯控制系统，这种通讯控制系统能够使人类大脑脱离对脑部外周神经和肌肉组成的大脑正常输出通路的依赖。近些年来，由于计算机科学、生物医学、信号处理以及人工智能等众多技术的广泛发展，因此人们也逐渐的开始关注作为其交叉学科的脑-机接口技术。

脑-机接口技术包括五个主要部分：脑电信号采集、预处理、特征提取、特征分类与识别以及控制设备。其中，如何有效地提取脑电信号的特征以及提高脑电信号分类的正确率是目前技术研究的重点。至今，人们利用众多信号处理技术的方法如CSP(CommonSpatial Pattern，共空间模式)算法、稀疏自编码方法等来提取脑电信号的特征，使脑电信号的定量分析取得了重大进展，然而对于大脑中是否存在动力学行为一直是热门话题。以动力学的观点来看，人类大脑的神经元网络系统是由许多神经元连接形成的极度错综复杂的一种超高维的非线性系统。鉴于生命现象的多重耦合性质，连同生物材料的活性、大阻尼特性等多种非线性特性，这使得生命系统具有强烈的非线性，因而可以从这种非线性特性深入研究，这将对信号处理技术和生命现象研究都具有深远意义。

发明内容

本发明所要解决的技术问题是基于脑电信号的非线性振动理论，提供一种基于ITD模态参数识别的脑电信号特征提取方法，其能够从运动想象的脑电信号中提取出有效的特征，从而能够极大地提高后续特征分类的正确率。

本发明解决上述技术问题所采用的技术方案为：一种基于ITD模态参数识别的脑电信号特征提取方法，其特征在于包括以下步骤：

步骤1：按给定的采样频率在一个通道一次随机采集脑电信号；然后将脑电信号以向量形式进行表示，记为Y，Y＝[y₁…y_n…y_N]；其中，给定的采样频率为F赫兹，Y的维数为1×N，Y中的每个元素为一个样本值，N表示Y中的样本值的总个数，N≥5，N＝F×t，t表示采集一个通道的一条脑电信号的采集时间，t的单位为秒，n为正整数，n的初始值为1，1≤n≤N，y₁表示Y中的第1个样本值，y_n表示Y中的第n个样本值，y_N表示Y中的第N个样本值，符号“[]”为向量或矩阵的表示符号；

步骤2：根据脑电信号中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号各自的频率范围，设计各自对应的有限脉冲响应带通滤波器；然后利用脑电信号中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号各自对应的有限脉冲响应带通滤波器分别对Y进行滤波处理，对应提取得到Y中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号，以向量形式对应表示为Y_δ、Y_θ、Y_α、Y_β，Y_δ＝[y_δ,1…y_δ,n…y_δ,N]，Y_θ＝[y_θ,1…y_θ,n…y_θ,N]，Y_α＝[y_α,1…y_α,n…y_α,N]，Y_β＝[y_β,1…y_β,n…y_β,N]；其中，脑电信号中的δ波节律信号的频率范围为0.78～3.90赫兹，脑电信号中的θ波节律信号的频率范围为3.91～7.80赫兹，脑电信号中的α波节律信号的频率范围为7.81～13.28赫兹，脑电信号中的β波节律信号的频率范围为13.29～30.47赫兹，Y_δ、Y_θ、Y_α、Y_β的维数均为1×N，Y_δ中的每个元素为一个样本值，y_δ,1表示Y_δ中的第1个样本值，y_δ,n表示Y_δ中的第n个样本值，y_δ,N表示Y_δ中的第N个样本值，Y_θ中的每个元素为一个样本值，y_θ,1表示Y_θ中的第1个样本值，y_θ,n表示Y_θ中的第n个样本值，y_θ,N表示Y_θ中的第N个样本值，Y_α中的每个元素为一个样本值，y_α,1表示Y_α中的第1个样本值，y_α,n表示Y_α中的第n个样本值，y_α,N表示Y_α中的第N个样本值，Y_β中的每个元素为一个样本值，y_β,1表示Y_β中的第1个样本值，y_β,n表示Y_β中的第n个样本值，y_β,N表示Y_β中的第N个样本值；

步骤3：对Y_δ进行ITD模态参数识别，提取得到Y_δ的四个ITD模态参数，将Y_δ的四个ITD模态参数按序排列构成的向量作为Y_δ的特征向量，记为fea_δ；以相同的方法对Y_θ进行ITD模态参数识别，提取得到Y_θ的四个ITD模态参数，将Y_θ的四个ITD模态参数按序排列构成的向量作为Y_θ的特征向量，记为fea_θ；以相同的方法对Y_α进行ITD模态参数识别，提取得到Y_α的四个ITD模态参数，将Y_α的四个ITD模态参数按序排列构成的向量作为Y_α的特征向量，记为fea_α；以相同的方法对Y_β进行ITD模态参数识别，提取得到Y_β的四个ITD模态参数，将Y_β的四个ITD模态参数按序排列构成的向量作为Y_β的特征向量，记为fea_β；其中，fea_δ、fea_θ、fea_α、fea_β的维数均为1×4；

步骤4：将fea_δ、fea_θ、fea_α、fea_β按序排列构成的向量作为Y的特征向量，记为fea，fea＝[fea_δfea_θfea_αfea_β]；其中，fea的维数为1×16。

所述的步骤2的具体过程为：

根据脑电信号中的δ波节律信号的频率范围，设计第一有限脉冲响应带通滤波器，记为H1，H1的采样频率为F赫兹，H1的阶数为10，H1的阻带上边缘频率为0.7赫兹、阻带下边缘频率为4赫兹、通带上边缘频率为0.8赫兹、通带下边缘频率为3.5赫兹；然后利用H1对Y进行滤波处理，提取得到Y中的δ波节律信号；

同样，根据脑电信号中的θ波节律信号的频率范围，设计第二有限脉冲响应带通滤波器，记为H2，H2的采样频率为F赫兹，H2的阶数为10，H2的阻带上边缘频率为3.9赫兹、阻带下边缘频率为8赫兹、通带上边缘频率为4赫兹、通带下边缘频率为7.5赫兹；然后利用H2对Y进行滤波处理，提取得到Y中的θ波节律信号；

根据脑电信号中的α波节律信号的频率范围，设计第三有限脉冲响应带通滤波器，记为H3，H3的采样频率为F赫兹，H3的阶数为10，H3的阻带上边缘频率为7.8赫兹、阻带下边缘频率为13赫兹、通带上边缘频率为7.9赫兹、通带下边缘频率为12赫兹；然后利用H3对Y进行滤波处理，提取得到Y中的α波节律信号；

根据脑电信号中的β波节律信号的频率范围，设计第四有限脉冲响应带通滤波器，记为H4，H4的采样频率为F赫兹，H4的阶数为10，H4的阻带上边缘频率为13赫兹、阻带下边缘频率为31赫兹、通带上边缘频率为14赫兹、通带下边缘频率为30赫兹；然后利用H4对Y进行滤波处理，提取得到Y中的β波节律信号。

所述的步骤3中，fea_δ的获取过程为：

步骤3_1：将Y_δ＝[y_δ,1…y_δ,n…y_δ,N]离散信号看作自由响应振动信号，对Y_δ＝[y_δ,1…y_δ,n…y_δ,N]依次进行正常采样、d延时釆样、2d延时釆样，依次得到正常采样自由响应数据向量、d延时釆样自由响应数据向量、2d延时釆样自由响应数据向量，对应记为X_nor、X_{del_1}、X_{del_2}，X_nor＝[y_δ,1…y_δ,n…y_δ,m]，X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]，X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]；然后将X_nor＝[y_δ,1…y_δ,n…y_δ,m]重新表示为X_nor＝[x_nor,1…x_nor,n…x_nor,m]，将X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]重新表示为X_{del_1}＝[x_{del_1,1}…x_{del_1,n}…x_{del_1,m}]，将X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]重新表示为X_{del_2}＝[x_{del_2,1}…x_{del_2,n}…x_{del_2,m}]；其中，d为正整数，d表示采样间隔，

m为正整数，m表示X_nor、X_{del_1}、X_{del_2}的长度，1≤m≤N-2d，X_nor、X_{del_1}、X_{del_2}的维数均为1×m，y_δ,m、y_δ,d+1、y_δ,d+n、y_δ,d+m、y_δ,2d+1、y_δ,2d+n、y_δ,2d+m对应表示Y_δ中的第m个样本值、第d+1个样本值、第d+n个样本值、第d+m个样本值、第2d+1个样本值、第2d+n个样本值、第2d+m个样本值，x_nor,1、x_nor,n、x_nor,m对应表示X_nor中的第1个元素、第n个元素、第m个元素，x_{del_1,1}、x_{del_1,n}、x_{del_1,m}对应表示X_{del_1}中的第1个元素、第n个元素、第m个元素，x_{del_2,1}、x_{del_2,n}、x_{del_2,m}对应表示X_{del_2}中的第1个元素、第n个元素、第m个元素；

步骤3_2：根据X_nor和X_{del_1}构建第一增广矩阵，记为X_zg1，

并根据X_{del_1}和X_{del_2}构建第二增广矩阵，记为X_zg2，

然后令A＝X_zg2(X_zg1)^-1，并利用QR分解法对A＝X_zg2(X_zg1)^-1进行处理，得到A的特征值矩阵，记为Δ；其中，X_zg1和X_zg2的维数均为2×m，A表示系统特征矩阵，(X_zg1)^-1为X_zg1的伪逆，A的维数为2×2，Δ的维数为2×2；

步骤3_3：根据Δ计算Y_δ的固有频率和衰减系数，对应记为f_δ和λ_δ，

其中，

表示

的虚数部分，

表示

的实数部分，Δ_1,1表示Δ的第1行第1列元素，t'表示在Y_δ＝[y_δ,1…y_δ,n…y_δ,N]中采样d个样本值的时间间隔，

步骤3_4：应用自由响应理论将X_nor描述为：X_nor＝ψ_norE_nor，

其中，ψ_nor表示X_nor的模态矢量，ψ_nor的维数为1×2，E_nor为引入的中间变量，E_nor的维数为2×m，e表示自然基数，μ＝-λ_δ+j2πf_δ，μ^*为μ的共轭，μ^*＝-λ_δ-j2πf_δ，j表示虚数单位，

然后利用最小二乘法求解X_nor＝ψ_norE_nor，得到ψ_nor，ψ_nor＝X_norE_nor ^T(E_norE_nor ^T)^-1，其中，上标“T”表示向量或矩阵的转置，(E_norE_nor ^T)^-1为E_norE_nor ^T的逆；再计算ψ_nor的幅度和相位，对应记为ξ_nor和

ξ_nor＝abs(ψ_nor(1))+abs(ψ_nor(2))，

其中，abs()为取绝对值函数，ψ_nor(1)表示ψ_nor的第1个元素，ψ_nor(2)表示ψ_nor的第2个元素，

表示

的实数部分；最后将ξ_nor作为Y_δ的幅度，重新记为ξ_δ，将

作为Y_δ的相位，重新记为

步骤3_5：f_δ、λ_δ、ξ_δ、

为Y_δ的四个ITD模态参数，将f_δ、λ_δ、ξ_δ、

按序排列构成的向量作为Y_δ的特征向量fea_δ，

所述的步骤3中，fea_θ、fea_α、fea_β的获取为按照步骤3_1至步骤3_5的过程以相同的方式获取，

其中，f_θ、λ_θ、ξ_θ、

对应为Y_θ的固有频率、衰减系数、幅度、相位，f_α、λ_α、ξ_α、

对应为Y_α的固有频率、衰减系数、幅度、相位，f_β、λ_β、ξ_β、

对应为Y_β的固有频率、衰减系数、幅度、相位。

按给定的采样频率在同一个通道多次随机采集脑电信号，并按照步骤1至步骤4的过程以相同的方式获取在同一个通道每次采集的脑电信号的特征向量，将在同一个通道多次随机采集的脑电信号的特征向量构成一个特征矩阵，该特征矩阵中的每行向量代表一次随机采集的脑电信号的特征向量；然后对该特征矩阵进行分类，得到每次随机采集的脑电信号对应的运动想象类型。

针对多个不同通道，按给定的采样频率同时在多个通道采集脑电信号，并按照步骤1至步骤4的过程以相同的方式获取在每个通道采集的脑电信号的特征向量，将同时在多个通道采集的脑电信号的特征向量按序排列构成行向量；在进行多次采集后，得到多个行向量，将多个行向量按序排列构成一个特征矩阵，该特征矩阵中的每行向量代表同时在多个通道单次采集的脑电信号的特征向量按序排列构成的行向量；然后对该特征矩阵进行分类，得到每次采集的脑电信号对应的运动想象类型。

与现有技术相比，本发明的优点在于：

1)本发明方法从非线性振动理论出发，对脑电信号动力学模型进行分析，从振动利用工程学的角度对脑电信号加以研究，具有重要的科学理论价值，对于脑电信号特征提取具有重要意义。

2)利用本发明方法提取的脑电信号的特征用于分类，在现有的运动想象脑电信号数据库上进行仿真对比表明，本发明方法获取的振动模态特征具有更高的分类准确性，更好的稳定性。

附图说明

图1为本发明方法的总体实现框图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明提出的一种基于ITD模态参数识别的脑电信号特征提取方法，其总体实现框图如图1所示，其包括以下步骤：

步骤1：按给定的采样频率在一个通道一次随机采集脑电信号；然后将脑电信号以向量形式进行表示，记为Y，Y＝[y₁…y_n…y_N]；其中，给定的采样频率为F赫兹，给定的采样频率即采集系统的采样频率参数，Y的维数为1×N，Y中的每个元素为一个样本值，N表示Y中的样本值的总个数，N≥5，N＝F×t，t表示采集一个通道的一条脑电信号的采集时间，t的单位为秒，n为正整数，n的初始值为1，1≤n≤N，y₁表示Y中的第1个样本值，y_n表示Y中的第n个样本值，y_N表示Y中的第N个样本值，符号“[]”为向量或矩阵的表示符号。

步骤2：根据脑电信号中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号各自的频率范围(即频带)，设计各自对应的有限脉冲响应(FIR)带通滤波器；然后利用脑电信号中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号各自对应的有限脉冲响应带通滤波器分别对Y进行滤波处理，对应提取得到Y中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号，以向量形式对应表示为Y_δ、Y_θ、Y_α、Y_β，Y_δ＝[y_δ,1…y_δ,n…y_δ,N]，Y_θ＝[y_θ,1…y_θ,n…y_θ,N]，Y_α＝[y_α,1…y_α,n…y_α,N]，Y_β＝[y_β,1…y_β,n…y_β,N]；其中，脑电信号中的δ波节律信号的频率范围为0.78～3.90赫兹，脑电信号中的θ波节律信号的频率范围为3.91～7.80赫兹，脑电信号中的α波节律信号的频率范围为7.81～13.28赫兹，脑电信号中的β波节律信号的频率范围为13.29～30.47赫兹，Y_δ、Y_θ、Y_α、Y_β的维数均为1×N，Y_δ中的每个元素为一个样本值，y_δ,1表示Y_δ中的第1个样本值，y_δ,n表示Y_δ中的第n个样本值，y_δ,N表示Y_δ中的第N个样本值，Y_θ中的每个元素为一个样本值，y_θ,1表示Y_θ中的第1个样本值，y_θ,n表示Y_θ中的第n个样本值，y_θ,N表示Y_θ中的第N个样本值，Y_α中的每个元素为一个样本值，y_α,1表示Y_α中的第1个样本值，y_α,n表示Y_α中的第n个样本值，y_α,N表示Y_α中的第N个样本值，Y_β中的每个元素为一个样本值，y_β,1表示Y_β中的第1个样本值，y_β,n表示Y_β中的第n个样本值，y_β,N表示Y_β中的第N个样本值。

在本实施例中，步骤2的具体过程为：

根据脑电信号中的δ波节律信号的频率范围，设计第一有限脉冲响应带通滤波器，记为H1，H1的采样频率为F赫兹，H1的阶数为10，H1的阻带上边缘频率为0.7赫兹、阻带下边缘频率为4赫兹、通带上边缘频率为0.8赫兹、通带下边缘频率为3.5赫兹；然后利用H1对Y进行滤波处理，提取得到Y中的δ波节律信号。

同样，根据脑电信号中的θ波节律信号的频率范围，设计第二有限脉冲响应带通滤波器，记为H2，H2的采样频率为F赫兹，H2的阶数为10，H2的阻带上边缘频率为3.9赫兹、阻带下边缘频率为8赫兹、通带上边缘频率为4赫兹、通带下边缘频率为7.5赫兹；然后利用H2对Y进行滤波处理，提取得到Y中的θ波节律信号。

根据脑电信号中的α波节律信号的频率范围，设计第三有限脉冲响应带通滤波器，记为H3，H3的采样频率为F赫兹，H3的阶数为10，H3的阻带上边缘频率为7.8赫兹、阻带下边缘频率为13赫兹、通带上边缘频率为7.9赫兹、通带下边缘频率为12赫兹；然后利用H3对Y进行滤波处理，提取得到Y中的α波节律信号。

根据脑电信号中的β波节律信号的频率范围，设计第四有限脉冲响应带通滤波器，记为H4，H4的采样频率为F赫兹，H4的阶数为10，H4的阻带上边缘频率为13赫兹、阻带下边缘频率为31赫兹、通带上边缘频率为14赫兹、通带下边缘频率为30赫兹；然后利用H4对Y进行滤波处理，提取得到Y中的β波节律信号。

步骤3：对Y_δ进行ITD模态参数识别，提取得到Y_δ的四个ITD模态参数，将Y_δ的四个ITD模态参数按序排列构成的向量作为Y_δ的特征向量，记为fea_δ；以相同的方法对Y_θ进行ITD模态参数识别，提取得到Y_θ的四个ITD模态参数，将Y_θ的四个ITD模态参数按序排列构成的向量作为Y_θ的特征向量，记为fea_θ；以相同的方法对Y_α进行ITD模态参数识别，提取得到Y_α的四个ITD模态参数，将Y_α的四个ITD模态参数按序排列构成的向量作为Y_α的特征向量，记为fea_α；以相同的方法对Y_β进行ITD模态参数识别，提取得到Y_β的四个ITD模态参数，将Y_β的四个ITD模态参数按序排列构成的向量作为Y_β的特征向量，记为fea_β；其中，fea_δ、fea_θ、fea_α、fea_β的维数均为1×4。

在本实施例中，步骤3中，fea_δ的获取过程为：

步骤3_1：从运动想象的脑电信号的诱发过程来分析，当出现运动想象诱导信号刺激时，大脑开始响应这时大脑神经元突触传递信号总数急剧增加，对应的动力学方程表现为强非线性，这段时间采集到的脑电信号为瞬态非线性信号，具有混沌特性，非常难以分类处理；然后在适应了诱发刺激信号后，大脑神经元进入抑制阶段突触传递信号总数下降，对应的动力学方程中，不断减小，这段时间采集到的脑电信号为自由响应振动信号，动力学方程退化为弱非线性的振动方程。因此将Y_δ＝[y_δ,1…y_δ,n…y_δ,N]离散信号看作自由响应振动信号，对Y_δ＝[y_δ,1…y_δ,n…y_δ,N]依次进行正常采样、d延时釆样、2d延时釆样，依次得到正常采样自由响应数据向量、d延时釆样自由响应数据向量、2d延时釆样自由响应数据向量，对应记为X_nor、X_{del_1}、X_{del_2}，X_nor＝[y_δ,1…y_δ,n…y_δ,m]，X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]，X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]；然后将X_nor＝[y_δ,1…y_δ,n…y_δ,m]重新表示为X_nor＝[x_nor,1…x_nor,n…x_nor,m]，将X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]重新表示为X_{del_1}＝[x_{del_1,1}…x_{del_1,n}…x_{del_1,m}]，将X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]重新表示为X_{del_2}＝[x_{del_2,1}…x_{del_2,n}…x_{del_2,m}]；其中，d为正整数，d表示采样间隔，

m为正整数，m表示X_nor、X_{del_1}、X_{del_2}的长度，1≤m≤N-2d，X_nor、X_{del_1}、X_{del_2}的维数均为1×m，y_δ,m、y_δ,d+1、y_δ,d+n、y_δ,d+m、y_δ,2d+1、y_δ,2d+n、y_δ,2d+m对应表示Y_δ中的第m个样本值、第d+1个样本值、第d+n个样本值、第d+m个样本值、第2d+1个样本值、第2d+n个样本值、第2d+m个样本值，x_nor,1、x_nor,n、x_nor,m对应表示X_nor中的第1个元素、第n个元素、第m个元素，x_{del_1,1}、x_{del_1,n}、x_{del_1,m}对应表示X_{del_1}中的第1个元素、第n个元素、第m个元素，x_{del_2,1}、x_{del_2,n}、x_{del_2,m}对应表示X_{del_2}中的第1个元素、第n个元素、第m个元素。

步骤3_2：根据X_nor和X_{del_1}构建第一增广矩阵，记为X_zg1，

并根据X_{del_1}和X_{del_2}构建第二增广矩阵，记为X_zg2，

然后令A＝X_zg2(X_zg1)^-1，并利用QR(正交三角)分解法对A＝X_zg2(X_zg1)^-1进行处理，得到A的特征值矩阵，记为Δ；其中，X_zg1和X_zg2的维数均为2×m，A表示系统特征矩阵，(X_zg1)^-1为X_zg1的伪逆，A的维数为2×2，Δ的维数为2×2。

步骤3_3：根据Δ计算Y_δ的固有频率和衰减系数，对应记为f_δ和λ_δ，

其中，

表示

的虚数部分，

表示

的实数部分，Δ_1,1表示Δ的第1行第1列元素，t'表示在Y_δ＝[y_δ,1…y_δ,n…y_δ,N]中采样d个样本值的时间间隔，

步骤3_4：应用自由响应理论将X_nor描述为：X_nor＝ψ_norE_nor，

其中，ψ_nor表示X_nor的模态矢量，ψ_nor的维数为1×2，E_nor为引入的中间变量，E_nor的维数为2×m，e表示自然基数，e＝2.17…，μ＝-λ_δ+j2πf_δ，μ^*为μ的共轭，μ^*＝-λ_δ-j2πf_δ，j表示虚数单位，

然后利用最小二乘法求解X_nor＝ψ_norE_nor，得到ψ_nor，ψ_nor＝X_norE_nor ^T(E_norE_nor ^T)^-1，其中，上标“T”表示向量或矩阵的转置，(E_norE_nor ^T)^-1为E_norE_nor ^T的逆；再计算ψ_nor的幅度和相位，对应记为ξ_nor和

ξ_nor＝abs(ψ_nor(1))+abs(ψ_nor(2))，

其中，abs()为取绝对值函数，ψ_nor(1)表示ψ_nor的第1个元素，ψ_nor(2)表示ψ_nor的第2个元素，

表示

的实数部分；最后将ξ_nor作为Y_δ的幅度，重新记为ξ_δ，将

作为Y_δ的相位，重新记为

步骤3_5：f_δ、λ_δ、ξ_δ、

为Y_δ的四个ITD模态参数，将f_δ、λ_δ、ξ_δ、

按序排列构成的向量作为Y_δ的特征向量fea_δ，

在本实施例中，步骤3中，fea_θ、fea_α、fea_β的获取为按照步骤3_1至步骤3_5的过程以相同的方式获取，

其中，f_θ、λ_θ、ξ_θ、

对应为Y_θ的固有频率、衰减系数、幅度、相位，f_α、λ_α、ξ_α、

对应为Y_α的固有频率、衰减系数、幅度、相位，f_β、λ_β、ξ_β、

对应为Y_β的固有频率、衰减系数、幅度、相位。

步骤4：将fea_δ、fea_θ、fea_α、fea_β按序排列构成的向量作为Y的特征向量，记为fea，fea＝[fea_δ fea_θ fea_α fea_β]；其中，fea的维数为1×16。

在上述步骤1至步骤4的基础上，可进行后续的分类，如：按给定的采样频率在同一个通道多次随机采集脑电信号，并按照步骤1至步骤4的过程以相同的方式获取在同一个通道每次采集的脑电信号的特征向量，将在同一个通道多次随机采集的脑电信号的特征向量构成一个特征矩阵，该特征矩阵中的每行向量代表一次随机采集的脑电信号的特征向量；然后对该特征矩阵进行分类，得到每次随机采集的脑电信号对应的运动想象类型。

又如：针对多个不同通道，按给定的采样频率同时在多个通道采集脑电信号，并按照步骤1至步骤4的过程以相同的方式获取在每个通道采集的脑电信号的特征向量，将同时在多个通道采集的脑电信号的特征向量按序排列构成行向量；在进行多次采集后，得到多个行向量，将多个行向量按序排列构成一个特征矩阵，该特征矩阵中的每行向量代表同时在多个通道单次采集的脑电信号的特征向量按序排列构成的行向量；然后对该特征矩阵进行分类，得到每次采集的脑电信号对应的运动想象类型。

在本实施例中，多处涉及的按序排列在未给出具体的排列方式的情况下，表示在获取同类信息时按同一种自行定义的顺序排列即可，一般情况下可按先后顺序来排列。

令Y＝[9.72 9.38 6.01 7.46 7.23 5.07 7.92 5.73 5.18 5.71]，F为200，那么Y_δ＝[0.7029 1.4001 1.8677 2.4437 3.0069 3.4138 4.0195 4.4626 4.8542 5.2688]，Y_θ＝[0.5803 1.2786 1.8847 2.6101 3.3521 3.9654 4.6916 5.2505 5.6846 6.0344]，Y_α＝[-0.0364 0.3748 1.1926 2.1650 3.2642 4.3742 5.2636 5.9892 6.3732 6.3498]，Y_β＝[-1.2230 -2.5862 -2.9156 -2.1651 -0.4555 1.6578 2.6811 2.5778 1.55360.1579]；并令d等于2时，f_δ＝2.3755，λ_δ＝0.0828，ξ_δ＝7.8459，

f_θ＝5.4966，λ_θ＝-14.5290，ξ_θ＝3.2526，

f_α＝8.7462，λ_α＝-32.9419，ξ_α＝-1.9559，

f_β＝21.3646，λ_β＝2.9961，ξ_β＝-3.1203，

为验证本发明方法的可行性和有效性，对本发明方法进行实验。

利用本发明方法对国际标准脑电信号数据库目前最大的慢皮层电位(SlowCortical Potential，SCP)数据集中的CLA左右手运动想象数据集进行特征提取，将提取的特征矩阵通过adaboost算法进行分类，其分类准确率高达92.7％，与CSP(Common SpatialPattern，共空间模式)算法和稀疏自编码方法等大多数特征提取算法相比，其准确率明显提高了，这充分说明了本发明方法的可行性和有效性。

Claims (5)

1.一种基于ITD模态参数识别的脑电信号特征提取方法，其特征在于包括以下步骤：

步骤1：按给定的采样频率在一个通道一次随机采集脑电信号；然后将脑电信号以向量形式进行表示，记为Y，Y＝[y₁…y_n…y_N]；其中，给定的采样频率为F赫兹，Y的维数为1×N，Y中的每个元素为一个样本值，N表示Y中的样本值的总个数，N≥5，N＝F×t，t表示采集一个通道的一条脑电信号的采集时间，t的单位为秒，n为正整数，n的初始值为1，1≤n≤N，y₁表示Y中的第1个样本值，y_n表示Y中的第n个样本值，y_N表示Y中的第N个样本值，符号“[]”为向量或矩阵的表示符号；

步骤2：根据脑电信号中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号各自的频率范围，设计各自对应的有限脉冲响应带通滤波器；然后利用脑电信号中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号各自对应的有限脉冲响应带通滤波器分别对Y进行滤波处理，对应提取得到Y中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号，以向量形式对应表示为Y_δ、Y_θ、Y_α、Y_β，Y_δ＝[y_δ,1…y_δ,n…y_δ,N]，Y_θ＝[y_θ,1…y_θ,n…y_θ,N]，Y_α＝[y_α,1…y_α,n…y_α,N]，Y_β＝[y_β,1…y_β,n…y_β,N]；其中，脑电信号中的δ波节律信号的频率范围为0.78～3.90赫兹，脑电信号中的θ波节律信号的频率范围为3.91～7.80赫兹，脑电信号中的α波节律信号的频率范围为7.81～13.28赫兹，脑电信号中的β波节律信号的频率范围为13.29～30.47赫兹，Y_δ、Y_θ、Y_α、Y_β的维数均为1×N，Y_δ中的每个元素为一个样本值，y_δ,1表示Y_δ中的第1个样本值，y_δ,n表示Y_δ中的第n个样本值，y_δ,N表示Y_δ中的第N个样本值，Y_θ中的每个元素为一个样本值，y_θ,1表示Y_θ中的第1个样本值，y_θ,n表示Y_θ中的第n个样本值，y_θ,N表示Y_θ中的第N个样本值，Y_α中的每个元素为一个样本值，y_α,1表示Y_α中的第1个样本值，y_α,n表示Y_α中的第n个样本值，y_α,N表示Y_α中的第N个样本值，Y_β中的每个元素为一个样本值，y_β,1表示Y_β中的第1个样本值，y_β,n表示Y_β中的第n个样本值，y_β,N表示Y_β中的第N个样本值；

步骤3：对Y_δ进行ITD模态参数识别，提取得到Y_δ的四个ITD模态参数，将Y_δ的四个ITD模态参数按序排列构成的向量作为Y_δ的特征向量，记为fea_δ；以相同的方法对Y_θ进行ITD模态参数识别，提取得到Y_θ的四个ITD模态参数，将Y_θ的四个ITD模态参数按序排列构成的向量作为Y_θ的特征向量，记为fea_θ；以相同的方法对Y_α进行ITD模态参数识别，提取得到Y_α的四个ITD模态参数，将Y_α的四个ITD模态参数按序排列构成的向量作为Y_α的特征向量，记为fea_α；以相同的方法对Y_β进行ITD模态参数识别，提取得到Y_β的四个ITD模态参数，将Y_β的四个ITD模态参数按序排列构成的向量作为Y_β的特征向量，记为fea_β；其中，fea_δ、fea_θ、fea_α、fea_β的维数均为1×4；

步骤4：将fea_δ、fea_θ、fea_α、fea_β按序排列构成的向量作为Y的特征向量，记为fea，fea＝[fea_δ fea_θ fea_α fea_β]；其中，fea的维数为1×16。

2.根据权利要求1所述的一种基于ITD模态参数识别的脑电信号特征提取方法，其特征在于所述的步骤2的具体过程为：

根据脑电信号中的δ波节律信号的频率范围，设计第一有限脉冲响应带通滤波器，记为H1，H1的采样频率为F赫兹，H1的阶数为10，H1的阻带上边缘频率为0.7赫兹、阻带下边缘频率为4赫兹、通带上边缘频率为0.8赫兹、通带下边缘频率为3.5赫兹；然后利用H1对Y进行滤波处理，提取得到Y中的δ波节律信号；

同样，根据脑电信号中的θ波节律信号的频率范围，设计第二有限脉冲响应带通滤波器，记为H2，H2的采样频率为F赫兹，H2的阶数为10，H2的阻带上边缘频率为3.9赫兹、阻带下边缘频率为8赫兹、通带上边缘频率为4赫兹、通带下边缘频率为7.5赫兹；然后利用H2对Y进行滤波处理，提取得到Y中的θ波节律信号；

根据脑电信号中的α波节律信号的频率范围，设计第三有限脉冲响应带通滤波器，记为H3，H3的采样频率为F赫兹，H3的阶数为10，H3的阻带上边缘频率为7.8赫兹、阻带下边缘频率为13赫兹、通带上边缘频率为7.9赫兹、通带下边缘频率为12赫兹；然后利用H3对Y进行滤波处理，提取得到Y中的α波节律信号；

根据脑电信号中的β波节律信号的频率范围，设计第四有限脉冲响应带通滤波器，记为H4，H4的采样频率为F赫兹，H4的阶数为10，H4的阻带上边缘频率为13赫兹、阻带下边缘频率为31赫兹、通带上边缘频率为14赫兹、通带下边缘频率为30赫兹；然后利用H4对Y进行滤波处理，提取得到Y中的β波节律信号。

3.根据权利要求1或2所述的一种基于ITD模态参数识别的脑电信号特征提取方法，其特征在于所述的步骤3中，fea_δ的获取过程为：

步骤3_1：将Y_δ＝[y_δ,1…y_δ,n…y_δ,N]离散信号看作自由响应振动信号，对Y_δ＝[y_δ,1…y_δ,n…y_δ,N]依次进行正常采样、d延时釆样、2d延时釆样，依次得到正常采样自由响应数据向量、d延时釆样自由响应数据向量、2d延时釆样自由响应数据向量，对应记为X_nor、X_{del_1}、X_{del_2}，X_nor＝[y_δ,1…y_δ,n…y_δ,m]，X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]，X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]；然后将X_nor＝[y_δ,1…y_δ,n…y_δ,m]重新表示为X_nor＝[x_nor,1…x_nor,n…x_nor,m]，将X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]重新表示为X_{del_1}＝[x_{del_1,1}…x_{del_1,n}…x_{del_1,m}]，将X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]重新表示为X_{del_2}＝[x_{del_2,1}…x_{del_2,n}…x_{del_2,m}]；其中，d为正整数，d表示采样间隔，

m为正整数，m表示X_nor、X_{del_1}、X_{del_2}的长度，1≤m≤N-2d，X_nor、X_{del_1}、X_{del_2}的维数均为1×m，y_δ,m、y_δ,d+1、y_δ,d+n、y_δ,d+m、y_δ,2d+1、y_δ,2d+n、y_δ,2d+m对应表示Y_δ中的第m个样本值、第d+1个样本值、第d+n个样本值、第d+m个样本值、第2d+1个样本值、第2d+n个样本值、第2d+m个样本值，x_nor,1、x_nor,n、x_nor,m对应表示X_nor中的第1个元素、第n个元素、第m个元素，x_{del_1,1}、x_{del_1,n}、x_{del_1,m}对应表示X_{del_1}中的第1个元素、第n个元素、第m个元素，x_{del_2,1}、x_{del_2,n}、x_{del_2,m}对应表示X_{del_2}中的第1个元素、第n个元素、第m个元素；

步骤3_2：根据X_nor和X_{del_1}构建第一增广矩阵，记为X_zg1，

并根据X_{del_1}和X_{del_2}构建第二增广矩阵，记为X_zg2，

然后令A＝X_zg2(X_zg1)^-1，并利用QR分解法对A＝X_zg2(X_zg1)^-1进行处理，得到A的特征值矩阵，记为Δ；其中，X_zg1和X_zg2的维数均为2×m，A表示系统特征矩阵，(X_zg1)^-1为X_zg1的伪逆，A的维数为2×2，Δ的维数为2×2；

步骤3_3：根据Δ计算Y_δ的固有频率和衰减系数，对应记为f_δ和λ_δ，

其中，

表示

的虚数部分，

表示

的实数部分，Δ_1,1表示Δ的第1行第1列元素，t'表示在Y_δ＝[y_δ,1…y_δ,n…y_δ,N]中采样d个样本值的时间间隔，

步骤3_4：应用自由响应理论将X_nor描述为：X_nor＝ψ_norE_nor，

其中，ψ_nor表示X_nor的模态矢量，ψ_nor的维数为1×2，E_nor为引入的中间变量，E_nor的维数为2×m，e表示自然基数，μ＝-λ_δ+j2πf_δ，μ^*为μ的共轭，μ^*＝-λ_δ-j2πf_δ，j表示虚数单位，

然后利用最小二乘法求解X_nor＝ψ_norE_nor，得到ψ_nor，ψ_nor＝X_norE_nor ^T(E_norE_nor ^T)^-1，其中，上标“T”表示向量或矩阵的转置，(E_norE_nor ^T)^-1为E_norE_nor ^T的逆；再计算ψ_nor的幅度和相位，对应记为ξ_nor和

ξ_nor＝abs(ψ_nor(1))+abs(ψ_nor(2))，

其中，abs()为取绝对值函数，ψ_nor(1)表示ψ_nor的第1个元素，ψ_nor(2)表示ψ_nor的第2个元素，

表示

的实数部分；最后将ξ_nor作为Y_δ的幅度，重新记为ξ_δ，将

作为Y_δ的相位，重新记为

步骤3_5：f_δ、λ_δ、ξ_δ、

为Y_δ的四个ITD模态参数，将f_δ、λ_δ、ξ_δ、

按序排列构成的向量作为Y_δ的特征向量fea_δ，

所述的步骤3中，fea_θ、fea_α、fea_β的获取为按照步骤3_1至步骤3_5的过程以相同的方式获取，

其中，f_θ、λ_θ、ξ_θ、

对应为Y_θ的固有频率、衰减系数、幅度、相位，f_α、λ_α、ξ_α、

对应为Y_α的固有频率、衰减系数、幅度、相位，f_β、λ_β、ξ_β、

对应为Y_β的固有频率、衰减系数、幅度、相位。

4.根据权利要求3所述的一种基于ITD模态参数识别的脑电信号特征提取方法，其特征在于按给定的采样频率在同一个通道多次随机采集脑电信号，并按照步骤1至步骤4的过程以相同的方式获取在同一个通道每次采集的脑电信号的特征向量，将在同一个通道多次随机采集的脑电信号的特征向量构成一个特征矩阵，该特征矩阵中的每行向量代表一次随机采集的脑电信号的特征向量；然后对该特征矩阵进行分类，得到每次随机采集的脑电信号对应的运动想象类型。

5.根据权利要求4所述的一种基于ITD模态参数识别的脑电信号特征提取方法，其特征在于针对多个不同通道，按给定的采样频率同时在多个通道采集脑电信号，并按照步骤1至步骤4的过程以相同的方式获取在每个通道采集的脑电信号的特征向量，将同时在多个通道采集的脑电信号的特征向量按序排列构成行向量；在进行多次采集后，得到多个行向量，将多个行向量按序排列构成一个特征矩阵，该特征矩阵中的每行向量代表同时在多个通道单次采集的脑电信号的特征向量按序排列构成的行向量；然后对该特征矩阵进行分类，得到每次采集的脑电信号对应的运动想象类型。

Images