CN105662338B - 一种刺激序列的去噪性能判断处理方法及系统 - Google Patents

Abstract

本发明公开了一种刺激序列的去噪性能判断处理方法及系统，该系统包括获取单元、第一计算处理单元、第二计算处理单元、第三计算处理单元、第四计算处理单元以及判断处理单元。该方法包括：获取段序列，并且计算得出该段序列的噪音放大性能系数，从而根据计算得到的噪音放大性能系数来对刺激序列的去噪性能进行判断。通过使用本发明能够快速且准确地选取合适的刺激序列，这样不仅能减少恢复AEP信号的失真程度，而且还能避免选取过程中的主观性以及大大提高选取工作的效率和准确度。本发明可广泛应用于听觉诱发电位检测领域中。

Description

一种刺激序列的去噪性能判断处理方法及系统

技术领域

本发明涉及生物医学信号处理技术，尤其涉及一种在快速刺激下听觉诱发电位记录过程中关于刺激序列的去噪性能判断处理方法及系统。

背景技术

听觉诱发电位(Auditory evoked potential,AEP)为大脑对声音刺激所产生的电势变化，其反映了听觉神经通路和相关中枢神经的活动情况，而对听觉诱发电位进行记录分析的设备则通常被称为听觉诱发电位仪。在进行听觉诱发电位记录分析时，一般会采用短声或短音等持续时间很短的瞬态声音信号来实现声音刺激，然而，常见的AEP的幅度十分微弱，通常小于1μV，远小于背景数十毫伏的脑电和噪声的电位水平，因此，需要采用重复性刺激，通过叠加平均技术来提高信噪比。而对于刺激声的重复率，其则受到AEP持续时间的限制而不能过高，以免造成每个刺激单元诱发反应的重叠。

目前，为了打破刺激重复率的限制，专家们提出了一种可有效消除诱发电位重叠而带来的失真的技术方案。对于这一技术方案，其实现原理为：采用一个特别的刺激串单元来代替单个刺激并周期性重复，其中，所述的刺激串单元包含了同样的瞬态刺激声，而瞬态刺激声之间的间隔是不相等的，并且可比常规重复率限定的短，在一定的信号采样率情况下，可以将刺激串单元用一个二值信号表示，即用“1”表示刺激声的出现时刻，而刺激声之间的间隔则用“0”填充，这样一个二值信号串称为一个周期段的刺激序列，简称段序列；通过采用这样的表示方法，则可将所述刺激序列，即段序列，对应的脑电反应表示为该刺激序列和需要恢复的AEP信号的循环卷积，而由于这种卷积模型的求解计算上可在频域实现，因此则可将该刺激序列的傅里叶变换的倒数作为一个实现去卷积功能的恢复滤波器，这样将采集到的AEP信号经所述恢复滤波器处理后便能恢复瞬态的AEP信号，而该恢复得到的AEP信号会更加贴近实际，大大减少失真的程度。然而，由于所述刺激序列的傅里叶变换数值的模在某些频率点很小，因此会导致恢复滤波器的模在某些频率点上可能异常增大，从而过度强化了该频率成分的噪声能量，这样当通过傅里叶反变换恢复时域AEP信号时则会产生严重的失真。因此由此可得，采用这种去卷积滤波的方法的前提是需要选取一个合适的刺激序列，这样才能有效地抑制噪音，从而更有效地减少AEP信号的失真。

发明内容

为了解决上述的技术问题，本发明的目的是提供一种刺激序列的去噪性能判断处理方法。

本发明的另一目的是提供一种刺激序列的去噪性能判断处理系统。

本发明所采用的技术方案是：一种刺激序列的去噪性能判断处理方法，该方法包括：

获取段序列；

根据傅里叶变换的方法，从而计算段序列的频域形式表达式S(k△f)，其中，k△f表示为频率点，△f＝1/T，T为段序列的时长，k为正整数；

对背景噪声的功率谱在双对数坐标下进行线性拟合，从而获得拟合函数中的斜率参数；

利用段序列的频率范围，从而计算正整数k的值；

利用得到的正整数k的最大值和最小值、斜率参数以及段序列的频域形式表达式，从而计算得出该段序列的噪音放大性能系数；

根据计算得到的噪音放大性能系数，从而对刺激序列的去噪性能进行判断。

进一步，所述段序列的频域形式表达式S(k△f)，其具体为：

其中，t_i＝{t₀,t₁,……,t_N-1}为段序列中刺激声出现的时刻。

进一步，所述的噪音放大性能系数，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，w_k表示为每个频率点k△f对应的归一化权重因子，而对于所述的归一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从而计算得到的。

进一步，对于所述的归一化权重因子，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，α表示为斜率参数。

进一步，所述拟合函数的具体表达式为log₁₀(P)＝-α×log₁₀(f)+C，其中，α表示为斜率参数。

本发明所采用的另一技术方案是：一种刺激序列的去噪性能判断处理系统，该系统包括：

获取单元，用于获取段序列；

第一计算处理单元，用于根据傅里叶变换的方法，从而计算段序列的频域形式表达式S(k△f)，其中，k△f表示为频率点，△f＝1/T，T为段序列的时长，k为正整数；

第二计算处理单元，用于对背景噪声的功率谱在双对数坐标下进行线性拟合，从而获得拟合函数中的斜率参数；

第三计算处理单元，用于利用段序列的频率范围，从而计算正整数k的值；

第四计算处理单元，用于利用得到的正整数k的最大值和最小值、斜率参数以及段序列的频域形式表达式，从而计算得出该段序列的噪音放大性能系数；

判断处理单元，用于根据计算得到的噪音放大性能系数，从而对刺激序列的去噪性能进行判断。

进一步，所述段序列的频域形式表达式S(k△f)，其具体为：

其中，t_i＝{t₀,t₁,……,t_N-1}为段序列中刺激声出现的时刻。

进一步，所述的噪音放大性能系数，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，w_k表示为每个频率点k△f对应的归一化权重因子，而对于所述的归一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从而计算得到的。

进一步，对于所述的归一化权重因子，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，α表示为斜率参数。

进一步，所述拟合函数的具体表达式为log₁₀(P)＝-α×log₁₀(f)+C，其中，α表示为斜率参数。

本发明的有益效果是：本发明的方法是通过计算段序列的噪音放大性能系数从而判断段序列的去噪性能，这样便能根据判断结果从而直接选择出合适的刺激序列，或者进行刺激序列优化工作后得出合适的刺激序列。由此可得，本发明的方法不仅能避免刺激序列选择过程中的主观性，而且还具有步骤简单、易于实现等优点，因此，通过使用本发明的方法便能够快速地对合适的刺激序列进行选取，大大提高选取工作的效率和准确率。

本发明的另一有益效果是：本发明的系统是通过对段序列的噪音放大性能系数进行计算从而判断段序列的去噪性能，这样便能根据判断结果从而直接选择出合适的刺激序列，或者进行刺激序列优化工作后得出合适的刺激序列。由此可得，本发明的系统不仅能避免刺激序列选择过程中的主观性，而且还具有结构简单、易于实现等优点，因此，通过使用本发明的方法便能够快速地对合适的刺激序列进行选取，大大提高选取工作的效率和精确率。

附图说明

下面结合附图对本发明的具体实施方式作进一步说明：

图1是在快速刺激下获取瞬态听觉诱发反应的原理示意图；

图2是本发明一种刺激序列的去噪性能判断处理方法的步骤流程图；

图3是对背景噪声的功率谱在双对数坐标下进行线性拟合的一具体实施例原理示意图；

图4是恢复滤波器的一具体实施例实现原理示意图；

图5是两组恢复瞬态听觉诱发反应的效果对比示意图；

图6是本发明一种刺激序列的去噪性能判断处理系统的结构框图。

具体实施方式

如图1所示，其显示了在快速刺激下获取瞬态听觉诱发电位的试验过程。图1中的(A)为真实的瞬态反应，而图1中的(B)为作用于听觉系统的由“0”和“1”组成的段序列。由于相邻刺激间隔小于瞬态反应时长，会引起前后两个反应之间的重叠，因此这个过程可以用卷积模型来描述，即可得到图1中(C)所示的重叠反应。而考虑到图1中(D)所示的加性背景脑电等噪声干扰伪迹作用，得到实际记录的听觉诱发反应，如图1中的(E)所示，而这个过程可以由下式表述：

其中，y(t)表示为实际记录到的听觉诱发反应，x(t)表示为瞬态反应，s(t)表示为段序列，表示为循环卷积运算，n(t)表示为加性背景脑电噪音。而利用傅里叶变换，将上式变为频域形式的表达式，具体如下：

Y(k△f)＝X(k△f)S(k△f)+N(k△f)

故在频域，瞬态反应的估计为：

上式中的S(k△f)^-1即为图1中(F)所示的由段序列生成的恢复滤波器，因此，实际记录到的听觉诱发反应经过恢复滤波器的去卷积运算后，便得到恢复的瞬态反应，即图1中(G)所示的。由上述可知，恢复滤波器S(k△f)^-1对N(k△f)具有放大或抑制的作用，若段序列在频域出现奇异值，即在个别频率点上幅值过小甚至趋于零，那么对应频率点的噪声则将被异常放大，从而严重影响瞬态反应的恢复质量。因此由此可得，选择一个合适的段序列，这才能有效地抑制噪音，从而更有效地减少恢复后的AEP信号的失真。而本发明的目的则是提供一种刺激序列的去噪性能判断处理方案，这样根据刺激序列的去噪性能评断结果，便能快速且准确地选取一个合适的刺激序列，以实现后续工作的优化。

如图2所示，一种刺激序列的去噪性能判断处理方法，其包括：

获取段序列；

根据傅里叶变换的方法，从而计算段序列的频域形式表达式S(k△f)，其中，k△f表示为频率点，△f＝1/T，T为段序列的时长，k为正整数；

对背景噪声的功率谱在双对数坐标下进行线性拟合，从而获得拟合函数中的斜率参数；

利用段序列的频率范围，从而计算正整数k的值；

利用得到的正整数k的最大值和最小值、斜率参数以及段序列的频域形式表达式，从而计算得出该段序列的噪音放大性能系数；

根据计算得到的噪音放大性能系数，从而对刺激序列的去噪性能进行判断。

进一步作为优选的实施方式，所述段序列的频域形式表达式S(k△f)，其具体为：

其中，t_i＝{t₀,t₁,……,t_N-1}为段序列中刺激声出现的时刻。

进一步作为优选的实施方式，所述的噪音放大性能系数，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，w_k表示为每个频率点k△f对应的归一化权重因子，而对于所述的归一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从而计算得到的。

进一步作为优选的实施方式，对于所述的归一化权重因子，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，α表示为斜率参数。

进一步作为优选的实施方式，所述拟合函数的具体表达式为log₁₀(P)＝-α×log₁₀(f)+C，其中，α表示为斜率参数。

本发明方法的第一具体实施例

一种刺激序列的去噪性能判断处理方法，其具体包括：

S1、获取段序列；

例如，获取的段序列为s(t)＝[1,0,0,...,1,0,0,...,1,0,0,...]，“1”表示刺激声出现的时刻，而相邻刺激声之间的间隔则用“0”表示，其个数并不相等，且代表刺激声出现间隔的长度；并且当信号的采样步长为△t时，则段序列的时长为T＝L△t，其中，L表示为s(t)数据的总数；

S2、由于段序列相当于完整刺激序列的一个周期段，因此，根据傅里叶变换的方法来计算段序列的频域形式表达式S(k△f)，可得到离散的复指数的线性叠加，即对于计算得到的段序列的频域形式表达式S(k△f)，其具体为：

其中，t_i＝{t₀,t₁,……,t_N-1}为段序列中刺激声出现的时刻，k△f表示为频率点，△f＝1/T，T为段序列的时长，k为正整数；

S3、对背景噪声的功率谱在双对数坐标下进行线性拟合，从而获得拟合函数中的斜率参数；

而对于这一步骤S3，其具体为：对背景噪声的功率谱在双对数坐标下进行最小二乘准则下的线性拟合，从而获得拟合函数中的斜率参数；

对于上述进行线性拟合后得到的拟合函数，其具体表达式为log₁₀(P)＝-α×log₁₀(f)+C，其中，α即表示为斜率参数；

S4、利用段序列的频率范围，从而计算正整数k的值；

对于步骤S4，其具体为：利用段序列的频率范围[f_L,f_H]，从而计算正整数k的值，使之满足条件k△f∈[f_L,f_H]；

S5、利用得到的正整数k的最大值k_H和最小值k_L以及得到的斜率参数α，从而计算每个频率点k△f对应的归一化权重因子，而对于所述的归一化权重因子，其具体公式为：

S6、利用得到的归一化权重因子和段序列的频域形式表达式，从而计算得出该段序列的噪音放大性能系数，而所述的噪音放大性能系数，其具体计算公式为：

S7、根据计算得到的噪音放大性能系数，从而对与该段序列对应的刺激序列的去噪性能进行判断，当γ小于1，则表示与该段序列对应的刺激序列对噪声有抑制效果，而当γ大于1时，则表示与该段序列对应的刺激序列会放大脑电噪声。

本发明方法的第二具体实施例

在上述第一具体实施例的描述基础上，对上述步骤S3进行进一步的描述。在本实施中，如图3所示，图3中的(A)描述了单个段序列时长的背景噪声的时域波形，而图3中的(B)则绘制了单个段序列时长的噪声在感兴趣的[10,1000]Hz通带内的频域波形，即背景噪声的功率谱，背景噪声的功率谱在直角坐标系下，其横坐标为频率f，纵坐标为功率P。从图中看到，背景噪声的绝大部分能量集中在频带为[10,300]Hz的低频区间，且功率谱在直角坐标系下表现出类似1/f的分布特点，幅值在低频段较大，随频率增加而减少。进一步，在双对数坐标尺度下用最小二乘法拟合得到回归直线具体为log₁₀(P)＝-1.98×log₁₀(f)+10.14，即P＝10^10.14/f^-1.98，如图3中的(C)所示。由此实施例可反映出，在获取脑电信号过程中，以自发脑电为主的背景噪声的功率谱呈现近似1/f的分布特点。

本发明方法的第三实施例

在本实施例中，该刺激序列在周期为204.8ms的段序列中包括8个刺激脉冲，出现时间为{0，27.2，64.0，100.8，121.6，153.6，172.8，188.8}ms。图4中的(B)表示为该刺激序列对应的恢复滤波器S(k△f)^-1，在通带[10，300]Hz范围内的频谱，(B)中显示了恢复滤波器的频谱峰值位于170Hz附近；图4中的(A)，其表示为，当α＝2时，每个频率点所对应的归一化权重因子，由图可以看出其大小与频率成反比，随频率的增加而降低；图4中的(C)则是加权运算后恢复滤波器的频谱，与(B)所示的频率相比，加权后恢复滤波器的频谱更突出了在低频段的权重。

本发明方法的第四实施例

如图5所示，本实施例所采用的刺激段序列分别为S1和S2，周期均为204.8ms，瞬态刺激声出现时间分别为{0,27.2,64.0,100.8,121.6,153.6,172.8,188.8}ms和{0,16.0,35.2,51.2,72.0,99.2,131.2,168.0}ms，对应的噪声放大指标分别为0.7474和1.9358。仿真测试中，S1和S2刺激诱发的稳态反应如图5中(A)所示。图5(B)所显示的虚线，其表示为原始的瞬态反应，而实线则表示为恢复的瞬态反应。图中清楚显示了，采用不同刺激序列，其信号恢复的效果存在很大差异，如图所示，噪音放大性能系数γ较小的刺激序列S1的恢复效果较好，而且对比原始的瞬态反应和恢复的瞬态反应波形可见，主要特征波V，No，Po和Na基本不存在失真情况，而Pa虽然幅度有所不同，但波形保持良好。但是，噪音放大性能系数γ较大的刺激序列S2的恢复效果则较差，原始的瞬态反应和恢复的瞬态反应存在非常明显的差异，并导致Pa波缺失。因此由此可清楚地表明了噪音放大性能系数γ能很好的反映刺激序列对瞬态信号的恢复效果。

综上所述，通过采用本发明便能快速且准确地选取合适的段序列来减少恢复的AEP信号的失真程度，而且，还能避免选择过程中的主观性，因此，能够大大提高工作效率和准确度。另外，本发明的方法还具有步骤简单、易于实现、便于后续优化工作实现等优点。

上述方法实施例中的描述均适用于以下系统实施例中。

如图6所示，一种刺激序列的去噪性能判断处理系统，该系统包括：

获取单元，用于获取段序列；

第一计算处理单元，用于根据傅里叶变换的方法，从而计算段序列的频域形式表达式S(k△f)，其中，k△f表示为频率点，△f＝1/T，T为段序列的时长，k为正整数；

第二计算处理单元，用于对背景噪声的功率谱在双对数坐标下进行线性拟合，从而获得拟合函数中的斜率参数；

第三计算处理单元，用于利用段序列的频率范围，从而计算正整数k的值；

第四计算处理单元，用于利用得到的正整数k的最大值和最小值、斜率参数以及段序列的频域形式表达式，从而计算得出该段序列的噪音放大性能系数；

判断处理单元，用于根据计算得到的噪音放大性能系数，从而对刺激序列的去噪性能进行判断。

进一步作为优选的实施方式，所述段序列的频域形式表达式S(k△f)，其具体为：

其中，t_i＝{t₀,t₁,……,t_N-1}为段序列中刺激声出现的时刻。

进一步作为优选的实施方式，所述的噪音放大性能系数，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，w_k表示为每个频率点k△f对应的归一化权重因子，而对于所述的归一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从而计算得到的。

进一步作为优选的实施方式，对于所述的归一化权重因子，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，α表示为斜率参数。

进一步作为优选的实施方式，所述拟合函数的具体表达式为log₁₀(P)＝-α×log₁₀(f)+C，其中，α表示为斜率参数。

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

Claims (6)

1.一种刺激序列的去噪性能判断处理方法，其特征在于：该方法包括：

获取段序列；

根据傅里叶变换的方法，从而计算段序列的频域形式表达式S(kΔf)，其中，kΔf表示为频率点，Δf＝1/T，T为段序列的时长，k为正整数；

对背景噪声的功率谱在双对数坐标下进行线性拟合，从而获得拟合函数中的斜率参数；

利用段序列的频率范围，从而计算正整数k的值；

利用得到的正整数k的最大值和最小值、斜率参数以及段序列的频域形式表达式，从而计算得出该段序列的噪音放大性能系数；

所述的噪音放大性能系数，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，w_k表示为每个频率点kΔf对应的归一化权重因子，而对于所述的归一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从而计算得到的；具体地，对于所述的归一化权重因子，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，α表示为斜率参数；

根据计算得到的噪音放大性能系数，从而对刺激序列的去噪性能进行判断。

2.根据权利要求1所述一种刺激序列的去噪性能判断处理方法，其特征在于：所述段序列的频域形式表达式S(kΔf)，其具体为：

其中，t_i＝{t₀,t₁,……,t_N-1}为段序列中刺激声出现的时刻；N为段序列中刺激声出现的时刻的总个数。

3.根据权利要求1或2所述一种刺激序列的去噪性能判断处理方法，其特征在于：所述拟合函数的具体表达式为log₁₀(P)＝-α×log₁₀(f)+C，其中，α表示为斜率参数，f表示为频率，P表示为功率，C表示为常数。

4.一种刺激序列的去噪性能判断处理系统，其特征在于：该系统包括：

获取单元，用于获取段序列；

第一计算处理单元，用于根据傅里叶变换的方法，从而计算段序列的频域形式表达式S(kΔf)，其中，kΔf表示为频率点，Δf＝1/T，T为段序列的时长，k为正整数；

第二计算处理单元，用于对背景噪声的功率谱在双对数坐标下进行线性拟合，从而获得拟合函数中的斜率参数；

第三计算处理单元，用于利用段序列的频率范围，从而计算正整数k的值；

第四计算处理单元，用于利用得到的正整数k的最大值和最小值、斜率参数以及段序列的频域形式表达式，从而计算得出该段序列的噪音放大性能系数；

判断处理单元，用于根据计算得到的噪音放大性能系数，从而对刺激序列的去噪性能进行判断；

所述的噪音放大性能系数，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，w_k表示为每个频率点kΔf对应的归一化权重因子，而对于所述的归一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从而计算得到的；具体地，对于所述的归一化权重因子，其具体计算公式为：

其中，k_L表示为正整数k的最小值，k_H表示为正整数k的最大值，α表示为斜率参数。

5.根据权利要求4所述一种刺激序列的去噪性能判断处理系统，其特征在于：所述段序列的频域形式表达式S(kΔf)，其具体为：

其中，t_i＝{t₀,t₁,……,t_N-1}为段序列中刺激声出现的时刻；N为段序列中刺激声出现的时刻的总个数。

6.根据权利要求4或5所述一种刺激序列的去噪性能判断处理系统，其特征在于：所述拟合函数的具体表达式为log₁₀(P)＝-α×log₁₀(f)+C，其中，α表示为斜率参数，f表示为频率，P表示为功率，C表示为常数。