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

【问题标题】：Bandpass butterworth filter frequencies in scipyscipy中的带通巴特沃斯滤波器频率
【发布时间】：2014-03-18 17:42:51
【问题描述】：

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

【问题讨论】：

标签： python numpy scipy filtering

【解决方案1】：

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

【讨论】：

谢谢。问题是freqz和lfilter似乎都需要a、b输入，而不是z、p、k，所以问题是：我如何实际使用z、p、k来过滤信号？跨度>
是的，freqz 需要更新以处理高阶过滤器。我不认为有办法使用它。对于过滤和绘制频率响应（与 freqz 相同），请参阅我发布的链接中的 butter_sos_example.py
作为更新，scipy 0.16.0 添加了“sos”作为黄油的输出选项。 freqz_zpk 和 sosfreqz 现在也可用。

【解决方案2】：

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

Github

【讨论】：

此错误已通过使用 ZPK 转换和添加新功能进行 SOS 过滤得到修复

【解决方案3】：

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

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

第二种选择是重新采样到更接近滤波器频率的采样率。例如，在您的第二个示例中，您的采样率为 25，您的最高截止频率为 0.4。您可以使用低通抗混叠滤波器，然后以 10 倍抽取以达到 2.5 的采样率。使用较低的采样率，您的带通滤波器系数将对舍入误差不太敏感。如果你这样做，你必须确保抗锯齿过滤器没有同样的问题。

【讨论】：

也许我误解了你的建议，但恐怕你会以这种方式引入别名。例如。你有 10Hz 的噪声，它低于原始采样的奈奎斯特频率，但如果你抽取到 2.5Hz，它会混叠
是的，您必须使用抗混叠滤波器作为抽取的一部分。我澄清了这一点。
你能举一个低通抗混叠滤波器的实际例子吗？
即使您正在使用 Matlab，也需要使用 zp2sos 函数将过滤器分成二阶部分。

【解决方案4】：

脚本中创建的带通 (BP) 滤波器的阶数实际上是图中所示的阶数的两倍。回想一下，滤波器的阶数是传递函数分母中多项式的阶数。规范带通滤波器始终是偶数阶。

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

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

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

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

【讨论】：