scipy中的带通巴特沃斯滤波器频率

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更新之后，情节看起来仍然很糟糕。怎么了？

Answer 1

显然，问题是一个已知的错误：

Github

Answer 2

1. 不要对高阶滤波器使用b, a = butter，无论是在Matlab中，还是在SciPy或Octave中。传递函数格式has numerical stability problems，因为有些系数很大，而有些系数很小。这就是我们更改过滤器设计函数to use zpk format internally的原因。要了解它的好处，您需要使用z, p, k = butter(output='zpk')，然后使用极和零，而不是分子和分母。
2. 不要在一个阶段做高阶数字滤波器。这是一个坏主意，不管你在什么软件或硬件上实现它们。通常最好把它们分成second-order sections。在Matlab中，可以使用zp2sos自动生成这些文件。在SciPy中，可以使用sos = butter(output='sos')，然后使用^{}或^{}进行过滤。对于大多数应用程序，这是推荐的筛选方法。

Answer 3

这是数字滤波器中常见的问题。由于浮点数精度的限制，截止频率远低于奈奎斯特频率的高阶滤波器往往具有不稳定系数。上一次我检查（不可否认几年前）Matlab在保存精度方面比scipy做得好得多，尽管它仍然会给足够极端的滤波器带来问题。

如果你不能使用matlab，有几个选择。第一种方法是将滤波器分解成级联的二阶部分。基本上你要计算你想要的极点和零点，把它们分解成复杂的共轭对，然后计算每对的传递函数。

第二个选项是重新采样到与滤波器频率更相似的采样率。例如，在第二个示例中，采样率为25，最高截止频率为.4。您可以使用低通抗混叠滤波器，然后按10的倍数进行抽取，采样率为2.5。采样率越低，带通滤波器系数对舍入误差的敏感度就越低。如果这样做，则必须确保抗锯齿过滤器没有相同的问题。