带通的butterworth过滤器频率scipy

Question

我正在按照食谱设计scipy中的带通滤波器.但是,如果我过多地降低滤波频率,我会在高阶滤波器中使用垃圾.我究竟做错了什么？

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

1. 不要b, a = butter用于高阶滤波器，无论是在 Matlab、SciPy 还是 Octave 中。传递函数格式存在数值稳定性问题，因为有的系数很大，有的很小。这就是为什么我们在内部将滤波器设计功能更改为使用 zpk 格式。要看到这样做的好处，您需要使用z, p, k = butter(output='zpk')然后使用极点和零点而不是分子和分母。
2. 不要单级做高阶数字滤波器。无论您在什么软件或硬件上实施它们，这都是一个坏主意。通常最好将它们分解为二阶部分。在 Matlab 中，您可以使用zp2sos自动生成这些。在SciPy的，你可以使用sos = butter(output='sos')，然后筛选使用sosfilt()或sosfiltfilt()。这是对大多数应用程序推荐的过滤方式。

Answer 2

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

吉图布