在scipy中使用Butterworth带通滤波器的频率

Question

我正在使用scipy设计带通滤波器，参考了cookbook。然而，如果我将滤波频率减小太多，高阶滤波器会出现垃圾值。我做错了什么？

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

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 = 25
    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [1, 3, 5, 6, 9]:
        b, a = butter_bandpass(0.5, 4, fs, order=order)
        w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000))
        plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)  
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    plt.figure(2)
    plt.clf()
    for order in [1, 3, 5, 6, 9]:
        b, a = butter_bandpass(0.05, 0.4, fs, order=order)
        w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000))
        plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)  
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    plt.show()

更新：该问题已经在 Scipy 0.14 上得到讨论和解决。然而，Scipy 更新后绘图仍然看起来非常糟糕。出了什么问题？

- Fra

Answer 1

1. 在Matlab、SciPy或Octave中，不要使用b, a = butter进行高阶滤波器设计。该传递函数格式存在数值稳定性问题，因为一些系数非常大而另一些系数则非常小。这就是为什么我们改变了滤波器设计函数内部使用zpk格式的原因。为了看到这种好处，您需要使用z, p, k = butter(output='zpk')，然后使用极点和零点而不是分子和分母。
2. 不要在单个阶段中进行高阶数字滤波器处理，无论您使用的是哪种软件或硬件。这通常是一个坏主意。通常最好将它们分成二阶节。在Matlab中，您可以使用zp2sos自动生成它们。在SciPy中，您可以使用sos = butter(output='sos')，然后使用sosfilt()或sosfiltfilt()进行过滤。这是大多数应用程序推荐的过滤方法。

- endolith

---

谢谢。问题在于freqz和lfilter都需要a、b输入，而不是z、p、k，所以问题是：我该如何使用z、p、k来过滤信号？ - Fra

需要更新freqz以处理高阶滤波器。我认为没有办法使用它。关于过滤和绘制频率响应（和freqz相同的东西），请查看我在链接中发布的butter_sos_example.py。 - endolith

作为更新，scipy 0.16.0 将 'sos' 添加为 butter 的输出选项。freqz_zpk 和 sosfreqz 现在也已经可用。 - Static Void

Answer 2

显然，这个问题是一个已知的bug:

Github

- Fra

---

这个 bug 已经通过使用 ZPK 转换和添加新的函数来进行 SOS 过滤来修复。 - endolith

Answer 3

这是数字滤波器中的一个常见问题。截止频率远低于奈奎斯特频率的高阶滤波器由于浮点数精度有限而倾向于具有不稳定系数。据我所知（虽然是几年前的事情），Matlab在保留精度方面做得比scipy好得多，尽管它仍然会在极端滤波器的情况下出现问题。

如果您不能使用matlab，则有几个选项。第一个选项是将您的滤波器分解成级联的二阶部分。基本上，您计算所需的极点和零点，将它们分解成共轭复对，并为每个对计算传递函数。

第二个选项是重新采样到与滤波器频率更相似的采样率。例如，在您的第二个示例中，您的采样率为25，最高截止频率为0.4。您可以使用低通抗混叠滤波器，然后按10的倍数降采样到2.5的采样率。通过较低的采样率，您的带通滤波器系数将对四舍五入误差的影响更小。如果这样做，您必须确保抗混叠滤波器没有同样的问题。

- Evan

---

也许我误解了你的建议，但我担心这样会引入别名。例如，你在10Hz处有噪声，它低于原始采样的奈奎斯特频率，但如果你将其降采样到2.5Hz，它将被混叠。 - Fra

是的，作为抽取的一部分，您必须使用抗混叠滤波器。我已经澄清了这一点。 - Evan

你能举个关于低通抗混叠滤波器的实际例子吗？ - Fra

即使您正在使用Matlab，也需要使用zp2sos函数将您的滤波器分解为二阶节。 - endolith

Answer 4

发生的情况是，脚本中创建的带通（BP）滤波器的顺序实际上是图中显示顺序的两倍。请记住，滤波器的顺序是传递函数分母中多项式的顺序。标准的带通滤波器始终是偶数阶。

所显示的数字是低通（LP）原型的顺序（通常归一化为1 rad/s的截止频率），该原型用于应用LP到BP变换，从而使滤波器的阶数加倍。因此，例如，如果我们从一阶LP开始，则最终得到二阶带通滤波器：

1/(S+1) => LP-2-BP transf. => k.s/(s^2+a.s+b)

其中k、a和b是常数。标准带通滤波器的分子是k.s^(N/2)，因此滤波器的阶数N必须是偶数。

这个与带通滤波器（也适用于陷波或带阻滤波器）相关的顺序问题在 SciPy 文档 中没有提到。实际上，如果你在 plt.show() 之前使用 print(len(a)) 打印分母 a 的长度，你会发现它有19个系数，对应于一个18阶的多项式。

- user29690