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

Question

更新:

我找到了一个基于这个问题的Scipy Recipe!所以,对于任何有兴趣的人,请直接:

http://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html

---

我很难实现最初为一维numpy阵列(时间序列)实现Butterworth带通滤波器的简单任务.

我必须包括的参数是sample_rate,截止频率IN HERTZ和可能的顺序(其他参数,如衰减,固有频率等对我来说更加模糊,因此任何"默认"值都可以).

我现在拥有的是这个,它似乎是一个高通滤波器,但我不确定我是否做得对:

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"的表述中的表格.评论中的问号显示我只是在不了解正在发生的事情的情况下复制粘贴了一些示例.

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

谢谢你的帮助!

Answer 1

您可以跳过使用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()

以下是此脚本生成的图:

Answer 2

接受答案中的过滤器设计方法是正确的,但它有一个缺陷.使用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

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

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)

Answer 3

对于带通滤波器，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)

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