Second Order Butterworth Filter

Python. To grasp its fundamentals, I implemented the second order Butterworth filter for 500 Hz cutoff frequency by calculating its coefficients (pen and paper). Of course, it is so easy to use scipy butter function to get coefficients, particularly for higher orders. 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 cutoff frequency (500Hz), and fs is sampling rate (44100Hz):

(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*}$

(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
    #By pen and paper calculation for 500 hz:
    # 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] += 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')

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

#By pen and paper calculation for 500 hz:

# 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] += 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)