Second Order Butterworth Filter

The Laplace transform and z-transform are mathematical ways of breaking an impulse response into sinusoids and decaying exponentials. This is done by expressing the system’s characteristics as one complex polynomial divided by another complex polynomial. Practical Infinite-Impulse-Response (IIR) filters are usually based upon analogue equivalents (Butterworth, Chebyshev, etc.), using a transformation known as the bilinear transformation which maps the s-plane poles and zeros of the analogue filter into the z-plane. The bilinear transform is a mathematical transformation from the z-domain to the s-domain where T is sampling period, fc is cut-off frequency and fs is sampling rate:

(1) $\begin{equation*} s = \frac{2}{T}\frac{(1-z^{-1})}{(1+z^{-1})} \end{equation*}$

(2) $\begin{equation*} w_{c} = \frac{2}{T}tan\frac{\pi f_{_{c}}}{f_{s}} \end{equation*}$

(3) $\begin{equation*} G(s) = \frac{w{_{c}}^{2}}{s^{2} + s\sqrt{2}w_{c} + w{_{c}}^{2}} \end{equation*}$

def second_order_butterworth_coefficients(sample_rate: float, frequency: float):
    """
    Calculate second-order Butterworth filter coefficients.

    Args:
        sample_rate (float): The sampling rate.
        frequency (float): The cutoff frequency.

    Returns:
        list[float]: Coefficients [b0, b1, b2, a0, a1, a2].
    """
    # Constants
    pi = math.pi
    Q = 0.70711

    # Calculations
    n = 1 / math.tan(pi * frequency / sample_rate)
    n_squared = n * n
    inv_Q = 1 / Q
    c = 1 / (1 + inv_Q * n + n_squared)  # Scaling factor

    # Coefficients [b0, b1, b2, a0, a1, a2]
    return [
        c,                         # b0
        2 * c,                     # b1
        c,                         # b2
        1,                         # a0
        2 * c * (1 - n_squared),   # a1
        c * (1 - inv_Q * n + n_squared)  # a2
    ]

Let’s choose 500Hz cut-off frequency at 44100Hz sampling rate. Calculating Second Order Butterworth coefficients at real time is not computationally expensive when compared to higher orders.

(4) $\begin{equation*} G(z) = 0.00121\frac{{z^{2} + 2z + 1}}{{z^{2} - 1.89933z + 0.90416}} \end{equation*}$

import os
import numpy as np
import scipy.io.wavfile as wavfile

def ButterworthLowPass(fileName):
 
    #fs: sampling rate of the wav file, x: integer point array
        
    if (os.path.isfile(fileName) == False):					
        raise ValueError("Input file is not valid")
            
    fs, x = wavfile.read(fileName)
 
    #We only accept mono and 44.1khz audio files for this implementation  
    
    if (len(x.shape) != 1):                                  
            raise ValueError("Audio file is not mono")
    if (fs != 44100):                                        
            raise ValueError("Sampling rate of input sound is not 44100")   

    #The second order Butterworth filter coefficients
    #Precomputed values for 500 hz (hard-coding)
    # b: {0.00121, 0.00242, 0.00121}, a: {1, -1.89933, 0.90416}
    
    y = np.zeros_like(x)   
    
    for n in range(2, len(x)):
        y[n] = 0.00121 * x[n] + 0.00242 * x[n - 1] + 0.00121 * x[n - 2]
        y[n] = y[n] + 1.89933 * y[n - 1] - 0.90416 * y[n - 2]

    difference = x - y
    
    wavfile.write("Filtered-Audio.wav", 44100, y.astype(np.int16))
    wavfile.write("Difference-Audio.wav", 44100, difference.astype(np.int16))

if __name__ == "__main__":
    
    #Run it with your test audio file
    ButterworthLowPass('C:/Users/[your name]/Desktop/[test].wav')

[1] Butterworth Filter, Wikipedia

[2] The Scientist and Engineer’s Guide to Digital Signal Processing, Steven W. Smith

[3] Design of Digital Filters (Lecture Series, Oxford University)