如何使用Scipy.signal.butter实现带通Butterworth滤波器

Question

如何使用Scipy.signal.butter实现带通Butterworth滤波器

pythonscipysignal-processingdigital-filter

108

更新:

我在stackoverflow上找到了一篇基于此问题的Scipy食谱！所以，对于任何有兴趣的人，请直接访问：目录 » 信号处理 » 巴特沃斯带通滤波器

我现在正在尝试实现一个针对1-D 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

enter image description here

虽然文档和示例令人困惑和晦涩，但我想要实现在标记为“for bandpass”的评论中呈现的表单。注释中的问号表示我只是复制粘贴了一些例子，而并不理解正在发生的事情。

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

- heltonbiker

我在dsp.stackexchange尝试了一些东西，但他们过于关注（超出我的能力范围）工程的概念问题，而不是如何使用scipy函数。 - heltonbiker

3个回答

63

已被接受的答案中的滤波器设计方法是正确的，但它有一个缺陷。使用 b、a 设计的 SciPy 带通滤波器是不稳定的，可能会导致更高的滤波器阶数的错误滤波器。

相反，使用滤波器设计的 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

此外，您可以通过更改{{参数}}来绘制频率响应。该参数涉及{{IT技术相关内容}}。请注意，本文不会进行解释。

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

to

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

- user13107

8

在许多情况下，现在采用这种方法更好，因此我会给它点赞。正如对被接受答案的评论所述，还可以使用前后滤波器消除相位延迟。只需将sosfilt替换为sosfiltfilt即可。 - Mike

@Mike 和 user13107，高通和低通 Butterworth 滤波器是否也受到同样的 bug 影响？解决方案是否相同？ - dewarrn1

9

实际上称其为“bug”是不正确的；该算法已经被正确实现，但它本质上是不稳定的，因此选择这种算法只是一个糟糕的选择。但是，是的，它会影响到任何高阶滤波器-不仅仅是高通或低通滤波器，也包括其他滤波器，如切比雪夫滤波器等。总之，一般最好始终选择sos输出，因为这样可以始终避免不稳定性。除非您需要实时处理，否则应始终使用sosfiltfilt。 - Mike

1

抱歉我之前没有注意到这个答案！@user13107，是的，线性滤波器的传递函数（或“ba”）表示在滤波器阶数较大时存在一些严重的数值问题。即使是相对较低阶的滤波器，在所需带宽与采样频率相比较小的情况下也可能存在问题。我的原始答案是在SciPy添加SOS表示和在scipy.signal中的许多函数中添加fs参数之前编写的。这个答案已经过期需要更新了。 - Warren Weckesser

@Mike 但是你怎么知道未来数据的长度呢？假设我们有两个多正弦信号（相同频率分量），持续时间不同，我们使用 *filtfilt 进行离线滤波，即使我们使用相同的滤波器，由于持续时间不同，我们将得到不同的滤波数据？ - John

显示剩余6条评论

4

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

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

gpass是通带中的最大衰减（以分贝为单位），而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)

生成的滤波器长度取决于阻带深度和响应曲线的陡峭程度，后者由角频率和阻带频率之差决定。

- sizzzzlerz

我尝试实现它，但仍有一些缺失。其中一个问题是gpass=0.0会引发除以零的错误，所以我将其更改为0.1，错误停止了。此外，butter的文档说：“通带和阻带边缘频率从0到1进行归一化（1对应于每个样本的pi弧度）。”我不确定您的答案是否正确，因此我仍在努力，并将很快给出反馈。 - heltonbiker

虽然我的 ws 和 wp 每个都有两个元素，但滤波器仅执行低通或高通（通过 btype 参数），而不是带通。 - heltonbiker

1

根据http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.buttord.html上的文档，buttord设计低通、高通和带通滤波器。就gpass而言，我猜buttord不允许通过带中有0 dB的衰减。将其设置为非零值即可。 - sizzzzlerz

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Warren Weckesser · Accepted Answer

您可以跳过使用buttord函数，而是选择一个滤波器的顺序并查看它是否符合您的过滤条件。要为带通滤波器生成滤波器系数，请将filter order、截止频率Wn=[lowcut, highcut]、采样率fs（以与截止频率相同的单位表示）和波段类型btype="band"传递给butter()。

下面是一个脚本，定义了几个方便的函数来处理Butterworth带通滤波器。当作为脚本运行时，它会生成两个图表。其中一个显示了在相同的采样率和截止频率下，不同滤波器阶数的频率响应。另一个演示了滤波器（order=6）对样本时间序列的影响。

from scipy.signal import butter, lfilter

def butter_bandpass(lowcut, highcut, fs, order=5):
    return butter(order, [lowcut, highcut], fs=fs, btype='band')

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, fs=fs, worN=2000)
        plt.plot(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.arange(0, nsamples) / fs
    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()

这个脚本生成的图表如下：多个滤波器阶数的频率响应