我在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更新之后，剧情看起来还是很糟糕。怎么了？
4条答案
按热度按时间mzaanser1#
1.无论是在Matlab、SciPy还是Octave中，高阶滤波器都不要使用
b, a = butter。传递函数格式存在数值稳定性问题，因为有些系数非常大，而有些系数非常小。这就是我们更改滤波器设计函数to use zpk format internally的原因。要了解其优点，您需要使用z, p, k = butter(output='zpk')，然后处理极点和零点，而不是分子和分母。1.不要在单级中实现高阶数字滤波器。无论你在什么软件或硬件上实现它们，这都是一个坏主意。通常最好将它们分解为second-order sections。在Matlab中，你可以使用zp2sos来自动生成这些。在SciPy中，您可以使用
sos = butter(output='sos')，然后使用sosfilt()或sosfiltfilt()进行筛选。这是大多数应用程序的推荐筛选方式。b91juud32#
显然，该问题是一个已知的错误：
Github
6fe3ivhb3#
这是数字滤波器中的一个常见问题。由于浮点数的精度有限，截止频率远低于奈奎斯特频率的高阶滤波器往往具有不稳定的系数。上一次我检查（承认是几年前）Matlab在保持精度方面做得比scipy好得多，尽管它仍然会给足够极端的滤波器给予问题。
如果不能使用matlab，有几种方法可以选择。第一种方法是将滤波器分解为级联的二阶部分。基本上，计算所需的极点和零点，将它们分解为复共轭对，并计算每对的传递函数。
第二个选项是重新采样到与滤波器频率更接近的采样率。例如，在第二个示例中，采样率为25，最高截止频率为0.4。可以使用低通抗混叠滤波器，然后以10的系数抽取到采样率2.5。使用较低的采样率，带通滤波器系数对舍入误差的敏感度会降低。如果这样做，则必须确保抗混叠滤波器不会出现同样的问题。
thigvfpy4#
实际情况是，脚本中创建的带通（BP）滤波器的阶数实际上是图中所示阶数的两倍。请记住，滤波器的阶数是传递函数分母中多项式的阶数。正则带通滤波器始终为偶数阶。
所示数字是低通（LP）原型的阶数（通常归一化为1 rad/s的截止频率），用于应用LP至BP变换，使滤波器阶数加倍。因此，例如，如果我们以1阶LP开始，则以二阶带通结束：
其中，k、a 和 B 是常数，标准带通滤波器的分子是 k.s^（N/2），因此滤波器的阶数 N 必须是偶数。
带通滤波器的阶数问题（陷波滤波器或带阻滤波器也是如此）在SciPy文档中没有提及。事实上，如果你在
plt.show()之前打印分母a的长度（使用print(len(a))），你会看到它有19个系数，这对应于一个18阶多项式。