Learn how to derive and use the bilinear transform to convert analog systems into digital ones.

To design and implement a parametric filter, 4 steps are needed:

  1. Decide on the filter type you want to implement.
  2. Design an analog prototype.
  3. Digitize the analog prototype using the bilinear transform.
  4. Implement the digital filter in code.

In this article, we’ll explain the second step of that process: how to go from an analog prototype to a digital form of the parametric filter. This is shown on the diagram below.

Figure 1. In this article, we discuss the bilinear transform.

System Digitization Methods

There are many methods of digitalization (discretization) of analog (continuous) systems. The most popular ones are [OppenheimSchafer10]:

  • impulse-invariant transformation (IIT): sampling the impulse response of the continuous system,
  • bilinear transformation (BT): mapping the transfer function of the continuous system from the ss-plane to the zz-plane.

Because the impulse-invariant transformation may introduce aliasing, in the music domain, we usually use the bilinear transform. So how is it defined?

The Bilinear Transform

The bilinear transform is a change of variables, mapping the jωj \omega frequency axis of the ss-plane to the unit circle of the zz-plane.

Given an analog transfer function Ha(s)H_\text{a}(s) in the ss-domain, we obtain the discrete transfer function Hd(z)H_\text{d}(z) in the zz-domain by substituting

s=2T1z11+z1,(1)s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}, \quad (1)

where T=1/fsT = 1 / f_s is the sampling interval of the discrete system with fsf_s being the sampling rate. Therefore,

Hd(z)=Ha(2T1z11+z1).(2)H_\text{d}(z) = H_\text{a}\left(\frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}\right). \quad (2)

We can go back from the discrete transfer function to the analog transfer function by substituting

z=1sT/21+sT/2,(3)z = \frac{1-sT/2}{1+sT/2}, \quad (3)

what can be directly derived from Equation 1.

The bilinear transform formulas are summarized in Figure 2 (with c=2/Tc=2/T).

Figure 2. Bilinear transform formulas.

Derivation of the Bilinear Transform

How is the bilinear transform derived? How to come up with the mapping in Equation 1?

Below, I give you a simple, intuitive explanation. If you are not interested in it, you can skip to the properties of the bilinear transform.

The szs \rightarrow z mapping of the bilinear transform may be explained as a derivation of a digital representation of an integrator.

An integrator is a system at whose output the integral of the input signal may be observed

y(t)=0tx(τ)dτ,(4)y(t) = \int \limits_0^t x(\tau) d\tau, \quad (4)

where x(t)x(t) is the input signal over continuous time tt and y(t)y(t) is the output signal of the integrator.

The Laplace transform of Equation 4 yields

Y(s)=1sX(s),(5)Y(s) = \frac{1}{s}X(s), \quad (5)

where sCs \in \mathbb{C}. Therefore, the integrator can be shown on the diagram as in Figure 3.

Figure 3. A diagram of the analog integrator system.

We want to obtain a discrete system described by a discrete transfer function G(z)G(z) that behaves like 1s\frac{1}{s}. To this end, we observe the output of the continuous system at discrete time points t=kTt = kT, where kZk \in \mathbb{Z} and TT is the sampling interval (inverse of the sampling rate).

y(kT)=0kTx(τ)dτ=0(k1)Tx(τ)dτ+(k1)TkTx(τ)dτ=y((k1)T)+(k1)TkTx(τ)dτ.(6)y(kT) = \int \limits_{0}^{kT} x(\tau) d\tau = \int \limits_{0}^{(k-1)T} x(\tau) d \tau + \int \limits_{(k-1)T}^{kT} x(\tau) d \tau \newline = y((k-1)T) + \int \limits_{(k-1)T}^{kT} x(\tau) d \tau. \quad (6)

We can now approximate the last integral in Equation 6 via the trapezoidal rule

(k1)TkTx(τ)dτ(x((k1)T)+x(kT))T/2.(7)\int \limits_{(k-1)T}^{kT} x(\tau) d \tau \approx (x((k-1)T) + x(kT))T / 2. \quad (7)

Inserting Equation 7 into Equation 6 yields

y(kT)=y((k1)T)+(x((k1)T)+x(kT))T/2.(8)y(kT) = y((k-1)T) + (x((k-1)T) + x(kT))T / 2. \quad (8)

We now transform Equation 8 with the zz-transformation

Y(z)=z1Y(z)+(z1X(z)+X(z))T/2.(9)Y(z) = z^{-1} Y(z) + (z^{-1} X(z) + X(z))T / 2. \quad (9)

After rearranging Equation 9 we obtain the transfer function formulation

Y(z)=T21+z11z1X(z)=G(z)X(z).(10)Y(z) = \frac{T}{2}\frac{1 + z^{-1}}{1 - z^{-1}} X(z) = G(z) X(z). \quad (10)

Since G(z)G(z) was meant to approximate 1s\frac{1}{s}, we can write

1s=G(z)=T21+z11z1.(11)\frac{1}{s} = G(z) = \frac{T}{2}\frac{1 + z^{-1}}{1 - z^{-1}}. \quad (11)

which after inversion gives Equation 1, i.e., the bilinear transform. \Box

The bilinear transform can also be derived by representing a signal as a series of impulses, calculating its Laplace transform and zz-transform, equating the two, and approximating the resulting ss-to-zz mapping with a series expansion. But the integrator derivation seems more intuitive to me, especially the origin of the 2T\frac{2}{T} constant. What do you think? Let me know in the comments below the article. 😉

Properties of the Bilinear Transform

In the following, we assume that ωa=2πfa\omega_\text{a} = 2 \pi f_\text{a} (with faf_\text{a} in Hz) denotes the analog frequency in radians per second (value on the jωj\omega axis of the ss-plane) and ωd=2πfd\omega_\text{d} = 2 \pi f_\text{d} (with fdf_\text{d} in Hz) denotes the digital frequency in radians per second (corresponding to the angle of a point on the unit circle of the zz-plane).

The bilinear transform has a number of useful properties [ParksBurrus87].

First of all, it is called bilinear because the numerator and the denominator are linear in zz.

Second of all, the bilinear transform maps the left half-plane of the ss-plane into the interior of of the unit circle in the zz-plane. This is shown on Figure 4.

Figure 4. Bilinear transform maps the left half-plane of the ss-plane into the unit circle of the zz-plane. The infinitely long analog frequency axis jωj\omega becomes the finite-length circle. Axis bending in the middle is shown only for the visualization purposes: it is not what is meant by the mathematics.

As a consequence, the poles from the left half-plane of the ss-plane are mapped to the poles within the unit circle of the zz-plane. That means that stable analog filters are transformed into stable digital filters, what is a very desirable property in the context of musically useful parametric filters.

In the same way, the jωaj\omega_\text{a} analog frequency axis from the ss-plane is mapped to the z=ejωdTz = e^{j\omega_\text{d}T} unit circle of the zz-plane. That means that the analog system’s frequency response over an infinite domain gets squashed around a finite-length circle.

The above mappings mean that the bilinear transform is a nonlinear mapping of the phase and the frequency. The former results in the nonlinear phase of the resulting IIR filters [OppenheimSchafer10]. The latter manifests itself in the form of frequency warping as discussed later in the article.

Other properties of the bilinear transform are:

  • the obtained digital filter has the same order as the analog prototype,
  • optimality is preserved (e.g., maximally flat analog prototypes become maximally flat discrete filters),
  • a cascade of systems transformed with the bilinear transform is equivalent to the bilinear transform of a cascade of these systems.

Frequency Warping

The bilinear transform maps the full, infinite analog frequency axis to a finite-length discrete frequency axis (the unit circle). That results in frequency warping: “equal increments along the unit circle in the zz plane correspond to larger and larger bandwidths along the jωaj\omega_\text{a} axis in the ss-plane” [Smith07].

The mathematical relation between the analog frequencies ωa\omega_\text{a} and digital frequencies ωd\omega_\text{d} can be found by evaluating the transform along the analog frequency axis and the digital frequency axis. In other words, we need to insert s=jωas = j \omega_\text{a} and z=ejωdTz = e^{j\omega_\text{d} T} into Equation 1

jωa=2T1ejωdT1+ejωdT=2TejωdT/2(ejωdT/2ejωdT/2)ejωdT/2(e2jωdT/2+ejωdT/2)=2Tjsin(ωdT/2)cos(ωdT/2)=j2Ttan(ωdT/2).(12)j \omega_\text{a} = \frac{2}{T}\frac{1 - e^{-j\omega_\text{d}T}}{1 + e^{-j\omega_\text{d}T}} = \frac{2}{T}\frac{e^{-j\omega_\text{d}T/2}(e^{j\omega_\text{d}T/2} - e^{-j\omega_\text{d}T/2})}{e^{-j\omega_\text{d}T/2}(e^{2j\omega_\text{d}T/2} + e^{-j\omega_\text{d}T/2})} \newline = \frac{2}{T} \frac{j \sin(\omega_\text{d} T / 2)}{\cos(\omega_\text{d}T / 2)} = j \frac{2}{T} \tan(\omega_\text{d}T/2). \quad (12)

Therefore,

ωa=2Ttan(ωdT/2).(13)\omega_\text{a} = \frac{2}{T} \tan(\omega_\text{d}T/2). \quad (13)

From Equation 13, we can obtain the inverse relation

ωd=2Tarctan(ωaT/2).(14)\omega_\text{d} = \frac{2}{T} \text{arctan}(\omega_\text{a} T / 2). \quad (14)

The visualization of Equation 14 on Figure 5 shows what is frequency warping.

Figure 5. Plot showing how analog frequencies are mapped to digital frequencies by the bilinear transform.

Prewarping

In order to account for frequency warping of the bilinear transform in filter design, we need to perform the so-called prewarping. The prototype filter frequency scale must be prewarped so that after the bilinear transform the critical frequencies are in the correct places.

In other words, we want to have

Ha(jωa0)=Hd(ωd0)=H,(15)H_\text{a}(j \omega_\text{a0}) = H_\text{d}(\omega_\text{d0}) = H^*, \quad (15)

where we choose ωa0\omega_\text{a0} and ωd0\omega_\text{d0} (desired critical ferquencies of the prototype and digital filters respectively).

This goal can be easily achieved by simply scaling the analog frequency axis with a proper scalar before the bilinear transform.

In the context of parametric filters, our critical frequency will typically be the cutoff frequency of a filter, because that’s the frequency for which we want to have a specific transfer function value. But let’s derive the general scaling.

Prewarping Factor Derivation

Given critical frequencies ωa0\omega_\text{a0} and ωd0\omega_\text{d0}, we want to find a scaling factor KRK \in \mathbb{R} such that

Ha(jωa0)=Hd(ωd0)=H,(16)H_\text{a}(j \omega_\text{a0}) = H_\text{d}(\omega_\text{d0}) = H^*, \quad (16)

and

Hd(z)=Ha(Ks)zC,(17)H_\text{d}(z) = H_\text{a}(K s) \forall z \in \mathbb{C}, \quad (17)

where ss is given by Equation 1 (the bilinear transform).

From Equation 13, we know that without scaling by KK

ωa0=2Ttan(ωd0T/2).(18)\omega_\text{a0} = \frac{2}{T} \tan(\omega_\text{d0}T/2). \quad (18)

Here ωa0\omega_\text{a0} and ωd0\omega_\text{d0} are coupled; we cannot change them independently. To change this, we introduce the scaling KK (Equation 17)

ωa0=K2Ttan(ωd0T/2).(19)\omega_\text{a0} = K \frac{2}{T} \tan(\omega_\text{d0}T/2). \quad (19)

Solving for KK

K=ωa0tan(ωd0T/2)T2(20)K = \frac{\omega_\text{a0}}{\tan(\omega_\text{d0}T/2)} \frac{T}{2} \quad (20)

allows us to independently change ωa0\omega_\text{a0} and ωd0\omega_\text{d0}. \Box

So the bilinear transform combined with prewarping gives the following formula for ss substitution

s=ωa0tan(ωd0T/2)1z11+z1.(21)s = \frac{\omega_\text{a0}}{\tan(\omega_\text{d0}T/2)} \frac{1 - z^{-1}}{1 + z^{-1}}. \quad (21)

Cutoff Frequency from Prewarping

As you might recall from the analog prototype design tutorial, we design our analog filters to have the cutoff frequency equal to 1. Since thanks to prewarping, we can independently change ωa0\omega_\text{a0} and ωd0\omega_\text{d0}, we can simply set ωa0=1\omega_\text{a0} = 1 and ωd0=ωdc\omega_\text{d0} = \omega_\text{dc} (cutoff frequency of the digital filter). Thus, Equation 21 (bilinear transform + prewarping) becomes

s=1tan(ωd0T/2)1z11+z1.(22)s = \frac{1}{\tan(\omega_\text{d0} T/2)} \frac{1 - z^{-1}}{1 + z^{-1}}. \quad (22)

So prewarping actally lets us set the cutoff frequency of the filter we design 🙂

Example: Digitization of the Butterworth Low-Pass

Finally, let’s put to practice all that we’ve learned.

Let us digitize the analog prototype Butterworth low-pass of order 2 that we designed in the previous article.

It’s transfer function in the ss-domain was

Ha(s)=1s2+s2+1.(23)H_a(s) = \frac{1}{s^2 + s \sqrt{2} + 1}. \quad (23)

The cutoff frequency was set to 1.

To obtain a digital filter with the digital cutoff frequency equal to ωc\omega_\text{c} from Equation 23, we perform the bilinear transform with prewarping by inserting Equation 22 into Equation 23.

After very simple but tiresome derivations, we obtain

H2(z)=W2+2W2z1+W2z21+W2+W2+2(W21)z1+(W2W2+1)z2,(24)H_2(z) = \frac{W^2 + 2W^2 z^{-1} + W^2z^{-2}}{1 + W \sqrt{2} + W^2 + 2(W^2 - 1)z^{-1} + (W^2 - W\sqrt{2} + 1)z^{-2}}, \quad (24)

where W=tan(ωcT/2)W = \tan(\omega_\text{c} T / 2) and ωc\omega_\text{c} is the desired cutoff frequency of the digital filter in radians per second. The above result matches the one in [Zölzer08].

We thus obtained a stable IIR filter of order 2 with controllable cutoff frequency. Yay!

Efficiency

The filter in Equation 24 is efficient during processing but setting its cutoff frequency requires calculating the value of 5 coefficients (after normalization). That parameter-to-coefficient mapping is not efficient at all. Therefore, we need to turn to allpass-based filters… But that will be a topic of a different article. 😉

Summary

In summary, in this article we learned that the bilinear transform is a great method that allows us to digitize analog prototype filters.

We learned how to derive the bilinear transform by approximating a continuous-time integrator. We learned the properties of the bilinear transform such as stability and optimality preservation. We also counteracted the frequency warping with prewarping.

We learned that with prewarping we can set the cutoff frequency of analog prototypes to 1 because prewarping alters it anyway.

Finally, we discretized the second-order analog Butterworth low-pass with the bilinear transform and prewarping and obtained an IIR filter.

Thanks for reading! If you have any questions, don’t hesitate to ask them in the comments below 🙂

Bibiliography

Maximilian Schäfer, Music Processing - Synthesis lecture notes, Friedrich-Alexander-Universität Erlangen-Nürnberg, Winter Semester 2020/2021.

[OppenheimSchafer10] Alan V Oppenheim, Ronald W. Schafer, Discrete-Time Signal Processing, 3rd Edition, Pearson 2010.

[ParksBurrus87] T. W. Parks, C. S. Burrus, Digital Filter Design, John Wiley & Sons, Inc., 1987.

[Smith07] Julius O. Smith, Introduction to Digital Filters with Audio Applications, http://ccrma.stanford.edu/~jos/filters/, online book, 2007 edition, accessed November 26, 2021.

[Zölzer08] Zölzer Udo, Digital Audio Signal Processing, 2nd ed., Helmut Schmidt University, Hamburg, Germany, John Wiley & Sons Ltd, 2008.

Links above may be affiliate links. That means that I may earn a commission if you decide to make a purchase. This does not incur any cost for you. Thank you.