如何使用 Scipy.signal.butter 实现带通巴特沃斯滤波器

Posted 2023-02-14

【中文标题】如何使用 Scipy.signal.butter 实现带通巴特沃斯滤波器

【英文标题】：How to implement band-pass Butterworth filter with Scipy.signal.butter

【发布时间】：2012-08-19 02:02:17

【问题描述】：

更新：

我根据这个问题找到了一个 Scipy 食谱！所以，有兴趣的可以直接去：Contents » Signal processing » Butterworth Bandpass

---

我很难实现最初看似简单的任务，即为一维 numpy 数组（时间序列）实现巴特沃斯带通滤波器。

我必须包含的参数是采样率、以赫兹为单位的截止频率和可能的顺序（其他参数，如衰减、固有频率等对我来说比较模糊，所以任何“默认”值都可以）。

我现在拥有的是这个，它似乎可以用作高通滤波​​器，但我不确定我是否做得对：

def butter_highpass(interval, sampling_rate, cutoff, order=5):
    nyq = sampling_rate * 0.5

    stopfreq = float(cutoff)
    cornerfreq = 0.4 * stopfreq  # (?)

    ws = cornerfreq/nyq
    wp = stopfreq/nyq

    # for bandpass:
    # wp = [0.2, 0.5], ws = [0.1, 0.6]

    N, wn = scipy.signal.buttord(wp, ws, 3, 16)   # (?)

    # for hardcoded order:
    # N = order

    b, a = scipy.signal.butter(N, wn, btype='high')   # should 'high' be here for bandpass?
    sf = scipy.signal.lfilter(b, a, interval)
    return sf

文档和示例令人困惑且晦涩难懂，但我想实现标为“for bandpass”的推荐中呈现的表单。 cmets 中的问号表示我只是复制粘贴了一些示例而没有理解发生了什么。

我不是电气工程师或科学家，只是需要对 EMG 信号执行一些相当简单的带通滤波的医疗设备设计师。

【问题讨论】：

我在 dsp.stackexchange 上尝试过一些东西，但他们过于关注工程的概念问题（我无法处理），而不是使用 scipy 函数。

【参考方案1】：

对于带通滤波器，ws 是一个包含上下角频率的元组。这些代表滤波器响应比通带小 3 dB 的数字频率。

wp 是一个包含阻带数字频率的元组。它们代表最大衰减开始的位置。

gpass 是通带中的最大衰减，以 dB 为单位，而 gstop 是阻带中的衰减。

例如，假设您想设计一个采样率为 8000 个样本/秒的滤波器，其转角频率为 300 和 3100 Hz。奈奎斯特频率是采样率除以二，在本例中为 4000 Hz。等效数字频率为 1.0。那么这两个拐角频率是 300/4000 和 3100/4000。

现在假设您希望阻带从拐角频率降低 30 dB +/- 100 Hz。因此，您的阻带将从 200 和 3200 Hz 开始，从而产生 200/4000 和 3200/4000 的数字频率。

要创建您的过滤器，您可以将 buttord 称为

fs = 8000.0
fso2 = fs/2
N,wn = scipy.signal.buttord(ws=[300/fso2,3100/fso2], wp=[200/fs02,3200/fs02],
   gpass=0.0, gstop=30.0)

生成的滤波器的长度将取决于阻带的深度和响应曲线的陡度，响应曲线由拐角频率和阻带频率之间的差异决定。

【讨论】：

我尝试实现它，但仍然缺少一些东西。一件事是gpass=0.0 引发了除以零错误，所以我将其更改为 0.1 并且错误停止了。除此之外，butter 的文档说：Passband and stopband edge frequencies, normalized from 0 to 1 (1 corresponds to pi radians / sample). 我怀疑你的答案是否计算正确，所以我仍在努力，很快就会给出一些反馈。

（另外，虽然我的ws 和wp 各有两个元素，但滤波器只执行低通或高通（通过btype 参数），而不是带通）

根据docs.scipy.org/doc/scipy/reference/generated/… 的文档，buttord 设计了低通、高通和带通滤波器。至于 gpass，我猜 buttord 不允许在通带中衰减 0 dB。然后将其设置为某个非零值。

【参考方案2】：

您可以跳过使用 buttord，而只需为过滤器选择一个顺序，看看它是否符合您的过滤条件。要生成带通滤波器的滤波器系数，请给 butter() 滤波器阶数、截止频率Wn=[low, high]（表示为奈奎斯特频率的分数，即采样频率的一半）和波段类型btype="band"。

这是一个脚本，它定义了几个使用 Butterworth 带通滤波器的便利函数。当作为脚本运行时，它会生成两个图。一张显示了相同采样率和截止频率下几个滤波器阶数的频率响应。另一个图展示了过滤器（order=6）对样本时间序列的影响。

from scipy.signal import butter, lfilter


def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y


if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import freqz

    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 5000.0
    lowcut = 500.0
    highcut = 1250.0

    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [3, 6, 9]:
        b, a = butter_bandpass(lowcut, highcut, fs, order=order)
        w, h = freqz(b, a, worN=2000)
        plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)

    plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)],
             '--', label='sqrt(0.5)')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    # Filter a noisy signal.
    T = 0.05
    nsamples = T * fs
    t = np.linspace(0, T, nsamples, endpoint=False)
    a = 0.02
    f0 = 600.0
    x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
    x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
    x += a * np.cos(2 * np.pi * f0 * t + .11)
    x += 0.03 * np.cos(2 * np.pi * 2000 * t)
    plt.figure(2)
    plt.clf()
    plt.plot(t, x, label='Noisy signal')

    y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)
    plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
    plt.xlabel('time (seconds)')
    plt.hlines([-a, a], 0, T, linestyles='--')
    plt.grid(True)
    plt.axis('tight')
    plt.legend(loc='upper left')

    plt.show()

以下是此脚本生成的图：

【讨论】：

你知道为什么过滤后的输出总是从零开始吗？是否可以与实际输入值x[0]匹配？我用 Cheby1 低通滤波器尝试了类似的东西，但我遇到了同样的问题。

@LWZ：使用函数scipy.signal.lfilter_zi，以及zi 参数到lfilter。有关详细信息，请参阅lfilter_zi 的文档字符串。 TL；博士？只需将y = lfilter(b, a, data) 更改为zi = lfilter_zi(b, a); y, zo = lfilter(b, a, data, zi=zi*data[0])。 （但这可能对带通或高通滤波器没有影响。）

我注意到scipy.signal.lfiter()的输出相对于原始信号和signal.filtfilt()的输出有180度的相移，这是为什么呢？如果时间对我很重要，我应该改用 filtfilt() 吗？

这是该频率下滤波器的相位延迟。通过巴特沃斯滤波器的正弦波的相位延迟非线性地取决于频率。对于零相位延迟，是的，您可以使用filtfilt()。我的回答 here 包括一个使用 filtfilt() 来避免过滤器引起的延迟的示例。

嗨 Jason，我建议您在 dsp.stackexchange.com 上提问有关信号处理理论的问题。如果您对自己编写的某些代码没有按预期工作有疑问，可以在 *** 上提出一个新问题。

【参考方案3】：

accepted answer中的过滤器设计方法是正确的，但是有一个缺陷。用 b 和 a 设计的 SciPy 带通滤波器是 unstable，可能会在 higher filter orders 产生 erroneous filters。

改为使用滤波器设计的 sos（二阶部分）输出。

from scipy.signal import butter, sosfilt, sosfreqz

def butter_bandpass(lowcut, highcut, fs, order=5):
        nyq = 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq
        sos = butter(order, [low, high], analog=False, btype='band', output='sos')
        return sos

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        y = sosfilt(sos, data)
        return y

此外，您可以通过更改来绘制频率响应

b, a = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = freqz(b, a, worN=2000)

到

sos = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = sosfreqz(sos, worN=2000)

【讨论】：

+1 因为现在在许多情况下这是更好的方法。与已接受答案的 cmets 一样，也可以使用前向后向滤波来消除相位延迟。只需将sosfilt 替换为sosfiltfilt。

@Mike 和 user13107 相同的错误是否也会影响高通和低通巴特沃斯滤波器？解决方案是一样的吗？

@dewarrn1 将其称为“错误”实际上是不正确的；该算法是正确实现的，但它本质上是不稳定的，所以它只是一个糟糕的算法选择。但是，是的，它会影响 any 高阶滤波器——不仅仅是高通或低通滤波器，也不仅仅是巴特沃斯滤波器，还有其他 like Chebyshev 等等。无论如何，一般来说，最好总是选择sos 输出，因为这样总能避免不稳定。除非您需要实时处理，否则您应该始终使用sosfiltfilt。

对不起，我很久以前没有注意到这个答案！ @ user13107，是的，当滤波器的阶数很大时，线性滤波器的传递函数（或“ba”）表示会出现一些严重的数值问题。当所需带宽与采样频率相比较小时，即使是相对低阶的滤波器也会出现问题。我的原始答案是在将 SOS 表示添加到 SciPy 之前编写的，并且在将 fs 参数添加到 scipy.signal 中的许多函数之前编写。答案早就应该更新了。

对此有什么帮助吗？ ***.com/q/60866193/5025009