CN103340623B - 一种在高刺激率下提取诱发电位的方法 - Google Patents

Abstract

本发明公开了一种在高刺激率下提取诱发电位的方法，包括：重复输出一个带抖动的刺激序列，并抽样采集被测试者在该刺激序列重复刺激下的脑电信号，获得观测信号；将观测信号按照进行叠加平均计算，抑制观测信号中的随机噪声，得到平均观测信号；对平均观测信号进行快速傅里叶变换，得到平均观测信号的频域信号，计算刺激序列的连续周期傅里叶变换，得到刺激序列的频域信号；计算得到暂态诱发电位的频域信号；进行傅里叶逆变换，得到时域形式的暂态诱发电位信号。本发明极大地降低了噪声干扰成分对信号恢复的影响，可精确而稳定地提取并恢复出暂态诱发电位。本发明可广泛应用于生物医学领域中。

Description

一种在高刺激率下提取诱发电位的方法

技术领域

本发明涉及生物医学信号处理领域，特别是利用连续傅里叶变换提取高刺激率下暂态诱发电位的方法。 

背景技术

诱发电位(evoked potential，EP)是大脑对外界刺激的电活动的综合反应。诱发电位在临床医学、电生理学、心理学和认知科学等多个学科中具有重要意义。目前诱发电位根据刺激方式不同主要分成三种模态：听觉诱发电位(Auditory Evoked Potential，AEP)，视觉诱发电位(Visual Evoked Potential，VEP)和体感诱发电位(Somatosensory Evoked Potential，SEP)。此外，采用磁、激光等其他刺激方式获取诱发反应的方法也在不断发展。 

以听觉上的瞬态刺激为例，对听觉系统采用的click声(典型持续时间0.1ms)刺激可以诱发出数百毫秒的暂态反应。典型AEP信号的幅度一般在0.1到10μV范围内，远远小于幅度约为20-150μV的背景脑电及各种伪迹干扰信号。所以在实际运用中需要多次重复刺激，利用叠加平均技术抑制和AEP无关的干扰成分，从而得到清晰AEP信号。 

用一个脉冲序列表示刺激序列，脉冲出现的时刻代表刺激出现的时刻，相邻刺激脉冲之间的间隔(interstimulus intervals，ISI)用于等待前一个刺激诱发反应的结束，以便开始下一个刺激。一般用刺激率， 即刺激间隔的倒数(1/ISI)表征刺激出现的频度。高刺激率情况下，等待时间不足，相邻的诱发反应会相互重叠形成所谓的稳态反应(steady-state response)。此种情况下，传统的叠加平均技术无法获得完整的暂态AEP信号。 

从信号处理的角度来看，这种高刺激率导致的重叠现象可以用刺激脉冲序列和暂态AEP的循环卷积模型表征。如果刺激序列中ISI在一定的范围内随机变化，即刺激序列是抖动的，并且刺激序列的这种小的抖动不足以导致每个刺激诱发的暂态AEP波形的显著差别，那么可以利用傅立叶变换的方法将时域卷积计算转换到频率域，利用频率域的逆滤波技术重建暂态AEP信号，从而实现反卷积过程。 

目前基于傅立叶变换的反卷积计算都是利用数字计算机采用离散的快速傅立叶变换(fast-fourier transform，FFT)实现针对高刺激率下诱发电位的恢复。这种计算方式要求刺激序列和观测信号使用相同的时间分辨率，或者说在时域中采用相同的抽样率。这样就会出现两种对优化刺激序列不利的情况：(1)如果序列的优化是在这种离散条件下实现的，抽样率决定了所有可能序列的个数，限制了获取最优序列的搜索范围，不利于获得满足噪声抑制条件的刺激序列；(2)如果已经获取了连续时间意义上的理想刺激序列，采用FFT计算时需要近似为给定抽样率下的离散序列。这种近似的离散序列往往又破坏了原序列在频域上的约束。因此，现有技术无法控制噪声干扰的影响，无法实现稳定而精确地提取并恢复出暂态诱发电位，限制了高刺激率技术的临床研究和应用。 

发明内容

为了解决上述的技术问题，本发明的目的是提供一种在高刺激率下提取诱发电位的方法。 

本发明解决其技术问题所采用的技术方案是： 

一种在高刺激率下提取诱发电位的方法，包括： 

S1、重复输出一个带抖动的刺激序列，并抽样采集被测试者在该刺激序列重复刺激下的脑电信号，获得观测信号； 

S2、将观测信号按照每个刺激序列的起始时刻进行叠加平均计算，从而抑制观测信号中的随机噪声，得到平均观测信号； 

S3、对平均观测信号进行快速傅里叶变换，得到平均观测信号的频域信号，同时计算刺激序列的连续周期傅里叶变换，得到刺激序列的频域信号； 

S4、根据下式重建暂态诱发电位的频域信号： 

$\hat{X}\left({\mathrm{kf}}_0\right)=Y\left({\mathrm{kf}}_0\right){S_T}^{-1}\left({\mathrm{kf}}_0\right)$

其中，为重建的暂态诱发电位的频域信号，Y(kf₀)为平均观测信号的频域信号，S_T(kf₀)为刺激序列的频域信号； 

S5、对暂态诱发电位的频域信号进行傅里叶逆变换，得到时域形式的暂态诱发电位信号。 

进一步，所述步骤S1中刺激序列采用周期为T的周期性脉冲串信号，其表达式为： 

$s_T\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\δ}(t-t_i)$

其中，N表示一个周期内刺激脉冲的个数，t_i表示一个周期内第 i个刺激脉冲出现的时刻。 

进一步，所述步骤S3中刺激序列的频域信号为： 

$S_T\left({\mathrm{kf}}_0\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{e}^{-\mathrm{jk}2\mathrm{\π}{f}_0{t}_i}$

其中，f₀为刺激序列的频率分辨率，且f₀＝1/T。 

进一步，所述刺激序列在频率范围为[f_L，f_H]内满足以下约束条件： 

$\left|S_T^{-1}\left({\mathrm{kf}}_0\right)\right|\≤1,{\mathrm{kf}}_0\∈[f_L,f_H].$

进一步，所述刺激序列在一个周期内的刺激脉冲的个数N为13，刺激序列的周期T为1625.8ms。 

本发明的有益效果是：本发明的一种在高刺激率下提取诱发电位的方法，通过采集被测试者在刺激序列重复刺激下的脑电信号后得到观测信号，并对该观测信号进行叠加平均计算，得到平均观测信号，进而对平均观测信号进行快速傅里叶变换，并计算刺激序列的连续傅里叶变换，最后进行频域反卷积计算得到暂态诱发电位的频域信号，从而对其进行傅里叶逆变换，得到相应时间域上的暂态诱发电位信号。本发明通过结合刺激序列的连续傅里叶变换来提取暂态发电位信号，无需对刺激序列抽样成和观测信号一样的离散形式，或者将刺激序列独立抽样成符合误差限制要求的离散信号，这样可以保证选定的刺激序列在计算时仍然满足噪声抑制的约束条件。极大地降低了噪声干扰成分对信号恢复的影响，可精确而稳定地提取并恢复出暂态诱发电位。 

附图说明

下面结合附图和实施例对本发明作进一步说明。 

图1是实施例1中高刺激率下的卷积及反卷积过程的示意图； 

图2是实施例2中刺激序列的抽样对频域性质的影响示意图； 

图3是实施例2中在不同噪声下不同抽样率的AEP信号恢复效果与本方法恢复效果的对比示意图。 

具体实施方式

本发明提供了一种在高刺激率下提取诱发电位的方法，包括： 

S1、重复输出一个带抖动的刺激序列，并抽样采集被测试者在该刺激序列重复刺激下的脑电信号，获得观测信号； 

S2、将观测信号按照每个刺激序列的起始时刻进行叠加平均计算，从而抑制观测信号中的随机噪声，得到平均观测信号； 

S3、对平均观测信号进行快速傅里叶变换，得到平均观测信号的频域信号，同时计算刺激序列的连续周期傅里叶变换，得到刺激序列的频域信号； 

S4、根据下式重建暂态诱发电位的频域信号： 

$\hat{X}\left({\mathrm{kf}}_0\right)=Y\left({\mathrm{kf}}_0\right){S_T}^{-1}\left({\mathrm{kf}}_0\right)$

其中，为重建的暂态诱发电位的频域信号，Y(kf₀)为平均观测信号的频域信号，S_T(kf₀)为刺激序列的频域信号； 

S5、对暂态诱发电位的频域信号进行傅里叶逆变换，得到时域形式的暂态诱发电位信号。 

进一步作为优选的实施方式，所述步骤S1中刺激序列采用周期为T的周期性脉冲串信号，其表达式为： 

$s_T\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\δ}(t-t_i)$

其中，N表示一个周期内刺激脉冲的个数，t_i表示一个周期内第i个刺激脉冲出现的时刻。 

进一步作为优选的实施方式，所述步骤S3中刺激序列的频域信号为： 

$S_T\left({\mathrm{kf}}_0\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{e}^{-\mathrm{jk}2\mathrm{\π}{f}_0{t}_i}$

其中，f₀为刺激序列的频率分辨率，且f₀＝1/T。 

进一步作为优选的实施方式，所述刺激序列在频率范围为[f_L，f_H]内满足以下约束条件： 

$\left|S_T^{-1}\left({\mathrm{kf}}_0\right)\right|\≤1,{\mathrm{kf}}_0\∈[f_L,f_H].$

进一步作为优选的实施方式，所述刺激序列在一个周期内的刺激脉冲的个数N为13，刺激序列的周期T为1625.8ms。 

下面结合附图通过两个例子对本发明做进一步说明。 

实施例1：结合仿真计算说明本发明如何提取暂态AEP(听觉诱发电位) 

采用能够根据自定义的刺激序列输出高时间精度的刺激声并同时记录下刺激输出时间信息的刺激发生器。脑电采集系统采用符合规范的常规系统即可。 

(1)电极放置：根据常规AEP记录方式，在被测试者头部前额正中发际处放置记录电极，同侧乳突放置参考电极，于眉心处放置地电极，电极阻抗小于5KΩ。 

(2)刺激声输出：首先刺激序列满足以下条件： 

刺激序列采用周期为T的周期性脉冲串信号，其表达式为： 

$s_T\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\δ}(t-t_i)$

其中，N表示一个周期内刺激脉冲的个数，t_i表示一个周期内第i个刺激脉冲出现的时刻。 

刺激序列的频域信号为： 

$S_T\left({\mathrm{kf}}_0\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{e}^{-\mathrm{jk}2\mathrm{\π}{f}_0{t}_i}$

其中，f₀为刺激序列的频率分辨率，且f₀＝1/T。 

刺激序列在选取的频率范围[f_L，f_H]内满足以下约束条件： 

$\left|S_T^{-1}\left({\mathrm{kf}}_0\right)\right|\≤1,{\mathrm{kf}}_0\∈[f_L,f_H].$

一个实际的刺激序列如图1(a)所示。该序列包括13个刺激脉冲，平均刺激率为8Hz。将刺激声按选定的刺激序列重复输出，一般选择短声(click)，通过专业耳机(如ER-3型插入式耳机)直接输出给受试者；同时记录刺激声输出的时间，或者通过发送一个同步触发信号至脑电采集系统，由采集系统记录该时间。 

(3)测试的目的就是使测试获得的信号尽可能地接近实际信号，因此，这里假设要测试并获得的AEP信号为图1(b)所示的中潜伏期AEP信号，此信号即为要提取的诱发电位信号x(t)。 

(4)背景噪声：在没有施加刺激情况下按照典型的AEP中潜伏期记录参数采集自发脑电信号作为背景噪声，获得的背景噪声为图1(c)中间所示的曲线e(t)。 

(5)刺激序列和暂态AEP的卷积过程：高刺激率下相邻刺激脉冲引发的暂态AEP出现重叠的现象可视为刺激序列和暂态AEP信号的卷积效应引起的。因此，在没有背景噪声的理想情况下，系统记录的信号如图1(c)上面的曲线所示，该信号为

(6)观测信号的噪声叠加：实际测试中，背景噪声是不可避免地存在的。因此采集到的观测信号可以看作是步骤(5)所得到的卷积信号和步骤(4)所得到的背景噪声的叠加。 

(7)平均去噪：由于实际情况下AEP的幅度和背景噪声相比很小。为了提高观测数据的信噪比，需要根据刺激序列的安排重复循环刺激，并根据叠加原理进行叠加平均处理从而抑制噪声的幅度。把采集到的观测信号按照每个刺激序列的起始时刻进行分段，然后进行叠加平均，从而可抑制信号中的随机噪声。得到的平均观察信号为重叠的诱发电位，可以视为一种准稳态诱发电位，如图1(c)下面的曲线y(t)所示，y(t)的表达式为： 

$y\left(t\right)=s_T\left(t\right)\⊗x\left(t\right)+e\left(t\right),t\∈[0,T].$

(8)傅立叶变换：对步骤(7)得到的平均观测信号进行快速傅里叶变换，得到平均观测信号的频域信号Y(kf₀)，同时计算刺激序列的连续周期傅里叶变换，得到刺激序列的频域信号S_T(kf₀)； 

(9)反卷积计算：根据下式重建暂态诱发电位的频域信号： 

$\hat{X}\left({\mathrm{kf}}_0\right)=Y\left({\mathrm{kf}}_0\right){S_T}^{-1}\left({\mathrm{kf}}_0\right)$

其中，为暂态诱发电位的频域信号； 

然后，对暂态诱发电位的频域信号进行傅里叶逆变换，得 到时域形式的暂态诱发电位信号如图1(d)所示。进行傅里叶逆变换时考虑到计算过程的舍入误差，结果可能是复数，这里仅取实部即可。 

在实际测试中无法避免背景噪声的出现同时也无法准确地计算背景噪声的大小，因此本方法需要尽可能地消除背景噪声的影响。根据步骤(7)中y(t)的表达式，可知，暂态诱发电位的频域信号的计算公式为： 

$\hat{X}\left({\mathrm{kf}}_0\right)=Y\left({\mathrm{kf}}_0\right){S_T}^{-1}\left({\mathrm{kf}}_0\right)+E\left({\mathrm{kf}}_0\right){S_T}^{-1}\left({\mathrm{kf}}_0\right)$

其中，E(kf₀)为背景噪声的频域信号，而根据步骤(2)的描述可知，当刺激序列的频域信号S_T(kf₀)满足以下条件： kf₀∈[f_L，f_H]时，E(kf₀)成分被抑制，上式右边第二项对重建结果的影响将会减小，重建的暂态AEP信号和真实的暂态AEP信号将会很接近。而且本方法还进行了叠加平均处理，进一步消除了背景噪声的影响。因此，当选用刺激序列的频域信号满足约束条件后，本步骤可采用这个表达式来重建暂态诱发电位的频域信号。这样可在很大程度上抑制背景噪声的影响，使得最后提取得到的暂态诱发电位信号较为准确。 

实施例2：通过与现有方法进行对比来说明本发明的效果 

采用一个在连续时间意义上获得的优化的刺激序列，如图2(a)所示，该刺激序列包括13个刺激脉冲，刺激脉冲出现的时间t_i＝{0，110.6，248.2，363.2，478.6，614.8，735.8，867.1，1006.5，1130.1，1240.7，1351.3，1490.8}ms，序列的长度为T＝1625.8ms，选取的频带范围 [f_L，f_H]内频谱的例数|S_c(kf₀)^-1|如图2(d)所示。 

为了说明现有技术通过离散傅里叶变换提取暂态诱发电位带来的误差，我们选择10个抽样率把该序列分别转化为离散序列，具体抽样率为{1，2，3，4，5，7，9，11，15，20}kHz。抽样率为20kHz时的离散序列如图2(c)所示，t_i′表示序列离散后各刺激脉冲所在位置。图2(b)为刺激序列中第四个刺激脉冲离散前后(t_i和t_i′)的对比图，可看出离散后第四个刺激脉冲出现的时刻t₄′相对于离散前时刻t₄约有0.3ms的延后。离散后刺激序列在选取的频带范围[f_L，f_H]内频谱的倒数|S_d(kf₀)^-1|如图2(e)所示。对比图2(d)和图2(e)可看出对刺激序列的抽样离散化导致刺激序列的频谱不再满足约束条件，如图2(e)中圈起来部分所示。 

如前所述，脑电信号包括确定性的诱发电位信号和随机性的背景噪声两个部分。诱发电位信号选用一段典型的中潜伏期AEP信号，背景噪声是在相同记录条件下采集的自发脑电信号。将AEP信号和刺激序列循环卷积并按照不同比例叠加上背景噪声，即可得到不同信噪比(Signal noise rate，SNR)的观察信号。通过计算提取得到的AEP信号与真实信号之间的均方根相对误差来对提取的AEP信号进行性能评估。 

本实施例中，对三种噪声条件下(SNR＝9.5，0.0，-6.0dB)选择不同抽样率恢复暂态AEP的误差进行了计算，如图3所示，图中纵坐标为恢复信号和原始信号的误差值(用百分比表示)，横坐标表示抽样率，横坐标最右边部分的c表示采用本发明提取得到的AEP信 号与原始信号的误差。根据图3可看出在不同信噪比环境下，对连续时间刺激序列抽样，会影响提取得到的诱发电位信号的质量，虽然随着抽样率的提高，误差值会降低，但是其误差值都大于采用本发明来提取诱发电位信号的误差值。 

实施例1及实施例2中对如何提取听觉诱发电位(AEP)信号的过程以及对比效果的描述，同样也适用于视觉诱发电位(VEP)信号和体感诱发电位(SEP)信号。 

以上是对本发明的较佳实施进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可作出种种的等同变形或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围内。 

Claims (4)

1.一种在高刺激率下提取诱发电位的方法，其特征在于，包括：

S1、重复输出一个带抖动的刺激序列，并抽样采集被测试者在该刺激序列重复刺激下的脑电信号，获得观测信号；

刺激序列采用周期为T的周期性脉冲串信号，其表达式为：

$s_T\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\δ}(t-t_i)$

其中，N表示一个周期内刺激脉冲的个数，t_i表示一个周期内第i个刺激脉冲出现的时刻；

S2、将观测信号按照每个刺激序列的起始时刻进行叠加平均计算，从而抑制观测信号中的随机噪声，得到平均观测信号；

S3、对平均观测信号进行快速傅里叶变换，得到平均观测信号的频域信号，同时计算刺激序列的连续周期傅里叶变换，得到刺激序列的频域信号；

S4、根据下式重建暂态诱发电位的频域信号：

$\hat{X}\left(k{f}_0\right)=Y\left(k{f}_0\right){S_T}^{-1}\left(k{f}_0\right)$

其中，为重建的暂态诱发电位的频域信号，Y(kf₀)为平均观测信号的频域信号，S_T(kf₀)为刺激序列的频域信号；

S5、对暂态诱发电位的频域信号进行傅里叶逆变换，得到时域形式的暂态诱发电位信号；

其中，T为刺激序列的周期，f₀为刺激序列的频率分辨率，k为整数。

2.根据权利要求1所述的一种在高刺激率下提取诱发电位的方法，其特征在于，所述步骤S3中刺激序列的频域信号为：

$S_T\left(k{f}_0\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{e}^{-\mathrm{jk}2\mathrm{\π}{f}_0{t}_i}$

其中，f₀为刺激序列的频率分辨率，且f₀＝1/T。

3.根据权利要求2所述的一种在高刺激率下提取诱发电位的方法，其特征在于，所述刺激序列在频率范围为[f_L,f_H]内满足以下约束条件： $\left|S_T^{-1}\left(k{f}_0\right)\right|\≤1,$ kf₀∈[f_L,f_H]。

4.根据权利要求3所述的一种在高刺激率下提取诱发电位的方法，其特征在于，所述刺激序列在一个周期内的刺激脉冲的个数N为13，刺激序列的周期T为1625.8ms。