How to design a digital Butterworth bandpass filter?

Question

I am looking into designing a Bandpass Butterworth filter in python, but, I was not sure I am designing my filter correctly. What I have are the following:

• High cutoff frequency = 200Hz
• Low cutoff frequency = 10Hz
• Sampling frequency = 1000Hz
• for my data, I used Filter order = 6

My code definition are below:

# section of my imports:
from scipy.signal import find_peaks, find_peaks_cwt, argrelextrema, welch, lfilter, butter, savgol_filter, medfilt, freqz, filtfilt
from scipy.signal import argrelextrema, filtfilt, butter, lfilter

def butter_bandpass(lowcut, highcut, fs, order):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='bandpass', output='ba')
    # sos = butter(order, [low, high], btype='bandpass', output='sos')
    return b, a
    # return sos

def butter_bandpass_filter(data, lowcut, highcut, fs, order):
    # sos = butter_bandpass(lowcut, highcut, fs, order=order)
    # y = signal.sosfilt(sos=sos, x=data)
    # y = signal.sosfiltfilt(sos=sos, x=data)
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = filtfilt(b=b, a=a, x=data)
    # y = lfilter(b, a, data)
    return y

How can I get the passband and stopband attenuation, also, where can I find the required equations to use in order for me to get my Butterworth filter design equation |H(w)|? Similar to the following link: (Bandpass and Bandstop Filter Design). I calculated the digital frequencies in radians per second:

• wh = 400π rad/sec
• wl = 20π rad/sec
• w(ah) ≈ 21.93 rad/sec
• w(al) ≈ 1.096 rad/sec
• W ≈ 20.84 rad/sec
• w^2 ≈ 578.53

Last steps are the prototype transformation from lowpass-to-bandpass and transforming the equation into Bilinear Transformation Technique (BLT) to get the digital filter are missing. So, what equation do I need to get the digital filter?

Answer 1

In python the direct command is scipy.signal.butter. This will return the filter coefficients (numerator and denominator) based on an array of critical frequencies as described here:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html

Once you have the numerator and denominator coefficients you can use sicpy.signal.freqz to evaluate the frequency response:

import scipy.signal as sig
w, h = sig.freqz(num, den)

(Then simply plot w vs h typically as $20\log_{10}(|h|)$ to view the magnitude in dB, along with the angle of h to view the phase.

freqz simply evaluates the frequency response of the filter (returning $H(z)$ when $z=e^{j\omega}$ (the unit circle on the z-plane). From which we can see the magnitude and phase versus frequency. To use the filter coefficients to filter a time domain signal $x$, use sicpy.signal.lfilter which will provide the convolution of the filter coefficients with the signal to return the filtered result. sig.signal.filtfilt is a "zero-phase" filter which will pass the signal through the filter implementation in both the forward and reverse directions, eliminating the phase component through cancellation in the time reversal, but then also doubling the magnitude response. filtfilt is a non-causal filter that is useful in post-processing applications when we want the output and input to be perfectly in alignment without having to compensate for filter delay between the input and output, but it is not a filter that can be implemented and provide such zero-phase, zero-delay (non-causal).

Note that I am of the opinion that digital filters when mapped from analog prototypes such as this are typically inferior to direct digital designs with FIR filters using optimized algorithms (such as that provided by scipy.signal.firls and scipy.signal.firpm), other than being useful exercise for educational purposes or when modelling an analog system. This point may be my own personal myth, so posted that specifically as another question here.