Written by Niels Moseley
on July 24, 2023

Bilinear transform

The procedure outlined here will transform an analog prototype filter, of which the pole location are know, into a digital version. Tables with pole locations for standard filters, such as Butterworth and Chebyshev responses, can be found on the internet. For example, see: Analog Devices.

If the real part of the pole locations given by the tables are positive numbers, they must be made negative. Furthermore, the prototype filters have a cutoff frequency of 1 rad/s. In order to move the cutoff to the desired frequency, multiply the real and imaginary parts of the pole positions by $2 \pi f_{cut}$, with $f_{cut}$ in Hz.

Note that, because of simplifications, the procedure only works for even-order all-pole IIR lowpass filters.

Generating a digital version

A digital version of the analog prototype filter is constructed using digital second-order sections. Due to the discrete-time nature, these digital filters have different coefficients from the analog one. The bilinear transform is a good method to obtain these coefficients.

The bilinear transforms the left-hand side of the s-plane into the unit circle (and its interior) of the z-domain. This ensures that a stable analog filter results in a stable digital filter.

The mapping is achieved by the following equation: $$z = \dfrac{2\cdot f_s + s}{2\cdot f_s - s},$$ where $s$ is the position in the s-plane, $z$ is the position in the z-plane and $f_s$ is the sample rate in samples per second.

Generating a digital version of the analog prototype filter is done by transforming the pole locations of the analog filter using the mapping function above.

The following Python code implements the bilinear transform mapping function and the cutoff frequency update:

def bilinear(s : complex, fcut : float, fs : float) -> complex:
    sc = s * 2.0*3.1415927*fcut
    z = (2.0*fs + sc) / (2.0*fs - sc)
    return z

where $s$ is the analog pole location, $fcut$ is the desired cutoff frequency in Hz and $fs$ is the sample rate in samples per second.

Second-order digital IIR filter

The second-order digital IIR filter is defined by the complex pole $p_d$ and its complex conjugate ${p_d}^*$ and has the following z-domain transfer function: $$ H(z) = \dfrac{1 + 2 \cdot z^{-1} + 1 \cdot z^{-2} }{1 -2 \Re(p_d) \cdot z^{-1} + |p_d|^2 \cdot z^{-2} },$$ where $p_d$ is the pole location of the second-order filter.

DC response

The gain of the digital filter at DC is equal to $H(z)|_{z=1}$, so: $$ H(1) = \dfrac{4}{1-2 \Re(p_d)+|p_d|^2} $$

Nyquist response

The gain of the digital filter at Nyquist is equal to $H(z)|_{z=-1}$, so: $$ H(-1) = \dfrac{0}{1+2 \Re(p_d)+|p_d|^2} = 0 $$

Example

#!/usr/bin/python3
# Using the bilinear transform to design digital filters.

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

def bilinear(s : complex, fcut : float, fs : float) -> complex:
    sc = s * 2.0*3.1415927*fcut
    z = (2.0*fs + sc) / (2.0*fs - sc)
    return z

# The analog Chebyshev IIR filter has only poles
poles = [-0.4501 + 0.3840j, -0.1865 + 0.9272j]    # 0.25 dB ripple

fc = 3500   # cutoff frequency in Hz
fs = 48000  # sample rate in sample per second

biquads = []
sec = 1
for pa in poles:
    pd = bilinear(pa, fc, fs)

    # biquad denominator coefficients
    a1 = -2*np.real(pd)
    a2 = np.abs(pd)*np.abs(pd)
    a = np.array([1, a1, a2])

    # biquad numerator coefficients
    b = np.array([1.0,2.0,1.0])
    
    # calculate gain of each section
    # and normalize DC to 0 dB
    gain = np.sum(b) / np.sum(a)
    b = b / gain

    biquads.append(b)
    biquads.append(a)

    print("Section", sec)
    print("  Numerator   : ", b)
    print("  Denominator : ", a)
    sec = sec + 1


## plot the response of the filter

w1,h1 = signal.freqz(biquads[0], biquads[1], fs=fs, worN=4096)
w2,h2 = signal.freqz(biquads[2], biquads[3], fs=fs, worN=4096)

plt.figure(1)
plt.semilogx(w1, 20.0*np.log10(np.abs(np.multiply(h1,h2))), [fc,fc], [-100, 10])
plt.title("Chebyshev IIR filter response")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Magnitude [dB]")
plt.xlim([0, fs/2])
plt.ylim([-20, 10])
plt.grid()

plt.figure(2)
plt.semilogx(w1, 20.0*np.log10(np.abs(np.multiply(h1,h2))), [fc,fc], [-100, 10])
plt.title("Chebyshev IIR filter response [zoom]")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Magnitude [dB]")
plt.xlim([0, fs/2])
plt.ylim([-1, 1])
plt.grid()
plt.show()

Output from the Python script


Section 1
  Numerator   :  [0.01499998 0.02999997 0.01499998]
  Denominator :  [ 1.         -1.603211    0.66321093]
Section 2
  Numerator   :  [0.04145069 0.08290138 0.04145069]
  Denominator :  [ 1.         -1.68328578  0.84908853]

← → Top