This is an example on how to design a filter in the analog domain, and then use the bilinear transform to transform it to the digital domain, while preserving the cut-off frequency.

We'll be using formulas derived on the Bilinear Transform and Butterworth Filters pages.

Design criteria

In this example, we'll design a digital fourth order Butterworth low-pass filter, with a sample frequency of 360Hz and a cut-off frequency of 45Hz.

Frequency Pre-Warping

As discussed in the page on the Bilinear Transform, we have to apply pre-warping to the cut-off frequency before designing a filter. If we don't the cut-off frequency will shift to an incorrect frequency when we discretize the filter.

First, let's calculate the normalized digital frequency ωc,d, using the cut-off frequency fc and the sample frequency fs: fc=45Hzfs=360HzTs=1fsωc=2πfc282.7rads1ωc,d=ωcfs=π4radsample0.7854radsample The Nyquist-Shannon sampling theorem tells us that we can never sample frequencies higher than fs/2 without losing information. This also means that the cut-off frequency can never be higher than half of the sample frequency. Or in other words, all normalized frequencies will be in the interval [0,π].

Next, we'll use the pre-warping formula we derived in the page on the Bilinear Transform, in order to calculate the analog design frequency ωc,a: ωc,a=2Tstan(ωc,d2)=720tan(π8)298.2rads1 Note that this frequency is relatively close to ωc, but it is not the same. The higher the cut-off frequency (relative to the sample frequency), the larger the error between ωc and ωc,a.

Designing the Butterworth filter in the Analog Domain

Now that we know the pre-warped analog cut-off frequency, we can start designing the analog filter.
We'll use the formula for the Butterworth low-pass filter derived in the page on Butterworth Filters: (1)H4(s)=1B4(s)wheressωc,a B4(s)=(s22cos(2π4+144)s+1)(s22cos(2π2+4+144)s+1) Defining these constants will make the calculations much easier: (2)α2cos(5π8)=22(3)β2cos(7π8)=2+2 (4)B4(s)=(s2+αs+1)(s2+βs+1)=s4+s3(α+β)+s2(αβ+2)+s(α+β)+1

Discretizing the Analog Filter

We can now just apply the Bilinear Transform to the analog transfer function, by substituting s=2Tsz1z+1. Therefore: s=2fsωa,cz1z+1 Again, we'll introduce a constant to simplify the expression: (5)γ2fsωa,c=2fs2fstan(ωc,d2)=cot(πfcfs)=cot(π8)=1+2 (6)s=γz1z+1 What follows is just rearranging the expression of B4(s) from Equation 4, using the substitution of Equation 6.
Finally, we end up with an expression for the transfer function, using Equation 1, and we can determine the coefficients using the constants defined in Equations 2, 3 & 5. B4(s)=s4+s3(α+β)+s2(αβ+2)+s(α+β)+1B4(z)=γ4(z1)4(z+1)4+γ3(z1)3(z+1)(z+1)4(α+β)+γ2(z1)2(z+1)2(z+1)4(αβ+2)+γ(z1)(z+1)3(z+1)4(α+β)+(z+1)4(z+1)4=1(z+1)4[γ4(z1)4+γ3(z1)3(z+1)(α+β)+γ2(z1)2(z+1)2(αβ+2)+γ(z1)(z+1)3(α+β)+(z+1)4]=1(z+1)4[γ4(z44z3+6z24z+1)+γ3(z42z3+2z1)(α+β)+γ2(z42z2+1)(αβ+2)+γ(z4+2z32z1)(α+β)+(z4+4z3+6z2+4z+1)]=1(z+1)4[(γ4+γ3(α+β)+γ2(αβ+2)+γ(α+β)+1)z4+(4γ42γ3(α+β)+2γ(α+β)+4)z3+(6γ42γ2(αβ+2)+6)z2+(4γ4+2γ3(α+β)2γ(α+β)+4)z+(γ4γ3(α+β)+γ2(αβ+2)γ(α+β)+1)]H4(z)=1B4(z)=(z+1)4[(γ4+γ3(α+β)+γ2(αβ+2)+γ(α+β)+1)z4+(4γ42γ3(α+β)+2γ(α+β)+4)z3+(6γ42γ2(αβ+2)+6)z2+(4γ4+2γ3(α+β)2γ(α+β)+4)z+(γ4γ3(α+β)+γ2(αβ+2)γ(α+β)+1)]b0+b1z1+b2z2+b3z3+b4z4a0+a1z1+a2z2+a3z3+a4z4 b0=1b1=4b2=6b3=4b4=1a0=γ4+γ3(α+β)+γ2(αβ+2)+γ(α+β)+197.95a1=4γ42γ3(α+β)+2γ(α+β)+4192.8a2=6γ42γ2(αβ+2)+6170.0a3=4γ4+2γ3(α+β)2γ(α+β)+470.96a4=γ4γ3(α+β)+γ2(αβ+2)γ(α+β)+111.79

Frequency Response & Pole-Zero Map

We can check the filter's frequency response to make sure that we didn't make any mistakes. As mentioned in other pages, the frequency response of a digital system can be a obtained by evaluating the transfer function H(z) along the unit circle (z=ejω). We'll plot the magnitude in decibels. A(ω)=20log10|H(ejω)| We can also plot the phase angle of the response: ϕ(ω)=H(ejω) You can see that the corner frequency lies around 45Hz. We can check this mathematically: A(ωc,d)=3.01dB

MATLAB & GNU Octave

If you have to design many different filters, you don't want to calculate them all by hand. Luckily, MATLAB and GNU Octave come with a command to calculate the coefficients of Butterworth filters.

f_s = 360;    % Sample frequency in Hz
f_c = 45;     % Cut-off frequency in Hz
order = 4;    % Order of the butterworth filter

omega_c = 2 * pi * f_c;       % Cut-off angular frequency
omega_c_d = omega_c / f_s;    % Normalized cut-off frequency (digital)

[b, a] = butter(order, omega_c_d / pi);    % Design the Butterworth filter
disp('a = '); disp(a);                     % Print the coefficients
disp('b = '); disp(b);
H = tf(b, a, 1 / f_s);                     % Create a transfer function
bode(H);                                   % Show the Bode plot

Note that MATLAB expects the normalized frequency as a number from 0 to 1, so we have to divide by π before passing it to the butter function.

Python

A similar function is available in the SciPy signal package: butter.

from scipy.signal import butter, freqz, freqs
import matplotlib.pyplot as plt
from math import pi
import numpy as np

f_s = 360    # Sample frequency in Hz
f_c = 45     # Cut-off frequency in Hz
order = 4    # Order of the butterworth filter

omega_c = 2 * pi * f_c       # Cut-off angular frequency
omega_c_d = omega_c / f_s    # Normalized cut-off frequency (digital)

# Design the digital Butterworth filter
b, a = butter(order, omega_c_d / pi)    
print('Coefficients')
print("b =", b)                           # Print the coefficients
print("a =", a)

w, H = freqz(b, a, 4096)                  # Calculate the frequency response
w *= f_s / (2 * pi)                       # Convert from rad/sample to Hz

# Plot the amplitude response
plt.subplot(2, 1, 1)            
plt.suptitle('Bode Plot')
H_dB = 20 * np.log10(abs(H))              # Convert modulus of H to dB
plt.plot(w, H_dB)
plt.ylabel('Magnitude [dB]')
plt.xlim(0, f_s / 2)
plt.ylim(-80, 6)
plt.axvline(f_c, color='red')
plt.axhline(-3, linewidth=0.8, color='black', linestyle=':')

# Plot the phase response
plt.subplot(2, 1, 2)
phi = np.angle(H)                         # Argument of H
phi = np.unwrap(phi)                      # Remove discontinuities 
phi *= 180 / pi                           # and convert to degrees
plt.plot(w, phi)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [°]')
plt.xlim(0, f_s / 2)
plt.ylim(-360, 0)
plt.yticks([-360, -270, -180, -90, 0])
plt.axvline(f_c, color='red')

plt.show()

Comparison Between the Analog and Digital Filter

We can easily plot the Bode plots of the two filters on top of each other, in order to compare their properties.

from scipy.signal import butter, freqz, freqs
import matplotlib.pyplot as plt
from math import pi
import numpy as np

f_s = 360    # Sample frequency in Hz
f_c = 45     # Cut-off frequency in Hz
order = 4    # Order of the butterworth filter

omega_c = 2 * pi * f_c       # Cut-off angular frequency
omega_c_d = omega_c / f_s    # Normalized cut-off frequency (digital)

# Design the digital Butterworth filter
b_d, a_d = butter(order, omega_c_d / pi)    
print('Digital Coefficients')
print("b =", b_d)                         # Print the coefficients
print("a =", a_d)

w, H_d = freqz(b_d, a_d, 4096)            # Calculate the frequency response
w *= f_s / (2 * pi)                       # Convert from rad/sample to Hz

# Design the analog Butterworth filter
b_a, a_a = butter(order, f_c, analog=True)
print('Analog Coefficients')
print("b =", b_a)                         # Print the coefficients
print("a =", a_a)

w, H_a = freqs(b_a, a_a, w)               # Calculate the frequency response

# Plot the amplitude response
plt.subplot(2, 1, 1)            
plt.suptitle('Bode Plot')
H_d_dB = 20 * np.log10(abs(H_d))          # Convert modulus of H to dB
H_a_dB = 20 * np.log10(abs(H_a))
plt.plot(w, H_d_dB, color='blue', label='Digital')
plt.plot(w, H_a_dB, color='green', label='Analog')
plt.legend()
plt.ylabel('Magnitude [dB]')
plt.xlim(0, f_s / 2)
plt.ylim(-80, 6)
plt.axvline(f_c, color='red')
plt.axhline(-3, linewidth=0.8, color='black', linestyle=':')

# Plot the phase response
plt.subplot(2, 1, 2)
phi_d = np.angle(H_d)                     # Argument of H
phi_a = np.angle(H_a)
phi_d = np.unwrap(phi_d) * 180 / pi       # Remove discontinuities
phi_a = np.unwrap(phi_a) * 180 / pi       # and convert to degrees
plt.plot(w, phi_d, color='blue')
plt.plot(w, phi_a, color='green')  
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [°]')
plt.xlim(0, f_s / 2)
plt.ylim(-360, 0)
plt.yticks([-360, -270, -180, -90, 0])
plt.axvline(f_c, color='red')

plt.show()

Discretization using Second Order Sections (SOS)

For higher order filters, small quantization errors on the transfer function coefficients can result in large errors on the pole and zero locations.
A solution is to factor the transfer function into second order factors or sections.

Recall from the previous sections, and define αk, H2,k and H1 as follows: (7)γcot(πfcfs)(8)ssωc,a=γz1z+1(9)αk2cos(2π2k+n+14n)(10)H2,k(s)1B2,k(s)1s2αks+1(11)H1(s)1s+1(12)Hn(s)={k=0n21H2,k(s)even nH1(s)k=0n121H2,k(s)odd n

Second Order Sections

We'll use the same technique as before to substitute s into B2,k(s) using the pre-warped bilinear transform relation 8 to get the discrete-time Butterworth polynomial B2,k(z): B2,k(s)s2αks+1B2,k(z)=γ2(z1z+1)2αkγz1z+1+1=γ2(z1)2(z+1)2αkγ(z1)(z+1)(z+1)2+(z+1)2(z+1)2=γ2(z1)2αkγ(z1)(z+1)+(z+1)2(z+1)2=(γ2αkγ+1)z2(2γ22)z+(γ2+αkγ+1)(z+1)2H2,k(z)=z2+2z+1(γ2αkγ+1)z2(2γ22)z+(γ2+αkγ+1)=1+2z1+z2(γ2αkγ+1)(2γ22)z1+(γ2+αkγ+1)z2 The coefficients of the k-th factor of the discrete-time transfer function are thus: bk,0=1bk,1=2bk,2=1ak,0=γ2αkγ+1ak,1=2(1γ2)ak,2=γ2+αkγ+1

First Order Section

For odd orders n, we need n12 second order sections and a single first order section.

Again, we'll use the the pre-warped bilinear transform relation 8: H1(s)1s+1=1γz1z+1+1=z+1γ(z1)+(z+1)=z+1(γ+1)z+(1γ)=1+z1(γ+1)+(1γ)z1 This gives us the coefficients of the first order factor of the discrete-time transfer function: b0=1b1=1a0=γ+1a1=1γ