我正在按照 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 更新后,情节看起来仍然很糟糕。怎么了?
最佳答案
- 不要将
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()
.这是为大多数应用推荐的过滤方式。
关于python - scipy 中的带通 Butterworth 滤波器频率,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21862777/