2
$\begingroup$

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?

$\endgroup$
8
  • 1
    $\begingroup$ Please edit your code to show the relevant imports (presumably scipy.signal and wherever you're importing butter from). Somewhere in scipy.signal there's a function that'll give you the frequency response of a filter -- you can use that to graph the response. $\endgroup$
    – TimWescott
    Dec 1 '21 at 19:20
  • $\begingroup$ @TimWescott Thank you for the reply. I found in scipy function called freqz. So, I’ll have a look into that. In terms of getting an equation of filter, how do I calculate it? $\endgroup$
    – WDpad159
    Dec 1 '21 at 19:32
  • $\begingroup$ @WDpad159 why are you using filtfilt? Again, please explicitly add the imports you're using, if only to allow anyone to reproduce your results with simple copy and paste. $\endgroup$ Dec 1 '21 at 19:41
  • $\begingroup$ @MarcusMüller See my updated code for the imports $\endgroup$
    – WDpad159
    Dec 1 '21 at 20:00
  • $\begingroup$ "where can I find the required equations" -- do you mean how do you arrive at difference equations to implement the filter outside of Python/numpy, or do you mean that you want the design equations that butter uses? $\endgroup$
    – TimWescott
    Dec 1 '21 at 20:56
1
$\begingroup$

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.

$\endgroup$
13
  • $\begingroup$ There are good reasons to use IIR filters -- occasionally. Usually FIR filters make more sense, but it's good to know how to do it the good old fashioned way. I think that if I were writing IIR filters for something I'd get the poles & zeros and build up my $2^{nd}$ order sections, rather than starting from polynomials and trying to factor them accurately. $\endgroup$
    – TimWescott
    Dec 2 '21 at 1:44
  • $\begingroup$ Yup no disagreement but I don’t know if any good reasons to implement a butterworth vs least squares except for the reasons I stated- got any? $\endgroup$ Dec 2 '21 at 1:47
  • 1
    $\begingroup$ Sorry -- forgot to mention that. I thought you were dissing IIR filters in general, not mapping from the classics (which were basically formulated because they kinda worked, and had closed-form solutions -- even if I were designing a filter with caps and resistors and op-amps I'd fit to the frequency spec, not just pull a Chebychev out of my hat). $\endgroup$
    – TimWescott
    Dec 2 '21 at 5:13
  • 1
    $\begingroup$ @WDpad159 freqz is just for the evaluation. The filter equation is described by the numerator and denominator polynomials expressed as vectors of their coefficients (and for freqz these vectors are in decreasing negative powers of z, as in 1+z^-1+z^-2 etc). To use it to filter a time domain signal x, use sig.lfilter(num, den, x). I'll add that last point to the answer in case that was core to your question. $\endgroup$ Dec 2 '21 at 13:38
  • 1
    $\begingroup$ Yes it is fine to use. It will not give you the frequency magnitude response of the filter but will give you the response you would get if you passed it through your filter twice. That may in fact be fine for you. If you want to see what an IIR filter with the coefficients you derived would actually do to your signal (accurately) use lfilter not filtfilt. $\endgroup$ Dec 2 '21 at 14:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.