python中快速和窄带通数字滤波器的实现

Question

我试图为采样率为 1000Hz 的数据创建 [0.5Hz, 2.0Hz] 带通滤波器。 但是，当我创建 FIR 滤波器 (scipy.signal.firls) 时，numtaps 太大了 (4001)。 当我将此 FIR 与 scipy.signal.filtfilt 一起使用时，需要 20 到 30 秒。 这是没用的。 所以我正在尝试使用 IIR 滤波器，但我遇到了问题，因为滤波器响应看起来很奇怪。

我的代码是：

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

narrowtaps = 4001

fs = 1000
dt = 1 / fs
nyq = fs * 0.5
fl = 0.5
fu = 2.0
transwidth = 0.2
f_band = np.array([fl, fu]) / nyq
f_bands = [0, fl - transwidth, fl, fu, fu + transwidth, nyq]
wp = [fl, fu]
ws = [fl - transwidth, fu + transwidth]
desired = [0, 0, 1, 1, 0, 0]

att = signal.kaiser_atten(narrowtaps, transwidth / nyq)
beta = signal.kaiser_beta(att)
bf1 = signal.firwin(numtaps=narrowtaps, cutoff=f_band, window=('kaiser', beta), pass_zero='bandpass')
bf2 = signal.firwin2(numtaps=narrowtaps, freq=f_bands, window=('kaiser', beta), gain=desired, fs=fs)
bl = signal.firls(numtaps=narrowtaps, bands=f_bands, desired=desired, weight=[1, 1, 2], fs=fs)
br = signal.remez(numtaps=narrowtaps, bands=f_bands, desired=desired[::2], weight=[1, 1, 2], fs=fs)

n, wn = signal.buttord(wp, ws, gpass, gstop, fs=fs)
bi1, ai1 = signal.butter(n, wn, 'bandpass', output='ba', fs=fs)
n, wn = signal.cheb1ord(wp, ws, gpass, gstop, fs=fs)
bi2, ai2 = signal.cheby1(n, 3, wn, 'bandpass', output='ba', fs=fs)
n, wn = signal.cheb2ord(wp, ws, gpass, gstop, fs=fs)
bi3, ai3 = signal.cheby2(n, 30, wn, 'bandpass', output='ba', fs=fs)
n, wn = signal.ellipord(wp, ws, gpass, gstop, fs=fs)
bi4, ai4 = signal.ellip(n, 3, 30, wn, 'bandpass', output='ba', fs=fs)

# filter comparison
w1, h1 = signal.freqz(bf1)
w2, h2 = signal.freqz(bf2)
w3, h3 = signal.freqz(bl)
w4, h4 = signal.freqz(br)

w5, h5 = signal.freqz(bi1, ai1)
w6, h6 = signal.freqz(bi2, ai2)
w7, h7 = signal.freqz(bi3, ai3)
w8, h8 = signal.freqz(bi4, ai4)

x_min, x_max = 0, 4
y_min, y_max = -40, 0.5
plt.plot(w1 / np.pi * nyq, 20 * np.log10(np.abs(h1)), label='firwin')
plt.plot(w2 / np.pi * nyq, 20 * np.log10(np.abs(h2)), label='firwin2')
plt.plot(w3 / np.pi * nyq, 20 * np.log10(np.abs(h3)), label='firls')
plt.plot(w4 / np.pi * nyq, 20 * np.log10(np.abs(h4)), label='remez')
plt.plot(w5 / np.pi * nyq, 20 * np.log10(np.abs(h5)), label='butter')
plt.plot(w6 / np.pi * nyq, 20 * np.log10(np.abs(h6)), label='cheby1')
plt.plot(w7 / np.pi * nyq, 20 * np.log10(np.abs(h7)), label='cheby2')
plt.plot(w8 / np.pi * nyq, 20 * np.log10(np.abs(h8)), label='ellip')
plt.grid()
plt.legend()
plt.xlim([x_min, x_max])
plt.ylim([y_min, y_max])
plt.show()

请告诉我scipy.signal.freqz 是否像 MATLAB 一样有问题，或者我的代码是否很奇怪。

Answer 1

与其说是scipy.freqz 或 MATLAB 有问题，不如说是“BA”表示对系数量化误差相当敏感，尤其是当系数数量增加时。

切换到“SOS”表示可以避免这个问题：

gpass = 1.5
gstop = 30
n, wn = signal.buttord(wp, ws, gpass, gstop, fs=fs)
sos1 = signal.butter(n, wn, 'bandpass', output='sos', fs=fs)
n, wn = signal.cheb1ord(wp, ws, gpass, gstop, fs=fs)
sos2 = signal.cheby1(n, 3, wn, 'bandpass', output='sos', fs=fs)
n, wn = signal.cheb2ord(wp, ws, gpass, gstop, fs=fs)
sos3 = signal.cheby2(n, 30, wn, 'bandpass', output='sos', fs=fs)
n, wn = signal.ellipord(wp, ws, gpass, gstop, fs=fs)
sos4 = signal.ellip(n, 3, 30, wn, 'bandpass', output='sos', fs=fs)

N = 16384
w9,  h9  = signal.sosfreqz(sos1, worN=N)
w10, h10 = signal.sosfreqz(sos2, worN=N)
w11, h11 = signal.sosfreqz(sos3, worN=N)
w12, h12 = signal.sosfreqz(sos4, worN=N)

然后您可以将过滤器与scipy.signal.sosfilt 或scipy.signal.sosfiltfilt 一起使用。