Fast and narrow bandpass digital filter implementation in python

I was trying to create [0.5Hz, 2.0Hz] bandpass filter to data with at a sampling rate of 1000Hz. However, when I create FIR filter (scipy.signal.firls), numtaps is too huge (4001). And when I used this FIR with scipy.signal.filtfilt, it took 20 to 30 seconds. It was useless. So I am of trying using an IIR filter, but I'm having trouble because the filter response looks weird.

My code is:

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()

Please tell me if scipy.signal.freqz is buggy like MATLAB, or if my code is weird.

It isn't so much that scipy.freqz or MATLAB is buggy, but rather that the "BA" representation is fairly sensitive to coefficient quantization errors, and especially so as the number of coefficient increases.

Switching to the "SOS" representation avoids this issue:

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)

You could then use the filter with either scipy.signal.sosfilt or scipy.signal.sosfiltfilt.