注意
转到末尾下载完整示例代码。
滤波背景知识#
这里我们介绍一些关于滤波的通用背景知识，以及 MNE-Python 中具体是如何实现的。关于数字滤波器设计的实际应用，推荐阅读 Parks & Burrus (1987) [ 1 ] 和 Ifeachor & Jervis (2002) [ 2 ]；关于 M/EEG 场景下的滤波，我们推荐阅读 Widmann et al. (2015) [ 3 ]。
注意
本教程相当深入地探讨了滤波的数学原理，以及选择滤波器时涉及的设计决策。如果你只想知道如何在 MNE-Python 中把默认滤波器应用到数据上，可以跳过本教程，改读滤波与重采样教程（不过总有一天，你还是应该回来读一读这篇教程 🙂）。
问题描述#
关于对电生理数据进行滤波的实际问题，Widmann et al. (2012) [ 4 ] 中有详细讨论，他们的结论如下：
正如 VanRullen (2011) [ 5 ] 所展示的，滤波可能使信号的时间进程（以及幅值）产生相当大的畸变。因此，不应轻率地使用滤波。不过，如果谨慎考虑滤波的影响并将滤波伪迹降到最低，仍可对滤波后电生理数据的时间动态做出有效解释，而且滤波还能检测出原本可能被遗漏的信号。
换言之，滤波可以提高信噪比（SNR），但如果不谨慎使用，也可能使数据发生畸变。我们希望在这里介绍一些滤波基础知识，使用户能更好地理解滤波中的权衡，以及 MNE-Python 为何选择某些特定的默认设置。
滤波基础#
先来了解一些基础数学。在频域中，数字滤波器具有一个传递函数，其形式为：
在时域中，可以用分子系数 \(b_k\) 和分母系数 \(a_k\)，由输入数据 \(x(n)\) 得到输出数据 \(y(n)\)，即：
换言之，时刻 \(n\) 的输出由对以下各项的求和决定：
分子系数 \(b_k\)，它们与之前的输入值 \(x(n-k)\) 相乘；以及
分母系数 \(a_k\)，它们与之前的输出值 \(y(n-k)\) 相乘。
注意，这两个求和分别对应 (1) 加权滑动平均和 (2) 自回归。
根据这些系数，滤波器分为两类：FIR（有限脉冲响应）和 IIR（无限脉冲响应）滤波器。FIR 滤波器使用有限个分子系数 \(b_k\)（其中 \(a_k=0\qquad \forall k\)），因此每个输出值 \(y(n)\) 只取决于 \(M\) 个先前的输入值。IIR 滤波器则依赖先前的输入和输出值，因此其脉冲响应实际上可以无限长。
如 Parks & Burrus (1987) [ 1 ] 所述，FIR 和 IIR 滤波器各有利弊：
因果 FIR 滤波器可以是线性相位的——即在所有频率上具有相同的时间延迟——而因果 IIR 滤波器则不能。FIR 滤波器的相位和群延迟特性通常也更好。
一般而言，在相同阶数下，IIR 滤波器可以具有比 FIR 滤波器更陡峭的截止特性。
IIR 滤波器通常数值稳定性较差，部分原因是（递归计算导致的）误差累积。
在 MNE-Python 中，我们默认使用 FIR 滤波。正如 Widmann et al. (2015) [ 3 ] 所指出的：
尽管 IIR 滤波器常被认为计算效率更高，但只有在需要高吞吐量和陡峭截止时才推荐使用（Ifeachor and Jervis, 2002 [ 2 ]，第 321 页）。FIR 滤波器更易于控制，始终稳定，具有明确的通带，无需额外计算即可校正为零相位，并且可以转换为最小相位。因此，在电生理数据分析的大多数场合，我们推荐使用 FIR 滤波器。
设计滤波器（FIR 或 IIR）时，总有一些权衡需要考虑，包括但不限于：
通带内的纹波
阻带的衰减
滚降的陡峭程度
滤波器阶数（即 FIR 滤波器的长度）
时域振铃
一般来说，在频率上越尖锐，在时间上就越宽，反之亦然。这是基本的时频权衡，下面的内容会体现这一点。
FIR 滤波器#
首先，我们将重点介绍 FIR 滤波器，这是 MNE-Python 默认使用的滤波器。
设计 FIR 滤波器#
这里我们将尝试设计一个低通滤波器，并考察时域与频域滤波器特性方面的权衡。稍后，在应用 FIR 滤波器一节中，我们将看看这类滤波器在实际使用时会对信号产生什么影响。
首先导入一些有用的滤波工具，并为我们的数据设置一些对 M/EEG 而言合理的默认值。
import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import fft, fftfreq
from scipy import signal
import mne
from mne.time_frequency.tfr import morlet
from mne.viz import plot_filter, plot_ideal_filter
sfreq = 1000.0
f_p = 40.0
# limits for plotting
flim = (1.0, sfreq / 2.0) # frequencies
dlim = (-0.2, 0.2) # delays
以理想低通滤波器为例，它在通带内（直到频率 \(f_p\)）的幅值响应为 1，在阻带内（低至频率 \(f_s\)）的幅值响应为 0，此处 \(f_p=f_s=40\) Hz（为简便起见，图中只显示到 -60 dB 的下限）：
nyq = sfreq / 2.0 # the Nyquist frequency is half our sample rate
freq = [0, f_p, f_p, nyq]
gain = [1, 1, 0, 0]
third_height = np.array(plt.rcParams["figure.figsize"]) * [1, 1.0 / 3.0]
ax = plt.subplots(1, figsize=third_height, layout="constrained")[1]
plot_ideal_filter(freq, gain, ax, title=f"Ideal {f_p} Hz lowpass", flim=flim)
假设上，这种滤波器在频域中具有零纹波、完美衰减和完美陡峭度。然而，由于频率响应存在不连续，要实现该滤波器就需要时域中的无限振铃（即无限阶）。换个角度理解：频域中的矩形窗实际上对应时域中的 sinc 函数，要表示它需要无限多个采样点（因而需要无限长的时间）。因此，尽管这种滤波器具有理想的频率抑制，其时域特性却很差。
我们试着用最朴素的方式构造一个长度为 0.1 s 的砖墙（sinc）滤波器，并观察滤波器本身在时域和频域中的表现：
效果不太好！把滤波器加长 10 倍（1 s）可以略微改善阻带抑制，但时域中仍有大量振铃。注意这里的 x 轴长了一个数量级，滤波器的群延迟也相应长得多（同样等于滤波器长度的一半，即 0.5 秒）：
我们再用更长的滤波器（10 s）让阻带更紧，x 轴也随之更大：
现在频率抑制非常尖锐，但滤波器在整整 10 秒内都在振铃。所以这种朴素方法大概不是构造低通滤波器的好办法。
幸运的是，有多种成熟的方法可以根据期望的响应特性来设计 FIR 滤波器，包括：
窗函数 FIR 设计（
scipy.signal.firwin2()、scipy.signal.firwin()）最小二乘设计（
scipy.signal.firls()）频域设计（在傅里叶域构造滤波器，并用 IFFT 反变换）
注意
当频率响应中存在“无需关心”（do not care）的区域时，Remez 和最小二乘设计具有优势。然而，我们希望在所有频率区域都得到良好控制的响应。当需要任意响应时，频域构造法很好，但对于较直接的滤波器应用，它通常不如窗函数方法干净（由于采样问题）。由于我们的滤波器（低通、高通、带通、带阻）相当简单，且需要精确控制所有频率区域，我们将主要使用并探讨窗函数 FIR 设计。
如果稍微放宽频域滤波器要求，我们就可以用这些函数构造一个低通滤波器，它转而具有一个过渡带，即通带频率 \(f_p\) 与阻带频率 \(f_s\) 之间的区域：
trans_bandwidth = 10 # 10 Hz transition band
f_s = f_p + trans_bandwidth # = 50 Hz
freq = [0.0, f_p, f_s, nyq]
gain = [1.0, 1.0, 0.0, 0.0]
ax = plt.subplots(1, figsize=third_height, layout="constrained")[1]
title = f"{f_p} Hz lowpass with a {trans_bandwidth} Hz transition"
plot_ideal_filter(freq, gain, ax, title=title, flim=flim)
接受频域中较平缓的滚降，可以让我们的时域响应可能好得多。过渡区的斜率变得更缓，但时域信号干净得多。这里仍以 1 s 滤波器为例：
由于我们的低通滤波器在 40 Hz 附近、过渡带宽为 10 Hz，我们其实可以使用更短的滤波器（10 Hz 的 5 个周期 = 0.5 s），仍能获得可接受的阻带衰减：
但如果把滤波器缩得太短（10 Hz 的 2 个周期 = 0.2 s），有效阻带频率就会被推到 60 Hz 以外：
如果我们想要一个只有 0.1 秒长的滤波器，大概应使用接近 25 Hz 的过渡带（0.2 s = 25 Hz 的 5 个周期）：
trans_bandwidth = 25
f_s = f_p + trans_bandwidth
freq = [0, f_p, f_s, nyq]
h = signal.firwin2(n, freq, gain, fs=sfreq)
plot_filter(
h, sfreq, freq, gain, "Windowed 50 Hz transition (0.2 s)", compensate=True, **kwargs
)
到目前为止，我们只讨论了非因果滤波，这意味着时刻 \(t\) 的每个采样点都是利用当前时刻 \(t\) 之后（\(t + \Delta t\)）和之前（\(t - \Delta t\)）的采样点来滤波的。从这个意义上说，每个采样点都受其前后采样点的影响。这在许多情况下都很有用，尤其是因为它不会延迟事件的时间。
然而，有时使用因果滤波会更有益，即每个采样点 \(t\) 只利用在其之后的时间点进行滤波。
注意，这种延迟是可变的（而线性/零相位滤波器为常数），但在通带内很小。零相位滤波器需要把线性相位滤波级的输出向前做时间平移（因而变成非因果），与之不同，最小相位滤波器无需任何补偿即可在通带内获得很小的延迟。
我们可以从现有的线性相位滤波器构造一个最小相位滤波器，注意其下降沿没有那么陡。这里我们使用函数 scipy.signal.minimum_phase(..., half=False) 来完成。
h_min = signal.minimum_phase(h, half=False)
plot_filter(h_min, sfreq, freq, gain, "Minimum-phase", **kwargs)
注意，与未补偿的因果线性相位滤波器相比，这个最小相位滤波器的延迟小得多，如下所示。这就是为什么 MNE-Python 目前只提供零相位非因果滤波器（即延迟补偿滤波器）或最小相位因果滤波器。线性相位滤波器具有恒定且较大的延迟：
应用 FIR 滤波器#
现在让我们把这些滤波器应用到一些数据上，看看它们的实际效果。
我们构造一个高斯窗正弦波（即 Morlet 小波的虚部）加上噪声（随机噪声和工频噪声）。注意，原始干净信号在我们的低通滤波器的通带和过渡带内都含有频率成分。
dur = 10.0
center = 2.0
morlet_freq = f_p
tlim = [center - 0.2, center + 0.2]
tticks = [tlim[0], center, tlim[1]]
flim = [20, 70]
x = np.zeros(int(sfreq * dur) + 1)
blip = morlet(sfreq, [morlet_freq], n_cycles=7)[0].imag / 20.0
n_onset = int(center * sfreq) - len(blip) // 2
x[n_onset : n_onset + len(blip)] += blip
x_orig = x.copy()
rng = np.random.default_rng(0)
x += rng.normal(scale=1e-3, size=len(x))
x += np.sin(2.0 * np.pi * 60.0 * np.arange(len(x)) / sfreq) / 2000.0
我们用截止较平缓的线性相位 FIR 滤波器对它滤波（这样可以补偿恒定的滤波器延迟）：
transition_band = 0.25 * f_p
f_s = f_p + transition_band
freq = [0.0, f_p, f_s, sfreq / 2.0]
gain = [1.0, 1.0, 0.0, 0.0]
# This would be equivalent:
h = mne.filter.create_filter(
x, sfreq, l_freq=None, h_freq=f_p, fir_design="firwin", verbose=True
)
x_v16 = np.convolve(h, x)
# this is the linear->zero phase, causal-to-non-causal conversion / shift
x_v16 = x_v16[len(h) // 2 :]
plot_filter(h, sfreq, freq, gain, "MNE-Python ≥0.16 default", compensate=True, **kwargs)
Setting up low-pass filter at 40 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 331 samples (0.331 s)
现在换一种设计方法 fir_design='firwin2' 来滤波，同样补偿恒定的滤波器延迟。与 fir_design='firwin' 相比，尽管长度是其两倍，这种方法的过渡并没有那么陡：
transition_band = 0.25 * f_p
f_s = f_p + transition_band
freq = [0.0, f_p, f_s, sfreq / 2.0]
gain = [1.0, 1.0, 0.0, 0.0]
# This would be equivalent:
# filter_dur = 6.6 / transition_band # s
# n = int(sfreq * filter_dur)
# h = signal.firwin2(n, freq, gain, nyq=sfreq / 2.)
h = mne.filter.create_filter(
x, sfreq, l_freq=None, h_freq=f_p, fir_design="firwin2", verbose=True
)
x_v14 = np.convolve(h, x)[len(h) // 2 :]
plot_filter(h, sfreq, freq, gain, "MNE-Python 0.14 default", compensate=True, **kwargs)
Setting up low-pass filter at 40 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed frequency-domain design (firwin2) method
- Hamming window
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 661 samples (0.661 s)
我们再用 MNE-Python 0.13 的默认设置滤波，它是一个长时长、陡峭截止的 FIR 滤波器，会被应用两次：
transition_band = 0.5 # Hz
f_s = f_p + transition_band
filter_dur = 10.0 # s
freq = [0.0, f_p, f_s, sfreq / 2.0]
gain = [1.0, 1.0, 0.0, 0.0]
# This would be equivalent
# n = int(sfreq * filter_dur)
# h = signal.firwin2(n, freq, gain, nyq=sfreq / 2.)
h = mne.filter.create_filter(
x,
sfreq,
l_freq=None,
h_freq=f_p,
h_trans_bandwidth=transition_band,
filter_length=f"{filter_dur}s",
fir_design="firwin2",
verbose=True,
)
x_v13 = np.convolve(np.convolve(h, x)[::-1], h)[::-1][len(h) - 1 : -len(h) - 1]
# the effective h is one that is applied to the time-reversed version of itself
h_eff = np.convolve(h, h[::-1])
plot_filter(
h_eff, sfreq, freq, gain, "MNE-Python ≤0.13 default", compensate=True, **kwargs
)
Setting up low-pass filter at 40 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed frequency-domain design (firwin2) method
- Hamming window
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 40.25 Hz)
- Filter length: 10001 samples (10.001 s)
我们再用 MNE-C 的默认设置滤波，它是一个采用频域技术设计的长时长、陡斜率 FIR 滤波器：
h = mne.filter.design_mne_c_filter(sfreq, l_freq=None, h_freq=f_p + 2.5)
x_mne_c = np.convolve(h, x)[len(h) // 2 :]
transition_band = 5 # Hz (default in MNE-C)
f_s = f_p + transition_band
freq = [0.0, f_p, f_s, sfreq / 2.0]
gain = [1.0, 1.0, 0.0, 0.0]
plot_filter(h, sfreq, freq, gain, "MNE-C default", compensate=True, **kwargs)
filter : 0.000 ... 42.5 Hz bins : 0 ... 348 of 4097 hpw : 3 lpw : 20
下面是一个最小相位滤波器的例子：
h = mne.filter.create_filter(
x,
sfreq,
l_freq=None,
h_freq=f_p,
phase="minimum",
fir_design="firwin",
verbose=True,
)
x_min = np.convolve(h, x)
transition_band = 0.25 * f_p
f_s = f_p + transition_band
filter_dur = 6.6 / transition_band # s
n = int(sfreq * filter_dur)
freq = [0.0, f_p, f_s, sfreq / 2.0]
gain = [1.0, 1.0, 0.0, 0.0]
plot_filter(h, sfreq, freq, gain, "Minimum-phase filter", **kwargs)
Setting up low-pass filter at 40 Hz
FIR filter parameters
---------------------
Designing a one-pass, non-linear phase, causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper transition bandwidth: 10.00 Hz
- Filter length: 331 samples (0.331 s)
MNE-Python 0.13 和 MNE-C 的滤波器都具有出色的频率衰减，但代价是时域中可能出现振铃（持续较久的纹波）。陡峭的滤波器会产生振铃，尤其是当信号的频率成分位于过渡带附近时。我们的 Morlet 小波信号在过渡带内也有能量，因此陡斜率、长时长的滤波器比更短、斜率更平缓的滤波器表现出更明显的时域振铃：
axes = plt.subplots(1, 2, layout="constrained")[1]
def plot_signal(x, offset):
"""Plot a signal."""
t = np.arange(len(x)) / sfreq
axes[0].plot(t, x + offset)
axes[0].set(xlabel="Time (s)", xlim=t[[0, -1]])
X = fft(x)
freqs = fftfreq(len(x), 1.0 / sfreq)
mask = freqs >= 0
X = X[mask]
freqs = freqs[mask]
axes[1].plot(freqs, 20 * np.log10(np.maximum(np.abs(X), 1e-16)))
axes[1].set(xlim=flim)
yscale = 30
yticklabels = [
"Original",
"Noisy",
"FIR-firwin (≥0.16)",
"FIR-firwin2 (0.14)",
"FIR-steep (≤0.13)",
"FIR-steep (MNE-C)",
"Minimum-phase",
]
yticks = -np.arange(len(yticklabels)) / yscale
plot_signal(x_orig, offset=yticks[0])
plot_signal(x, offset=yticks[1])
plot_signal(x_v16, offset=yticks[2])
plot_signal(x_v14, offset=yticks[3])
plot_signal(x_v13, offset=yticks[4])
plot_signal(x_mne_c, offset=yticks[5])
plot_signal(x_min, offset=yticks[6])
axes[0].set(
xlim=tlim,
title=f"FIR, Lowpass={f_p} Hz",
xticks=tticks,
ylim=[-len(yticks) / yscale, 1.0 / yscale],
yticks=yticks,
yticklabels=yticklabels,
)
for text in axes[0].get_yticklabels():
text.set(rotation=45, size=8)
axes[1].set(xlim=flim, ylim=(-60, 10), xlabel="Frequency (Hz)", ylabel="Magnitude (dB)")
plt.show()
IIR 滤波器#
MNE-Python 还提供基于 scipy.signal 中方法的 IIR 滤波功能。具体来说，我们使用通用函数 scipy.signal.iirfilter() 和 scipy.signal.iirdesign()，它们为 IIR 滤波器设计提供了统一的接口。
设计 IIR 滤波器#
我们继续设计 40 Hz 低通滤波器，并考察不同 IIR 滤波器的一些权衡。
通常默认的 IIR 滤波器是 Butterworth 滤波器，它的设计目标是具有最大平坦的通带。我们来看几个不同的滤波器阶数，即使用不同数量的系数、因而陡峭程度不同的滤波器：
注意
注意，下面 IIR 滤波器的群延迟（与相位有关）不是恒定的。对 FIR 而言，我们可以设计所谓的线性相位滤波器，其群延迟恒定，因而必要时可以补偿延迟（使滤波器变为非因果）。IIR 滤波器无法做到这一点，因为它们具有非线性相位（群延迟不恒定）。随着滤波器阶数增加，过渡带附近和过渡带内的相位畸变会加重。不过，如果可以使用非因果（前向-后向）滤波，例如借助 scipy.signal.filtfilt()，理论上可以缓解这些相位问题。
sos = signal.iirfilter(2, f_p / nyq, btype="low", ftype="butter", output="sos")
plot_filter(
dict(sos=sos), sfreq, freq, gain, "Butterworth order=2", compensate=True, **kwargs
)
x_shallow = signal.sosfiltfilt(sos, x)
del sos
这个滤波器的下降沿不是很陡。
注意
这里我们使用了二阶节（SOS）：通过 scipy.signal.sosfilt()，并且在向 scipy.signal.iirfilter() 传入 output='sos' 关键字参数时，底层调用了 scipy.signal.zpk2sos()。上面给出的滤波器定义使用多项式分子/分母（有时称为“tf”）形式 (b, a)，理论上与这里使用的 SOS 形式等价。然而在实际中，由于数值精度问题，SOS 形式能得到好得多的结果（示例见 scipy.signal.sosfilt()），因此应尽可能使用 SOS。
我们提高阶数，注意现在衰减更好，但脉冲响应也更长。同时改用 MNE 的滤波器设计函数，它简化了一些操作，并给出关于所得滤波器的一些信息：
iir_params = dict(order=8, ftype="butter")
filt = mne.filter.create_filter(
x, sfreq, l_freq=None, h_freq=f_p, method="iir", iir_params=iir_params, verbose=True
)
plot_filter(filt, sfreq, freq, gain, "Butterworth order=8", compensate=True, **kwargs)
x_steep = signal.sosfiltfilt(filt["sos"], x)
Setting up low-pass filter at 40 Hz
IIR filter parameters
---------------------
Butterworth lowpass zero-phase (two-pass forward and reverse) non-causal filter:
- Filter order 16 (effective, after forward-backward)
- Cutoff at 40.00 Hz: -6.02 dB
还有其他类型的 IIR 滤波器可以使用。完整列表请参阅 scipy.signal.iirdesign() 的文档。我们试试 Chebychev（I 型）滤波器，它以通带纹波为代价换取更好的阻带衰减：
iir_params.update(
ftype="cheby1",
rp=1.0, # dB of acceptable pass-band ripple
)
filt = mne.filter.create_filter(
x, sfreq, l_freq=None, h_freq=f_p, method="iir", iir_params=iir_params, verbose=True
)
plot_filter(
filt,
sfreq,
freq,
gain,
"Chebychev-1 order=8, ripple=1 dB",
compensate=True,
**kwargs,
)
Setting up low-pass filter at 40 Hz
IIR filter parameters
---------------------
Chebyshev I lowpass zero-phase (two-pass forward and reverse) non-causal filter:
- Filter order 16 (effective, after forward-backward)
- Cutoff at 40.00 Hz: -2.00 dB
如果能容忍更大的纹波，我们可以让它稍微更陡一些，但脉冲响应开始振铃明显更久（注意 x 轴尺度不同）：
iir_params["rp"] = 6.0
filt = mne.filter.create_filter(
x, sfreq, l_freq=None, h_freq=f_p, method="iir", iir_params=iir_params, verbose=True
)
plot_filter(
filt,
sfreq,
freq,
gain,
"Chebychev-1 order=8, ripple=6 dB",
compensate=True,
**kwargs,
)
Setting up low-pass filter at 40 Hz
IIR filter parameters
---------------------
Chebyshev I lowpass zero-phase (two-pass forward and reverse) non-causal filter:
- Filter order 16 (effective, after forward-backward)
- Cutoff at 40.00 Hz: -12.00 dB
与 FIR 滤波器类似，我们也可以定义因果 IIR 滤波器。
filt = mne.filter.create_filter(
x,
sfreq,
l_freq=None,
h_freq=f_p,
method="iir",
phase="forward",
iir_params=iir_params,
verbose=True,
)
plot_filter(
filt,
sfreq,
freq,
gain,
"Chebychev-1 order=8, ripple=6 dB",
compensate=False,
**kwargs,
)
Setting up low-pass filter at 40 Hz
IIR filter parameters
---------------------
Chebyshev I lowpass non-linear phase (one-pass forward) causal filter:
- Filter order 8 (forward)
- Cutoff at 40.00 Hz: -6.00 dB
应用 IIR 滤波器#
现在我们来看看之前那个 Morlet 信号在平缓和陡峭的 Butterworth IIR 滤波器下的表现：
axes = plt.subplots(1, 2, layout="constrained")[1]
yticks = np.arange(4) / -30.0
yticklabels = ["Original", "Noisy", "Butterworth-2", "Butterworth-8"]
plot_signal(x_orig, offset=yticks[0])
plot_signal(x, offset=yticks[1])
plot_signal(x_shallow, offset=yticks[2])
plot_signal(x_steep, offset=yticks[3])
axes[0].set(
xlim=tlim,
title=f"IIR, Lowpass={f_p} Hz",
xticks=tticks,
ylim=[-0.125, 0.025],
yticks=yticks,
yticklabels=yticklabels,
)
for text in axes[0].get_yticklabels():
text.set(rotation=45, size=8)
axes[1].set(xlim=flim, ylim=(-60, 10), xlabel="Frequency (Hz)", ylabel="Magnitude (dB)")
plt.show()
滤波的一些陷阱#
近期有多篇论文指出，误用滤波器可能导致错误推断的潜在风险。
低通滤波的问题#
一般来说，滤波器（尤其是非因果的零相位滤波器）可能使活动看起来比实际发生得更早或更晚。正如 VanRullen (2011) [ 5 ] 所述，对当时常用的低通滤波器的研究在应用于仿真数据时产生了伪迹。不过，在 Rousselet (2012) [ 6 ] 的许多真实数据示例中，这种有害影响很小。
也许更能说明问题的是，Widmann & Schröger (2012) [ 4 ] 指出，VanRullen (2011) [ 5 ] 中有问题的低通滤波器：
使用了最小二乘设计（类似
scipy.signal.firls()），其中包含“无需关心”的过渡区域，这可能导致不受控制的行为。滤波器长度与过渡带宽无关，这会造成过度振铃和信号畸变。
高通滤波的问题#
在高通滤波方面，Acunzo et al. (2012) [ 7 ] 发现，使用高于 0.1 Hz 的转角频率会：
“……产生系统性偏差，很容易导致对神经活动的误读。”
在一篇相关论文中，Widmann et al. (2015) [ 3 ] 也建议使用 0.1 Hz 的高通。Tanner et al. (2015) [ 8 ] 随后提供了更多此类畸变的证据。使用语义和句法加工的语言 ERP 研究数据（即 N400 和 P600），与未滤波数据相比，使用高于 0.3 Hz 的高通会导致显著效应出现过早，早得不合情理。据此，作者建议语言加工的最优高通值为 0.1 Hz。
我们可以重现 Tanner et al. (2015) [ 8 ] 中一个有问题的仿真：
“仿真成分是一个单周期余弦波，幅值为 5µV [原文如此]，在刺激后 500 ms 起始，持续 800 ms。为避免滤波边缘效应，仿真成分被嵌入 20 s 的零值之中。[…] 畸变 [是] 由 2 Hz 低通和高通滤波器造成的。[…] 使用 30 Hz 低通和 0.01 Hz 高通滤波器时，原始波形 [没有出现] 可见畸变。[…] 滤波频率对应于半幅（-6 dB）截止（12 dB/倍频程滚降）。”
注意
这个仿真信号不仅在通带内含有能量，在过渡带和阻带内也含有能量——最容易的理解方式或许是：该信号具有非零的直流值；同时它是一个经过加窗（这里乘以矩形窗）的平移余弦波，这会使余弦频率和直流频率扩散到其他频率（时域相乘对应频域卷积，因此在时域乘以矩形窗，意味着在频域将一个 sinc 函数与直流处和余弦频率处的冲激进行卷积）。
x = np.zeros(int(2 * sfreq))
t = np.arange(0, len(x)) / sfreq - 0.2
onset = np.where(t >= 0.5)[0][0]
cos_t = np.arange(0, int(sfreq * 0.8)) / sfreq
sig = 2.5 - 2.5 * np.cos(2 * np.pi * (1.0 / 0.8) * cos_t)
x[onset : onset + len(sig)] = sig
iir_lp_30 = signal.iirfilter(2, 30.0 / sfreq, btype="lowpass")
iir_hp_p1 = signal.iirfilter(2, 0.1 / sfreq, btype="highpass")
iir_lp_2 = signal.iirfilter(2, 2.0 / sfreq, btype="lowpass")
iir_hp_2 = signal.iirfilter(2, 2.0 / sfreq, btype="highpass")
x_lp_30 = signal.filtfilt(iir_lp_30[0], iir_lp_30[1], x, padlen=0)
x_hp_p1 = signal.filtfilt(iir_hp_p1[0], iir_hp_p1[1], x, padlen=0)
x_lp_2 = signal.filtfilt(iir_lp_2[0], iir_lp_2[1], x, padlen=0)
x_hp_2 = signal.filtfilt(iir_hp_2[0], iir_hp_2[1], x, padlen=0)
xlim = t[[0, -1]]
ylim = [-2, 6]
xlabel = "Time (s)"
ylabel = r"Amplitude ($\mu$V)"
tticks = [0, 0.5, 1.3, t[-1]]
axes = plt.subplots(2, 2, layout="constrained")[1].ravel()
for ax, x_f, title in zip(
axes,
[x_lp_2, x_lp_30, x_hp_2, x_hp_p1],
["LP$_2$", "LP$_{30}$", "HP$_2$", "LP$_{0.1}$"],
):
ax.plot(t, x, color="0.5")
ax.plot(t, x_f, color="k", linestyle="--")
ax.set(
ylim=ylim, xlim=xlim, xticks=tticks, title=title, xlabel=xlabel, ylabel=ylabel
)
plt.show()
类似地，在 Kappenman & Luck (2010) [ 9 ] 报告的 P300 范式中，他们发现施加 1 Hz 高通会降低在 N100 响应中发现显著差异的概率，原因可能是 P300 响应被高通滤波器在时间上抹开（并发生反转），从而倾向于抵消增强的 N100。不过他们仍指出，一定程度的高通滤波对于处理数据漂移仍然有用。
尽管这些论文通常建议高通使用 0.1 Hz 或更低的频率，但重要的是要记住（正如大多数作者所指出），滤波选择应取决于感兴趣信号的频率成分以及待抑制噪声的频率成分。例如，在一些使用 Sample 数据集的 MNE-Python 示例中，观察听觉或视觉 N100 响应时会使用约 1 Hz 的高通值，因为我们分析的是标准（而非偏差）试次，因而预期来自较晚或较慢成分的污染有限。
基线问题（还是解决方案？）#
在这场持续发展的讨论中，Tanner et al. (2015) [ 8 ] 建议使用基线校正来去除数据中的慢漂移。然而，Maess et al. (2016) [ 10 ] 认为，基线校正作为一种高通滤波形式，相比标准高通滤波并无显著优势。Tanner et al. (2016) [ 11 ] 反驳说，基线校正可以修正滤波带来的问题。
为了解他们的意思，我们再来看之前那个旧的仿真信号 x：
def baseline_plot(x):
fig, all_axes = plt.subplots(3, 2, layout="constrained")
for ri, (axes, freq) in enumerate(zip(all_axes, [0.1, 0.3, 0.5])):
for ci, ax in enumerate(axes):
if ci == 0:
iir_hp = signal.iirfilter(
4, freq / sfreq, btype="highpass", output="sos"
)
x_hp = signal.sosfiltfilt(iir_hp, x, padlen=0)
else:
x_hp -= x_hp[t < 0].mean()
ax.plot(t, x, color="0.5")
ax.plot(t, x_hp, color="k", linestyle="--")
if ri == 0:
ax.set(title=("No " if ci == 0 else "") + "Baseline Correction")
ax.set(xticks=tticks, ylim=ylim, xlim=xlim, xlabel=xlabel)
ax.set_ylabel(f"{freq:0.1f} Hz", rotation=0, horizontalalignment="right")
fig.suptitle(title)
plt.show()
baseline_plot(x)
作为回应，Maess et al. (2016) [ 12 ] 指出，这些仿真并未涉及各条件共有的刺激前活动，因为应用基线校正实际上会把基线期的拓扑复制到基线时段之外。如果我们给信号 x 加上一些一致的刺激前活动，就能看到这一点：一切都会变得很糟。
注意
看这些图时需要注意的重要一点是，它们针对的是单个仿真传感器。在多电极记录中，刺激前活动的拓扑（即空间模式）会泄漏到刺激后时段。由于平均后的刺激前空间模式会从各传感器的时间过程中被减去，这很可能造成时域信号的空间变化畸变。
在基线期加入一些活动：
两方似乎都承认，滤波截止频率的选择、乃至是否应用基线校正，都取决于所研究数据的特性，尤其是在以下方面：
潜在诱发活动的频率成分相对于滤波参数的关系。
基线期内不存在一致诱发活动这一假设的有效性。
因此，我们建议根据待分析数据的特性，谨慎地应用基线校正和/或高通值。
滤波默认设置#
MNE-Python 中的默认设置#
在 MNE-Python 中，滤波大多在 mne.io.Raw 层级进行，因此使用 mne.io.Raw.filter()。该函数在底层（除其他操作外）会调用 mne.filter.filter_data() 来实际对数据滤波，后者默认应用由 scipy.signal.firwin() 设计的零相位 FIR 滤波器。在 Widmann et al. (2015) [ 3 ] 中，他们建议了一组用于高通滤波的具体参数，包括：
“……提供下限通带边缘 25% 的过渡带宽，但在可能的情况下不低于 2 Hz，否则就取通带边缘到临界频率的距离。”
在实践中，这意味着如果采样率为 100 Hz（即 Nyquist 频率为 50 Hz），下表中每个高通值 l_freq 或低通值 h_freq 都会分别对应相应的 l_trans_bandwidth 或 h_trans_bandwidth：
l_freq 或 h_freq |
l_trans_bandwidth |
h_trans_bandwidth |
|---|---|---|
0.01 |
0.01 |
2.0 |
0.1 |
0.1 |
2.0 |
1.0 |
1.0 |
2.0 |
2.0 |
2.0 |
2.0 |
4.0 |
2.0 |
2.0 |
8.0 |
2.0 |
2.0 |
10.0 |
2.5 |
2.5 |
20.0 |
5.0 |
5.0 |
40.0 |
10.0 |
10.0 |
50.0 |
12.5 |
12.5 |
在使用 l_trans_bandwidth='auto' 和 h_trans_bandwidth='auto' 时，MNE-Python 已采用这一定义来选择高通（以及低通）的过渡带宽。
使用 filter_length='auto' 自动选择滤波器长度时，会使用最短过渡带宽的倒数，以确保在阻带频率处有足够的衰减。具体来说，该倒数（以采样点为单位）对于 Hann、Hamming 或 Blackman 窗分别乘以 3.1、3.3 或 5.0（窗由 fir_design='firwin' 时的 fir_window 参数选择），而在 fir_design='firwin2' 模式下则将这些数值翻倍。
注意
对于 fir_design='firwin2'，乘数比 Ifeachor & Jervis (2002) [ 2 ]（第 357 页）中给出的值翻了一倍，因为 scipy.signal.firwin2() 会对频率响应产生涂抹效应，我们通过增加滤波器长度来补偿。这就是 fir_design='firwin' 优于 fir_design='firwin2' 的原因。
在 0.14 版本中，我们在滤波器设计中默认使用 Hamming 窗，因为它可提供高达 53 dB 的阻带衰减，且通带纹波很小。
注意
在带通应用中，低通滤波器通常可以用比高通滤波器更少的采样点有效工作，因此在使用 fir_design='firwin2' 时，建议分别应用高通和低通。对于 fir_design='firwin' 设计模式，则无需分开操作，因为低通和高通部分会分别构造以满足过渡带要求。
关于如何在真实数据上使用 MNE-Python 的滤波函数，请参阅关于滤波与重采样的预处理教程。
MNE-C 中的默认设置#
MNE-C 默认使用：
低通滤波器的 5 Hz 过渡带。
高通滤波器的 3 个采样点过渡带。
8197 个采样点的滤波器长度。
滤波器在频域中设计，构造为线性相位滤波器，从而像 MNE-Python 的 phase='zero' 滤波选项那样对延迟进行补偿。
过渡区使用平方余弦斜坡。由于用它代替更平缓的（例如线性）过渡，在给定过渡宽度下，相比窗函数 FIR 设计中相同过渡宽度，会产生更多的时域振铃，但衰减也更迅速。
对于 M/EEG 数据中常见的采样率（例如 500-2000 Hz），默认滤波器长度通常具有出色的衰减，但振铃较长。
其他软件中的默认设置#
Widmann et al. (2015) [ 3 ] 中对各种软件包中的滤波做了很好的比较（可能有些过时）。简要来说：
- EEGLAB
- MNE-Python 0.14 的默认行为与 EEGLAB 非常相似（更多信息见 EEGLAB 滤波常见问题）。
- FieldTrip
- 默认情况下，FieldTrip 应用前向-后向 Butterworth IIR 滤波器，阶数为 4（带通和带阻滤波器）或 2（低通和高通滤波器）。在 MNE-Python 中使用
raw.filter(..., method='iir')滤波可以实现类似的滤波器（选项另见mne.filter.construct_iir_filter()）。更多信息请参阅 FieldTrip 带通文档。
报告滤波器参数#
在 Widmann et al. (2015) [ 3 ] 第 45 页，有一份便于使用的列表，列出了每次发表都应报告的重要滤波参数：
滤波器类型（高通、低通、带通、带阻、FIR、IIR）
截止频率（包括其定义）
滤波器阶数（或长度）
滚降或过渡带宽
通带纹波与阻带衰减
滤波器延迟（零相位、线性相位、非线性相位）与因果性
计算方向（单次前向/反向，或前向与反向两次）
下面我们将说明在 MNE 中如何处理这些参数。
滤波器类型#
根据所使用的函数或方法，可以指定滤波器类型。举例来说，在 mne.filter.create_filter() 中，相关参数为 l_freq、h_freq、method，如果方法为 FIR，则还有 fir_window 和 fir_design。
截止频率#
MNE 中 FIR 滤波器的截止频率定义为过渡带中点的半幅截止频率。也就是说，如果你构造一个 h_freq = 40 的低通 FIR 滤波器，滤波函数会给出一个取决于 h_trans_bandwidth 参数的过渡带宽。低通 FIR 滤波器期望的半幅截止频率则位于 h_freq + transition_bandwidth/2。
滤波器长度（阶数）与过渡带宽（滚降）#
在 MNE-Python 中的默认设置一节中，我们已经讨论了在未通过相应滤波函数参数指定自定义值时所使用的默认滤波器长度和过渡带宽。
如果你想了解通过“auto”设置所使用的滤波器长度和过渡带宽，可以用 mne.filter.create_filter() 再次打印这些设置：
# Use the same settings as when calling e.g., `raw.filter()`
fir_coefs = mne.filter.create_filter(
data=None, # data is only used for sanity checking, not strictly needed
sfreq=1000.0, # sfreq of your data in Hz
l_freq=None,
h_freq=40.0, # assuming a lowpass of 40 Hz
method="fir",
fir_window="hamming",
fir_design="firwin",
verbose=True,
)
# See the printed log for the transition bandwidth and filter length.
# Alternatively, get the filter length through:
filter_length = fir_coefs.shape[0]
No data specified. Sanity checks related to the length of the signal relative to the filter order will be skipped.
Setting up low-pass filter at 40 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal lowpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 331 samples (0.331 s)
注意
如果你使用的是 IIR 滤波器，mne.filter.create_filter() 不会在日志中打印滤波器长度和过渡带宽。你可以改用 iir_params 参数指定滚降，或保留默认值，即四阶（Butterworth）滤波器。
通带纹波与阻带衰减#
使用标准的 scipy.signal.firwin() 设计时（MNE 中的 FIR 滤波器即如此），通带纹波和阻带衰减取决于设计所用的窗。标准窗的数值列于下表（见 Ifeachor & Jervis (2002) [ 2 ]，第 357 页）：
| 窗函数名称 | 通带纹波 | 阻带衰减 |
|---|---|---|
Hann |
0.0545 dB |
44 dB |
Hamming |
0.0194 dB |
53 dB |
Blackman |
0.0017 dB |
74 dB |
滤波器延迟与计算方向#
要报告这些信息，阅读所应用滤波函数或方法的文档字符串也许就足够了。例如，在 mne.filter.create_filter 的文档字符串中，关于 phase 参数是这样说的：
滤波器的相位，仅在
method='fir'时使用。默认构造对称的线性相位 FIR 滤波器。如果phase='zero'（默认），则对该滤波器的延迟进行补偿。如果phase=='zero-double'，则该滤波器会被应用两次，一次前向、一次反向。如果为'minimum'，则使用最小相位因果滤波器。
小结#
滤波时总有一些权衡需要考虑。一个重要权衡是在时域特性（如振铃）与频域衰减特性（如有效过渡带宽）之间取得平衡。频率截止陡峭的滤波器作用于过渡带内含频率成分的信号时，其输出可能长时间振铃。因此一般来说，能容忍的过渡带越宽，滤波器在时域中的表现就越好。
参考文献#
脚本运行总耗时：（0 分钟 8.039 秒）