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.
