如何使用 Scipy.signal.butter 实现带通巴特沃斯滤波器 [英] How to implement band-pass Butterworth filter with Scipy.signal.butter

问题描述

更新:

我在这个问题中找到了一个 Scipy Recipe！因此，对于任何感兴趣的人，请直接访问:内容 » 信号处理 » Butterworth Bandpass

<小时>

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

我必须包含的参数是采样率、以赫兹为单位的截止频率以及可能的顺序(其他参数，如衰减、自然频率等对我来说更模糊，因此任何默认"值都可以).

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

def butter_highpass(interval, sampling_rate, cutoff, order=5):nyq = 采样率 * 0.5stopfreq = 浮动(截止)角频率 = 0.4 * 停止频率 # (?)ws =cornerfreq/nyqwp = stopfreq/nyq# 带通:# wp = [0.2, 0.5], ws = [0.1, 0.6]N, wn = scipy.signal.buttord(wp, ws, 3, 16) # (?)# 对于硬编码顺序:# N = 顺序b, a = scipy.signal.butter(N, wn, btype='high') # 'high' 应该在这里用于带通吗?sf = scipy.signal.lfilter(b, a, 间隔)返回SF

文档和示例令人困惑和晦涩，但我想实现标记为for bandpass"的推荐中显示的表单.评论中的问号表明我只是复制粘贴了一些示例而没有理解发生了什么.

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

解决方案

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

这是一个脚本，它定义了几个使用巴特沃斯带通滤波器的便利函数.当作为脚本运行时，它会绘制两个图.一个显示了相同采样率和截止频率下几个滤波器阶数的频率响应.另一幅图展示了滤波器(阶数=6)对样本时间序列的影响.

from scipy.signal 导入黄油，lfilterdef butter_bandpass(lowcut, highcut, fs, order=5):nyq = 0.5 * fs低 = 低切/nyq高 = 高切/nyqb, a = butter(order, [low, high], btype='band')返回 b, adef butter_bandpass_filter(data, lowcut, highcut, fs, order=5):b, a = butter_bandpass(lowcut, highcut, fs, order=order)y = lfilter(b, a, 数据)返回 y如果 __name__ == "__main__":将 numpy 导入为 np导入 matplotlib.pyplot 作为 plt从 scipy.signal 导入 freqz# 采样率和所需的截止频率(以赫兹为单位).fs = 5000.0低切 = 500.0高切 = 1250.0# 绘制几个不同阶数的频率响应.plt.figure(1)plt.clf()对于 [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('频率(Hz)')plt.ylabel('增益')plt.grid(真)plt.legend(loc='best')# 过滤噪声信号.T = 0.05nsamples = T * fst = np.linspace(0, T, nsamples, endpoint=False)a = 0.02f0 = 600.0x = 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='噪声信号')y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)plt.plot(t, y, label='过滤后的信号 (%g Hz)' % f0)plt.xlabel('时间(秒)')plt.hlines([-a, a], 0, T, linestyles='--')plt.grid(真)plt.axis('紧')plt.legend(loc='左上')plt.show()

以下是此脚本生成的图:

UPDATE:

I found a Scipy Recipe based in this question! So, for anyone interested, go straight to: Contents » Signal processing » Butterworth Bandpass

---

I'm having a hard time to achieve what seemed initially a simple task of implementing a Butterworth band-pass filter for 1-D numpy array (time-series).

The parameters I have to include are the sample_rate, cutoff frequencies IN HERTZ and possibly order (other parameters, like attenuation, natural frequency, etc. are more obscure to me, so any "default" value would do).

What I have now is this, which seems to work as a high-pass filter but I'm no way sure if I'm doing it right:

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

The docs and examples are confusing and obscure, but I'd like to implement the form presented in the commend marked as "for bandpass". The question marks in the comments show where I just copy-pasted some example without understanding what is happening.

I am no electrical engineering or scientist, just a medical equipment designer needing to perform some rather straightforward bandpass filtering on EMG signals.

解决方案

You could skip the use of buttord, and instead just pick an order for the filter and see if it meets your filtering criterion. To generate the filter coefficients for a bandpass filter, give butter() the filter order, the cutoff frequencies Wn=[low, high] (expressed as the fraction of the Nyquist frequency, which is half the sampling frequency) and the band type btype="band".

Here's a script that defines a couple convenience functions for working with a Butterworth bandpass filter. When run as a script, it makes two plots. One shows the frequency response at several filter orders for the same sampling rate and cutoff frequencies. The other plot demonstrates the effect of the filter (with order=6) on a sample time series.

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()

Here are the plots that are generated by this script:

这篇关于如何使用 Scipy.signal.butter 实现带通巴特沃斯滤波器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！