电生理数据的数字滤波器设计
背景：滤波是脑电图(EEG)和脑磁图(MEG)数据预处理中的一个普遍步骤。滤波除了衰弱噪声信号成分的预期效果之外，滤波还会导致各种意想不到的不利情况例如信号失真和滤波伪影。
方法：对滤波器响应(脉冲响应和频率响应)进行评价，对滤波器类型(高通/低通/带通/带阻；有限/无限脉冲响应，FIR/IIR)的和滤波器参数(截止频率，滤波器阶数和衰减，纹波，延迟和因果关系)选择提供了一些实用的指导，以优化信噪比，避免或减少信号失真。
结果：介绍并讨论了常用电生理软件中各种滤波器的实现方法。对滤波器响应进行比较和评估。
结论：我们提出了识别常见的不良滤波效应和滤波伪影的策略，并在实际例子中演示了它们。本文讨论了通过参数的选择和最佳的实践来限制和筛选替代方案。本文发表在Journal of Neuroscience Methods杂志。(可添加微信号siyingyxf或18983979082获取原文，另思影提供免费文献下载服务，如需要也可添加此微信号入群)
1.介绍
在脑电和脑磁数据的预处理中，滤波是必不可少的一步。滤波可能会严重地改变信号的外观，从而影响得到的结果，但是这是这个过程本身的性质。因此，有些人会选择不滤波或者谨慎的滤波。
在这篇文章中，我们主张使用滤波器，同时注意副作用。滤波是提高电生理数据信噪比的一个非常有用和必要的工具，使许多分析成为可能。滤波器提高了信噪比，但也引入了信号失真。重要的是，在许多情况下，这些信号失真不一定是滤波本身的结果，而是滤波器设计不好的结果。我们已经确定了滤波不大好的两个原因：
(a)缺乏对滤波参数含义的理解;
(b)没有调查数据滤波的后果。
总之，研究人员往往不知道使用滤波器所引起的信号的变化。因此，滤波并没有得到应有的作用(考虑到可能的失真)，因为它只是数据预处理的“一小部分”。我们强烈建议根据每个应用场景的需要选择滤波器并调整它们的参数，以实现所需的噪声衰减同时不使估计的参数产生偏差，并避免在没有仔细考虑的情况下重复使用以前应用过的滤波器。电生理学的研究人员应该调查滤波器对他们想要结果的影响。我们会鼓励每个人都使用滤波器，但要意识到它们就像锋利的刀子，一个非常有用的工具，但要小心处理。
在这里，我们的目的是规范化电生理数据相关的滤波器参数，并评估滤波器响应。为了保持这一结果使用的广泛性，我们缩小和简化了复杂的频率滤波的一般概念，使之成为电生理学滤波器设计中常应用的相关部分，并尽可能避免使用公式。我们假设对信号的频域表示和傅里叶变换有一个粗略的理解。在第二部分，我们演示了如何调整滤波器参数和突出的可能的注意事项。最后，我们研究了常见的信号失真，并介绍了如何识别和处理滤波器失真的启发式方法。一般来说，我们主要关注滤波器设计和滤波器实现的实际方面。
2.第1部分：滤波器设计
时域滤波或频域滤波(相对于空间滤波和其他类型的滤波)是指对特定频率(频带)的信号成分进行衰减。一般来说，滤波背后的基本原理是衰减信号中的噪声，同时保留感兴趣的信号。通常情况下，在同一频带内，信号成分和噪声成分会有重叠。本文所讨论的时域滤波器不能在同一频带内分离信号和噪声;它们只能减弱目标波段内的所有东西。因此，在选择和设计滤波器时必须注意调整滤波器参数，以实现提高信噪比。然而，要设计完全不改变信号的滤波器是不可能的。相反，通过选择适当的滤波器，根据研究人员的目标来适当改变信号，可以实现提高信号信噪比的作用。
系统性的设计和研究滤波器的方法是将每个滤波器看作一个双端口网络中的一个元件，滤波器就是一个有两个输入和两个输出的黑盒。本文不进行详细讨论，只是简单地说，这是一个强大的概念，用于分析模拟电子或信号处理中的复杂情况。这个黑盒子为我们提供了一种分析滤波器属性的标准化方法：将测试信号输入并对输出进行评估。所有的时域滤波器都可以用这种方法来研究。二端口网络模型的另一个有趣之处是，更复杂的滤波器(如带通滤波器)可以很容易地构建成一个低通滤波器和高通滤波器的链。
滤波器响应
按照惯例，用一个非常尖锐的脉冲信号作为测试信号来测试滤波器。这个输入所得到的滤波器响应就是所谓的脉冲响应。频率响应是脉冲响应经过傅里叶变换之后得到的结果，由幅值响应和相位响应两部分组成。所有这些响应都可以用来表征滤波器的特性。脉冲和频率响应则描述了滤波器在时域或频域的传递函数。也就是说，它们描述了一个滤波器对输入信号的影响，最终产生了一个经过过滤后的输出。因此，了解滤波器的响应对于良好的滤波器设计是至关重要的。
图1以线性相位低通滤波器为例显示了一组响应函数。时域响应给出了滤波器是如何改变信号的，它们清楚地显示了滤波器导致的延迟。频域响应显示了滤波器频谱分解的衰减细节，有助于评估和适当调整滤波器参数。图1A和B为时域响应的示例结果，输入分别为脉冲信号和阶跃信号。脉冲和阶跃响应反映了经过滤波后的脉冲或阶跃信号的输出。图中，（1）为截止频率，（2）为过渡带，（3）为通带，（4）为阻带，（5）为通带纹波，（6）为阻带衰减，（7）为滤波器延迟，（8）为信号失真。
如图1A所示，在时域中，用滤波器的脉冲响应来描述滤波器。图中的横坐标为样本数（或以秒为单位的时间），纵坐标为幅值。脉冲响应反映了滤波器在过滤脉冲信号后的输出(图1A中的黑色虚线为脉冲信号)。与脉冲响应类似，阶跃响应反映了滤波器在过滤阶跃信号后的输出(例如，一系列的0后跟一系列的1)。图1B显示的就是阶跃响应，黑色虚线为阶跃信号。注意，在电生理学中，直流偏移阶跃信号出现在信号不连续时，这种信号会导致直流滤波器出现伪影(通常出现在一段信号的开始或结束时或记录停顿时)。由于脉冲信号和阶跃信号都具有跨越整个频谱的能量，所以它们可以用于评估滤波器过滤宽频带复杂信号时是否存在明显的滤波器失真。
图2中，A-E为低通滤波器响应，F-J为高通滤波器响应。图A和图B分别显示的是低通滤波器的脉冲响应和阶跃响应，图F和G显示的是高通滤波器的脉冲响应和阶跃响应。图C和D分别是低通滤波器过滤脉冲信号和阶跃信号的幅值响应，图H和I分别是高通滤波器过滤脉冲信号和阶跃信号的幅值响应。图E和J分别是低通滤波器和高通滤波器的相位响应。阶跃响应最好地显示了信号失真现象，低通滤波器得出的结果明显平滑(较低的截止频率更平滑)，高通滤波器得出的结果明显粗糙(较低的截止频率也更粗糙)。具有较长脉冲响应的滤波器有一个更陡的滚减和一个较窄的过渡带。低通最小相位滤波器尽管具有“最小相位”特性，但仍有较大的延迟。这两种最小相位滤波器都会使信号严重失真。
截止频率
截止频率将滤波器的带通和阻带分开，始终位于过渡带中（如图1C所示）。这是在信号处理过程中应用滤波器时很常见的值，但这不足以描述滤波器的特性。截止频率的定义并不相同：-3dB（half-energy）截止常用于IIR滤波器，-6dB（half-amplitude）截止常用于FIR滤波器和双端IIR滤波器。因此，截止频率应始终与使用的定义相一致。为了避免不必要的信号失真，必须选择截止频率，使信号所需要的频谱成分不被衰减，但尽可能多的去除噪声。这可能会使一些特定的滤波器无法在某些ERP/F研究中使用。一些作者反对使用高通滤波器或者将高通滤波器的截止频率设置的特别低（<0.1Hz）；在估计潜伏期时有些作者也反对使用低通滤波器。我们认可他们的观点，这些顾虑是有必要的。但是，从另一方面考虑，假如滤波器真的提高了信噪比而没有造成系统性的偏差，过滤后的数据肯定是要比未过滤的数据更精确。
滚降、过渡带宽和滤波器阶数
如图1C所示，通带与阻带之间带上截止频率的这一段过渡区间为过渡带。FIR滤波器中的−6dB截止频率在过渡频带的中心。过渡频带的边缘分别由超过带通和阻带的纹波的幅值响应定义。过渡频带中幅值响应的斜率就是滚降。狭窄的过渡带就会带来一个斜率较大的滚降，宽过渡带就会得到一个斜率较小的滚降。拥有大斜率滚降的滤波器比斜率小的滚降的滤波器能更好地分离相邻频带的信号和噪声成分。带有窄过渡带或陡峭滚降的滤波器比带有宽过渡带或较缓滚减的滤波器具有更长的脉冲响应。具有陡峭滚降的滤波器容易产生更强的信号失真和伪影。卷积意味着当前滤波器的输出不仅依赖于当前输入，还依赖于过去的输入（因果滤波器，或非因果过滤器的过去和未来输入）。也就是说，脉冲响应越长，从其中计算当前输出的输入数据的范围就越广。这通常被解释为滤波器在频域的精度与在时域的精度成反比。在可能的情况下，较短脉冲响应的滤波器和较宽的过渡带是最好的。这是一个重要的论点，因为上述原因，有人反对使用带阻滤波器和高通滤波器，因为这两种滤波器往往需要一个非常陡峭的滚减。但是另一方面，由于低频噪声通常是电生理数据中最强的噪声来源，应用高通滤波器很可能会显著提高信噪比。然而，高通滤波器又会导致明显的信号失真，这就必须进一步讨论。
如图1B中(8)所示，滤波器的伪影出现在阶跃信号跃升的瞬态。为了避免振铃效应，过滤数据不应该出现不连续和直流偏移的情况。此外，信号必须适当地填充在信号边缘进行滤波。所需的数据填充量取决于滤波器阶数，这意味着滤波应该最好对连续而不是epoch的数据进行。如果对于特定的目标，可以使用高斯滤波器或贝塞尔滤波器来避免或减少振铃效应。
通带纹波/阻带衰减
实际实现的幅值响应通常偏离要求的幅值响应。这种偏差通常称为通带的通带纹波和阻带的阻带衰减，图1C和D可以看出。通带纹波是线性或对数单位的最大通带偏差。例如，当通带偏差为0.01时，滤波器输出在通带内不会放大或衰减超过1%的信号（在MATLAB中，通带纹波:rp=20log10((1+0.01)/(1−0.01))=0.174dB）。阻带衰减最常用的是对数单位。当阻带衰减为−60dB(或0.001)时，信号在阻带内衰减为1000倍。在大多数滤波器中，通带纹波和阻带衰减都可以得到很好的控制。但是，越少的通带纹波和越强的阻带衰减需要更长的脉冲响应。例如，通带纹波在0.002-0.001(0.2-0.1%)和阻带衰减为−54到−60dB对于许多ERP/F的应用已经是很合理的值了。对于高振幅低频噪声(近乎直流)，阻带衰减可能需要更大的，到100dB左右。由于数值很小，阻带衰减最好的评价方式是对数比例的幅度响应(图1D)，而通带纹波更好的评价方式是线性比例的幅度响应（图1C）。
延迟
每个滤波器的输入和输出之间都必然存在延迟。电生理应用中最相关的参数是群延迟（group delay），反映了各频率分量通过滤波器后在时间上的延迟情况。假如滤波器对各频率分量的延迟相同，就是线性相位滤波器。因此，具有通频带中所有频谱成分的信号不会改变其时间形状。线性相位滤波器具有完美对称的脉冲响应。所谓的具有非对称脉冲响应的非线性相位滤波器会在不同的频带引入不同的延迟。因此，即使所有的频谱成分都在通带内，非线性相位滤波器也会扭曲宽频带信号对如ERP造成影响，而且当分析时频相位或相幅耦合时，它们会干扰跨频相位关系。
线性相位滤波器的延迟可以通过将滤波器输出及时移回来进行校正，从而得到无延迟的零相位滤波器。由于进行了移位，经过滤波的输出信号中的每个结果都是从未经过滤波的输入信号的前(过去)和后(未来)数据计算出来的，因此，该过滤器被归类为非因果性的。在实践中，这就意味着零相位滤波器输出的信号可能在输入信号开始前就已经偏离了基线，可能系统地低估了低通滤波后的延迟，在高通滤波后会出现人为改变的成分，或者把刺激后的振荡放置在刺激前的区间中导致对刺激前阶段的错误解释。相比之下，因果滤波器只根据之前(过去)的输入来计算输出。如果需要因果滤波器，应考虑非线性最小相位滤波器，因为它在给定幅度响应的每个频率上只产生最小的可能延迟，但由于非线性而使信号发生失真。因果型高通最小相位(和其他非线性)滤波器引入相当小的延迟，而因果低通(带通和带阻)滤波器即使在最小相位特性下也会引入更大的延迟，这就是为什么在电生理学中不推荐它们的原因。
零相位延迟也可以用非线性滤波器实现，通过对滤波器输出的二次反向滤波(two-pass fifiltering)来补偿滤波器延迟。正反向滤波的结果是非因果滤波器会具有对称的脉冲响应。双通滤波(相当于在双端口模型中将同一个滤波器串联两次)使滤波器的阶数加倍，(有效)脉冲响应的长度加倍。双通滤波器还会增强通带波纹和阻带衰减。不同的软件会采用不同的策略来补偿滤波器的倍频和移位的截止频率。因此，为了可复制性，不仅要确定截止频率的定义，而且要确定截止频率是否适用于单通或双通滤波，包括可能的阶数和截止频率调整。重要的是，双通滤波器和等效双通滤波器(具有相同的有效脉冲响应长度)的幅值响应的形状可以显著不同，特别是在截止频率附近的频带。与相应的双通滤波器相比，单通线性相位滤波器(通过移位校正)可以获得相似的幅度响应。这使得线性相位单通滤波在许多应用中更可取。
IIR滤波器与FIR滤波器
数字滤波器有无限脉冲响应(IIR)或有限脉冲响应(FIR)两种类型。IIR滤波器具有非对称脉冲响应和非线性相位，截止频率通常为−3dB。IIR滤波应具有双重精度，并经常检查IIR滤波器的稳定性，特别是当滤波器具有极端截止频率时，例如在ERP/F数据分析的高通滤波中。在电生理学中通常使用Butterworth滤波器和椭圆IIR滤波器。Butterworth滤波器没有通带纹波和阻带纹波，相比于其他常用的滤波器例如椭圆IIR滤波器和Chebyshev滤波器在截止频率附近滚降较小。椭圆IIR滤波器有通带纹波和阻带纹波，在截止频率附近有最陡的滚降。
FIR滤波器可以具有（反）对称脉冲响应或非对称脉冲响应、（反）对称线性相位响应或非对称非线性相位响应。截止频率通常被指定为−6dB。脉冲响应和滤波器系数相同，FIR滤波器的滤波器长度就是脉冲响应的长度。在电生理学中几乎只使用奇数长度的对称的FIR滤波器(只有奇数长度的FIR滤波器可以通过左移校正到零相位延迟）。使用sinc窗口函数的FIR滤波器是基于正弦函数近似矩形幅度响应，因此，有时被称为“理想”滤波器。对于有限阶滤波器，脉冲响应必须通过加窗函数来减少通带纹波和阻带纹波。过渡带宽与滤波器阶数(滤波器长度-1)和窗口类型有关。研究中所要求的过渡带宽以及通带纹波和阻带纹波决定了滤波器的阶数。等纹波滤波器也被称为“最优”滤波器，因为它们对于给定的参数具有最小的阶数。
IIR滤波器和FIR滤波器的滤波器阶数无法比较。相反，必须比较产生的脉冲响应长度。尽管IIR滤波器通常被认为是计算效率更高的，但它们只在需要高吞量和快速截止时被推荐。对于离线数据分析，吞吐量和计算时间在现代计算机硬件上并不重要。因此，当需要因果滤波器时，应该考虑IIR滤波器。在某些特定情况下，因果滤波器是可取的。因果滤波器只会将早期的影响放到后期的延迟中。另一方面，由于数值精度原因，IIR滤波器的有效脉冲响应也是有限的，因此，所有相关特性也可以用FIR滤波器实现。综合起来，FIR滤波器更容易控制，总是稳定的，有一个明确的通带，可以纠正到零相位无需额外的计算，并可以转换为最小相位。因此，我们推荐在电生理数据分析中使用FIR滤波器。
第2部分：电生理学软件中的滤波器
EEGLAB 13 (MATLAB toolbox)
13版EEGlab的滤波器默认使用基本FIR滤波器（新），在EEGlab插件的基础上实现零相位汉明窗sinc FIR滤波器。基本FIR滤波器的通带边缘是指定的，而不是截止频率。如果未指定滤波器阶数，则滤波器阶数为默认值，选择低通带边缘的25%为过渡带宽，但在可能的情况下，不低于2Hz，否则过渡频带选择为通带边缘到临界频率的距离。我们建议直接使用加窗的sinc FIR滤波器，特别是有效的Kaiser窗口可以自由调整通带纹波/阻带衰减。截止频率(定义为半幅值，即−6dB)和滤波器阶数必须是指定好的。用基本的和加窗的sinc滤波器实现的非因果零相位滤波器也可任意转换为因果非线性最小相位滤波器。EEGlab还提供了一个零相位“Parks-McClellan (equiripple)”滤波器。对于这个滤波器，−6dB的截止频率和过渡带宽也是指定好的。频带权值和滤波器阶数可以根据过渡带宽、通带纹波和阻带纹波来指定或估计。这两种滤波器都包含一个可视化的滤波器响应。
还可以下载插件实现椭圆IIR滤波器。可以实现零相位(正向和反向滤波)和因果非线性相位(正向滤波)滤波器。在默认的基本FIR滤波器中，必须指定通带边缘而不是截止频率。在默认情况下，过渡带宽为1Hz，也可以在用户界面中指定过渡带宽。通带纹波默认为0.0025dB，阻带衰减为−40dB。纹波只能在命令行中编写才能调整为用户定义的值。带通(和带阻)滤波器的实现方法为单独的高通和低通滤波器的组合。
ERPlab（EEGlab的插件）
ERPlab是一个EEGlab中的插件，有自己的滤波器例程。IIR Butterworth、
汉明窗sinc FIR和Parks-McClellan FIR(陷波器)滤波器都可在此插件上实现为双通正向和反向非因果零相位滤波器。IIR滤波器和FIR滤波器的-6dB截止频率和滤波器阶数都是指定的。为了补偿双通滤波，滤波器的阶数被调整为IIR和FIR滤波器内部报告阶数的一半。对于FIR滤波器，截止频率经过隐性调整，经过两路滤波后，在指定截止频率处达到−6dB的衰减。FIR滤波器的阶数被限制为最大4096，因此它设计出来的高通和带通滤波器的截止频率很低,经常用于ERP/F分析。ERPlab中的滤波器不能指定过渡带宽。
Fieldtrip
Fieldtrip中提供了IIR Butterworth滤波器, 汉明窗sinc FIR滤波器和“firls”FIR滤波器这几种滤波器。Fieldtrip中默认的滤波器是一个六阶的Butterworth滤波器(高通，低通)或四阶(带通，带阻)的。FIR滤波器的默认阶数计算方法为3×fix(fs/fc)，fc为较低的截止频率。Fieldtrip中可以使用单通正向滤波器或单通反向滤波器，但不能对滤波器延迟进行校正。该滤波器被拟合为矩形频域函数。由于默认滤波器阶数过低，无法近似矩形频率响应，拟合可能会导致滤波器振铃效应过大、通带纹波过大、直流增益不统一等不利影响。我们强烈建议不要使用firls滤波器。双通滤波后的滤波器阶数和−3dB(IIR)或−6dB(FIR)截止频率是不知道，而且滤波器响应不是可视化的。过渡带宽不能指定，滤波器的阶数需要手动估计。在提交本稿件的修订版时，第一作者将EEGlab的sinc FIR滤波器移植并集成到Fieldtrip。Fieldtrip的一个即将发布的版本将允许使用者控制通带纹波和阻带衰减，根据过渡带宽估计滤波器阶数，以及使用firws滤波器时可以绘制滤波器响应。
BrainVision Analyzer
BrainVision Analyzer中提供了2、4或8阶零相位IIR Butterworth滤波器。-3dB的截止频率是指定的。零相位（从输入端来看的）是假定使用了一个半阶数的滤波器两次，即正向和反向。所应用的截止频率被调整以补偿双通滤波，并在指定的截止频率上保持了−3dB的衰减。使用的双通滤波和隐式调整截止频率和滤波器阶数的方法是没有明确的记载，这使得其滤波器难以用其他软件复制。
EEProbe
EEProbe提供了非因果零相位带窗口(Hamming, Hann，Blackman, Bartlett, Tukey和矩形窗口)的sinc和ParksMcClellan FIR滤波器。这些都是由xfir软件设计的。零相位是通过校正滤波器输出的延迟来实现的。线性和对数比例的幅值响应和阶跃响应可以显示。对于加窗sinc FIR滤波器，过渡带宽不能指定，也不报告，即必须估计所需的滤波器阶数。对于ParksMcClellan滤波器，可以指定过渡带宽、通带纹波和阻带衰减，并自动计算滤波器阶数。
图3显示的是ERPlab、Fieldtrip、BrainVision Analyzer实现的30Hz(采样频率fs=500hz)低通二阶IIR Butterworth滤波器的幅值响应。ERPLAB通过调整滤波器阶数来补偿双通的正向滤波和反向滤波，BrainVision Analyzer通过调整滤波器阶数和截止频率来补偿，Fieldtrip不补偿双通滤波(这意味着截止频率和阶数指的是单通滤波器)。为了比较，二阶IIR Butterworth滤波器的幅值响应为图中的灰色虚线。请注意，0-30Hz为带通频带。尽管滤波器参数相同，但这三种不同的实现会导致显著不同的幅度响应。
数据处理业务介绍：
图７显示了ERP数据中的高通滤波失真。我们可以看到0.75Hz高通滤波器(−6dB截止，零相位1100阶汉明窗sinc FIR滤波器)对标准音调产生的ERP影响很小。在异常音调产生的ERP中，此滤波器降低了P3的幅值，并人为地提高了N1/MMN和N2的幅值(相对于刺激前的基线)。而且，1.5Hz高频滤波器（参数一致）的失真明显更大。图G显示了滤波前后N2成分（190-210ms）的脑地形图和滤波前原始信号减去过滤后信号的脑地形图。需要注意的是，高通滤波器不仅会扭曲时间动力学，还会扭曲成分地形(相对于刺激前基线)。因果最小相位滤波器保留了N1/MMN和N2成分的地形和时间动态，但由于非线性效应而扭曲了P3成分的形态。