问题描述
我试图以 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,1,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,gain=desired,fs=fs)
bl = signal.firls(numtaps=narrowtaps,bands=f_bands,desired=desired,weight=[1,2],fs=fs)
br = signal.remez(numtaps=narrowtaps,desired=desired[::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,fs=fs)
bi2,ai2 = signal.cheby1(n,3,wn = signal.cheb2ord(wp,fs=fs)
bi3,ai3 = signal.cheby2(n,30,wn = signal.ellipord(wp,fs=fs)
bi4,ai4 = signal.ellip(n,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,fs=fs)
sos2 = signal.cheby1(n,3,wn = signal.cheb2ord(wp,fs=fs)
sos3 = signal.cheby2(n,30,wn = signal.ellipord(wp,fs=fs)
sos4 = signal.ellip(n,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 一起使用。