脑电信号的最优分数阶傅里叶变换

张相芬, 陀佳萍, 董柳吟, 袁非牛, 罗 阳

(上海师范大学 信息与机电工程学院, 上海200234)

脑电信号是头皮或大脑皮层脑神经细胞生理活动的反应, 它包含了许多病理信息, 可以为某些脑部疾病提供诊断依据[1]. 自19世纪20年代检测到脑电信号以来, 人们对其开展了大量的研究工作, 脑电信号的处理和分析至今仍然是一项十分困难但又非常重要的研究课题[2].

脑电信号的处理主要分为2个步骤[3]:首先是去除眼电伪迹、去直流和干扰,即将非脑电信号成分去除的预处理;其次是对脑电信号进行特征提取、模式识别等分析,以进行诊断分类或者识别.目前,脑电信号的预处理方法主要包括表面拉普拉斯参考法、主成分分析法及独立成分分析法(independent component analysis,ICA)等.表面拉普拉斯参考法可以在一定程度上去除冗余信号,但未考虑不同位置头皮传输特性的差异[4];主成分分析法可用于移除多通道脑电信号中的眼电干扰成分,此方法的缺点是去线性相关的约束不够[5];ICA主要是去除信号的独立成分干扰[6].

常用的脑电信号处理方法包括非线性动力学分析、人工神经网络分析、高阶谱分析以及时频分析等.随着非线性学科的发展,非线性动力学成为分析和处理脑电信号的有效工具.但是,非线性动力学分析方法对初值极其敏感,因此采用非线性动力学分析处理脑电信号存在困难[7].人工神经网络能高速寻找最优化解,但其泛化能力比较弱,较难找到通用模型,且容易陷入局部最优解[8].高阶谱分析在处理非线性脑电信号时能够有效抑制高斯噪声,但在真值附近起伏很大,受噪声干扰影响严重[9].

脑电信号较为复杂,从定量分析的角度来看,采用传统的时域或者频域分析方法往往性能不佳,因此,研究者想到采用时域与频域结合的方法,即时频分析方法.分数阶傅里叶变换是一种典型的时频分析方法,在脑电信号分析中受到了关注[10].脑电信号在最优分数阶傅里叶变换域有很好的能量聚集性[11].本文基于二维峰值搜索算法提出基于最优分数阶傅里叶变换对脑电信号进行分析处理,以为脑电信号特征提取及分类提供一种新的途径.

1 脑电信号及其预处理

1.1 脑电信号

脑电信号非常微弱,正常情况下,脑电信号的幅值在5～100 μV之间,其频率范围为0.5～35.0 Hz[12].按照频带划分,可将脑电信号分为4类[13].

1)δ波,极低的频带:0.5～3.0 Hz,幅值:10～20 μV,出现在睡眠中、人体缺氧或者大脑病变时.

2)θ波,低频带:4～7 Hz,幅值:20～40 μV,正常成年人在清醒状态下不会表现此波,是一种与幻想、睡眠或者回忆状态有关的波段信号.

3)α波,中频带:8～13 Hz,幅值:10～100 μV,是正常人脑电信号的基本节律,普通人的α波节律调节较好,幅度和频率变化都不明显.

4)β波,高频带:14～26 Hz,幅值:5～30 μV,当个体处于情绪激动、紧张或者亢奋的状态时,β波出现较多.

1.2 脑电信号预处理

脑电信号的预处理主要包括以下步骤[14].

1) 信号滤波.信号滤波即对原始的脑电信号进行低通及高通滤波[15].本文采用的脑电信号为MATLAB R2016a的EEGLAB工具箱自带的数据[16].脑电信号的频率范围一般在0.5～35.0 Hz之间,因此滤波时先对信号进行截止频率为0.5 Hz的高通滤波,再对其进行截止频率为35 Hz的低通滤波.

2) 独立成分分析.独立成分分析是从多元(多维)统计数据中寻找潜在因子或成分的一种方法.将独立成分分析应用在脑电信号分析中,可以较好地识别并去掉眼动和其他噪声[17].

3) 眼电伪迹去除.利用Adjust工具进行眼电伪迹半自动去除[18].去除眼电伪迹后的信号波形更加平稳.

4) 叠加平均.即把脑电信号的每一段32个通道的数据进行叠加平均.

2 分数阶傅里叶变换及二维峰值搜索算法

分数阶傅里叶变换是传统傅里叶变换的一种广义形式.分数阶傅里叶变换反映了信号在时域、频域的信息,适用于非平稳信号的处理,它的快速离散算法较为成熟,计算速度快[19].

分数阶傅里叶变换可以有若干种不同的定义,目前常用的是Ozakats从积分变换角度给出的定义[20],即在时间域的函数x(t)的p阶分数阶傅里叶变换是一个线性积分运算,

(1)

其中,Rp是分数阶傅里叶变换算子,核函数kp(t,v)的表达式为

(2)

式(1)可看作信号x(t)的p阶分数阶傅里叶变换,也可看作信号x(t)在α角度下的分数阶傅里叶变换.因此,分数阶傅里叶变换算子可记为Rα,其核函数也可记为kα(t,v).

本文采用的数值计算方法是Ozaktas采样型离散化算法[21],即首先对原始信号进行时域解调,然后利用香农公式对解调后的信号进行插值计算,最后,对分数阶域变量做离散化处理.

二维峰值搜索算法可以实现线性调频信号的检测与参数估计[22].基于二维峰值搜索算法寻找脑电信号的分数阶傅里叶变换的最优分数阶次可概括为:以变换阶次p作为变量进行扫描,对脑电信号进行连续的分数阶傅里叶变换,形成信号在(p,v)平面上的二维分布.在(p,v)平面进行阈值做峰值点的二维搜索,从而实现脑电信号的最优分数阶次参数估计[23].

设信号的时域表达式为

其中:A为脑电信号的幅值;φ0为初始相位,这里设为0;f0为起始频率;k0为脑电信号的调频率,其值为信号带宽与脉冲宽度的比.

信号的分数阶傅里叶变换为

其中,

(5)

满足k0+cot(pπ/2)=0的p为分数阶变换的最优变换阶次,记为p0,并且有

3 结果与分析

我们从脑电信号(共80多个数据段)中挑选出第3(epochs 3)、第7(epochs 7)、第20(epochs 20)、第33(epochs 33)和第52(epochs 52)这5个数据段.叠加平均后的各个数据段信号时域波形分别如图1所示.将叠加平均后的数据即脑电信号的幅值代入式(3),得到脑电信号的时域表达式.

图1 各个数据段叠加平均后的时域波形Fig.1 Time domain waveforms after each segment is superimposed on average

在参数(p,v)平面上,以0.001的步进值对叠加平均后的信号在0

图2为上述5个波段信号的二维峰值搜索及其一维峰值搜索的结果.由图2可见,二维搜索能得出唯一的峰值,而一维搜索有可能出现多个峰值,较难得出最优分数阶.

图2 各波段信号一维和二维峰值搜索对比图Fig.2 Peak search contrast map of 1D and 2D signals in different bands

对选取的各个数据段的脑电信号进行了最优分数阶傅里叶变换与普通阶次傅里叶变换,其结果如图3所示.由图3可以看出,与进行普通阶次变换的脑电信号相比,进行最优分数阶变换的脑电信号的时频域波形较为平稳,表明噪声得到了抑制.

基于最优分数阶傅里叶变换处理后的脑电信号可以进行大脑的生理状态判别.例如,图3a所示为ephochs 3进行最优分数阶傅里叶变换的时频域波形,显然,0～0.8 s之间的波形为θ波,表明此时个体处于睡眠或者回忆状态;0.8～1.2 s之间的波形为α波,表明此时个体处于安静或者闭眼状态;1.2～1.8 s之间的波形为β波,表明个体处于较激动亢奋的状态.而图3b所示的普通阶次变换下的脑电波形较难提取其节律波.

图3 脑电信号最优分数阶变换与普通阶次变换波形对比图

由此可见,本文基于最优分数阶傅里叶变换的脑电信号处理对后续脑电信号的研究具有重要的意义.

4 结 语

脑电信号包含了大量的生理和病理信息,脑电信号可以为某些脑疾病的临床诊断和治疗提供依据.本文的研究着眼于利用二维峰值搜索算法以及最优分数阶傅里叶变换对脑电信号进行处理和分析.

与一维搜索算法相比,二维搜索算法计算得到的最优阶次分数阶傅里叶变换具有更好的能量聚集性.对比普通阶次傅里叶变换的结果可知,基于二维峰值搜索算法进行脑电信号的最优分数阶傅里叶变换,能更好地去除信号的噪声,让信号波形更加平稳,为脑电信号的进一步研究打下了良好的基础.

本文基于二维峰值搜索算法确定分数阶傅里叶变换的最优阶次,并基于最优阶次的分数阶傅里叶变换对脑电信号进行分析.

本研究采用的是EEGLAB自带的脑电信号,信号成分较好,处理相对容易,因此该方法的适应性仍需进一步验证;在使用二维峰值搜索算法寻找变换阶次时,会存在计算速度与精度的矛盾.因此,最优阶次的搜索算法仍需进一步加以改进.