如何实现带通与Scipy.signal.butter巴特沃斯滤波器(How to implement

更新：

我发现，总部设在这个问题上SciPy的食谱！ 因此，对于任何有兴趣，直奔： 目录»信号处理»巴特沃思带通

---

我有一个很难实现的目标最初似乎实现巴特沃斯带通滤波器1-d numpy的阵列（时间序列）的一个简单的任务。

我必须包括的参数是SAMPLE_RATE，以赫兹并可能阶截止频率（其它参数，如衰减，固有频率等都是比较模糊的我，所以任何“默认”值会做）。

我现在是这样的，这似乎为高通滤波器的工作，但我没有办法知道我在做正确的事：

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

该文档和例子是混乱和模糊的，但我想实现在标记为“为带通”的表彰形式呈现。 在评论的问号表明，我只是复制粘贴一些示例不理解发生了什么。

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

Answer 1:

你可以跳过使用buttord的，而是只挑选过滤器的订单，看看它是否符合您的过滤标准。 为了产生用于一个带通滤波器的滤波器系数，得到黄油（）滤波器的阶数，截止频率Wn=[low, high] （表示为Nyquist频率，它是采样频率的一半的部分）和带类型btype="band" 。

下面是定义与巴特沃斯带通滤波器工作的一对夫妇的便利功能的脚本。 当作为一个脚本运行时，它两个地块。 一个示出了在几个滤波器阶对于相同采样率和截止频率的频率响应。 其他曲线表明了过滤器的上一个样本的时间序列的影响（与顺序= 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个分贝。

可湿性粉剂是含阻带数字频率的元组。 他们表示，其中最大衰减开始的位置。

GPASS是在以dB通带的最大attenutation而gstop是在阻带的attentuation。

说，例如，你想设计一个过滤器用于8000个采样/秒具有300和3100赫兹的角频率的采样率。 奈奎斯特频率是由两个分割的采样率，或在本例中，4000Hz的。 相应的数字频率为1.0。 这两个转折频率则四千分之三百和4000分之3100。

现在让我们说你想要的阻带将下降30分贝+/-从转角频率100赫兹。 因此，你的阻带将在导致4000分之200和四千分之三千二百数字频率200和3200赫兹开始。

要创建过滤器，你会打电话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)

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

文章来源: How to implement band-pass Butterworth filter with Scipy.signal.butter