技术领域
本发明涉及生物医学信号处理技术，尤其涉及一种在快速刺激下 听觉诱发电位记录过程中关于刺激序列的去噪性能判断处理方法及 系统。
背景技术
听觉诱发电位(Auditoryevokedpotential,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)，其具体为：
S ( k Δ f ) = Σ i = 0 N - 1 e - j 2 πkΔft i ]]>
其中，ti＝{t0,t1,……,tN-1}为段序列中刺激声出现的时刻。
进一步，所述的噪音放大性能系数，其具体计算公式为：
γ = ( Σ k = k L k H w k | S ( k Δ f ) | - 2 ) 1 / 2 ]]>
其中，kL表示为正整数k的最小值，kH表示为正整数k的最大值， wk表示为每个频率点k△f对应的归一化权重因子，而对于所述的归 一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从 而计算得到的。
进一步，对于所述的归一化权重因子，其具体计算公式为：
w k = k - α / 2 Σ k i = k L k H k i - α / 2 ]]>
其中，kL表示为正整数k的最小值，kH表示为正整数k的最大值， α表示为斜率参数。
进一步，所述拟合函数的具体表达式为 log10(P)＝-α×log10(f)+C，其中，α表示为斜率参数。
本发明所采用的另一技术方案是：一种刺激序列的去噪性能判 断处理系统，该系统包括：
获取单元，用于获取段序列；
第一计算处理单元，用于根据傅里叶变换的方法，从而计算段序 列的频域形式表达式S(k△f)，其中，k△f表示为频率点，△f＝1/T，T 为段序列的时长，k为正整数；
第二计算处理单元，用于对背景噪声的功率谱在双对数坐标下进 行线性拟合，从而获得拟合函数中的斜率参数；
第三计算处理单元，用于利用段序列的频率范围，从而计算正整 数k的值；
第四计算处理单元，用于利用得到的正整数k的最大值和最小值、 斜率参数以及段序列的频域形式表达式，从而计算得出该段序列的噪 音放大性能系数；
判断处理单元，用于根据计算得到的噪音放大性能系数，从而对 刺激序列的去噪性能进行判断。
进一步，所述段序列的频域形式表达式S(k△f)，其具体为：
S ( k Δ f ) = Σ i = 0 N - 1 e - j 2 πkΔft i ]]>
其中，ti＝{t0,t1,……,tN-1}为段序列中刺激声出现的时刻。
进一步，所述的噪音放大性能系数，其具体计算公式为：
γ = ( Σ k = k L k H w k | S ( k Δ f ) | - 2 ) 1 / 2 ]]>
其中，kL表示为正整数k的最小值，kH表示为正整数k的最大值， wk表示为每个频率点k△f对应的归一化权重因子，而对于所述的归 一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从 而计算得到的。
进一步，对于所述的归一化权重因子，其具体计算公式为：
w k = k - α / 2 Σ k i = k L k H k i - α / 2 ]]>
其中，kL表示为正整数k的最小值，kH表示为正整数k的最大值， α表示为斜率参数。
进一步，所述拟合函数的具体表达式为 log10(P)＝-α×log10(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(t)表示为实际记录到的听觉诱发反应，x(t)表示为瞬态反应，s(t)表示为段序列，表示为循环卷积运算，n(t)表示为加性背景脑电噪音。而利用傅里叶变换，将上式变为频域形式的表达式，具体如下：
Y(k△f)＝X(k△f)S(k△f)+N(k△f)
故在频域，瞬态反应的估计为：
Y ^ ( k Δ f ) = X ( k Δ f ) + N ( k Δ f ) S ( k Δ f ) - 1 ]]>
上式中的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)，其具体为：
S ( k Δ f ) = Σ i = 0 N - 1 e - j 2 πkΔft i ]]>
其中，ti＝{t0,t1,……,tN-1}为段序列中刺激声出现的时刻。
进一步作为优选的实施方式，所述的噪音放大性能系数，其具体 计算公式为：
γ = ( Σ k = k L k H w k | S ( k Δ f ) | - 2 ) 1 / 2 ]]>
其中，kL表示为正整数k的最小值，kH表示为正整数k的最大值， wk表示为每个频率点k△f对应的归一化权重因子，而对于所述的归 一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从 而计算得到的。
进一步作为优选的实施方式，对于所述的归一化权重因子，其具 体计算公式为：
w k = k - α / 2 Σ k i = k L k H k i - α / 2 ]]>
其中，kL表示为正整数k的最小值，kH表示为正整数k的最大值， α表示为斜率参数。
进一步作为优选的实施方式，所述拟合函数的具体表达式为 log10(P)＝-α×log10(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)，其具体为：
S ( k Δ f ) = Σ i = 0 N - 1 e - j 2 πkΔft i ]]>
其中，ti＝{t0,t1,……,tN-1}为段序列中刺激声出现的时刻，k△f表 示为频率点，△f＝1/T，T为段序列的时长，k为正整数；
S3、对背景噪声的功率谱在双对数坐标下进行线性拟合，从而获 得拟合函数中的斜率参数；
而对于这一步骤S3，其具体为：对背景噪声的功率谱在双对数 坐标下进行最小二乘准则下的线性拟合，从而获得拟合函数中的斜率 参数；
对于上述进行线性拟合后得到的拟合函数，其具体表达式为 log10(P)＝-α×log10(f)+C，其中，α即表示为斜率参数；
S4、利用段序列的频率范围，从而计算正整数k的值；
对于步骤S4，其具体为：利用段序列的频率范围[fL,fH]，从而 计算正整数k的值，使之满足条件k△f∈[fL,fH]；
S5、利用得到的正整数k的最大值kH和最小值kL以及得到的斜 率参数α，从而计算每个频率点k△f对应的归一化权重因子，而对于 所述的归一化权重因子，其具体公式为：
w k = k - α / 2 Σ k i = k L k H k i - α / 2 ]]>
S6、利用得到的归一化权重因子和段序列的频域形式表达式，从 而计算得出该段序列的噪音放大性能系数，而所述的噪音放大性能系 数，其具体计算公式为：
γ = ( Σ k = k L k H w k | S ( k Δ f ) | - 2 ) 1 / 2 ]]>
S7、根据计算得到的噪音放大性能系数，从而对与该段序列对应 的刺激序列的去噪性能进行判断，当γ小于1，则表示与该段序列对 应的刺激序列对噪声有抑制效果，而当γ大于1时，则表示与该段序 列对应的刺激序列会放大脑电噪声。
本发明方法的第二具体实施例
在上述第一具体实施例的描述基础上，对上述步骤S3进行进一 步的描述。在本实施中，如图3所示，图3中的(A)描述了单个段 序列时长的背景噪声的时域波形，而图3中的(B)则绘制了单个段 序列时长的噪声在感兴趣的[10,1000]Hz通带内的频域波形，即背景 噪声的功率谱，背景噪声的功率谱在直角坐标系下，其横坐标为频率 f，纵坐标为功率P。从图中看到，背景噪声的绝大部分能量集中在频 带为[10,300]Hz的低频区间，且功率谱在直角坐标系下表现出类似1/f 的分布特点，幅值在低频段较大，随频率增加而减少。进一步，在双 对数坐标尺度下用最小二乘法拟合得到回归直线具体为 log10(P)＝-1.98×log10(f)+10.14，即P＝1010.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)，其具体为：
S ( k Δ f ) = Σ i = 0 N - 1 e - j 2 πkΔft i ]]>
其中，ti＝{t0,t1,……,tN-1}为段序列中刺激声出现的时刻。
进一步作为优选的实施方式，所述的噪音放大性能系数，其具体 计算公式为：
γ = ( Σ k = k L k H w k | S ( k Δ f ) | - 2 ) 1 / 2 ]]>
其中，kL表示为正整数k的最小值，kH表示为正整数k的最大值， wk表示为每个频率点k△f对应的归一化权重因子，而对于所述的归 一化权重因子，其是根据正整数k的最大值和最小值以及斜率参数从 而计算得到的。
进一步作为优选的实施方式，对于所述的归一化权重因子，其具 体计算公式为：
w k = k - α / 2 Σ k i = k L k H k i - α / 2 ]]>
其中，kL表示为正整数k的最小值，kH表示为正整数k的最大值， α表示为斜率参数。
进一步作为优选的实施方式，所述拟合函数的具体表达式为 log10(P)＝-α×log10(f)+C，其中，α表示为斜率参数。
以上是对本发明的较佳实施进行了具体说明，但本发明创造并不 限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提 下还可做作出种种的等同变形或替换，这些等同的变形或替换均包含 在本申请权利要求所限定的范围内。