python 中的快速窄带通数字滤波器实现
Fast and narrow bandpass digital filter implementation in python
我试图为采样率为 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 一样有问题，或者我的代码是否很奇怪。
与其说 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 or scipy.signal.sosfiltfilt.
一起使用
我试图为采样率为 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 一样有问题，或者我的代码是否很奇怪。
与其说 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 or scipy.signal.sosfiltfilt.