The Bilinear Transform and Frequency Warping

The analogue filter catalogue is a century of accumulated design: Butterworth’s maximally flat response, Chebyshev’s equal-ripple passband, elliptic filters with ripples in both bands. Each has closed-form pole locations and a known order for a given specification. Reusing that catalogue for discrete-time processing requires a map from the $s$-plane to the $z$-plane that preserves the properties that matter: stability, the shape of the frequency response, and causality. The bilinear transform is the standard choice because it satisfies the first and third exactly and the second only after a correction.

The map itself

The substitution is

$$s = \frac{2}{T}\cdot\frac{z-1}{z+1}$$

where $T$ is the sample period, or, more commonly in normalized form, $s \leftarrow (1-z^{-1})/(1+z^{-1})$ with the factor of two absorbed into the frequency scaling. The map sends the left half of the $s$-plane onto the interior of the unit disk, the imaginary axis onto the unit circle, and the point at infinity to $z = -1$.

The stability preservation is exact and easy to see: a pole at $s = -\alpha + j\beta$ with $\alpha > 0$ maps to a point inside the circle, for every $\beta$. Repeated application to a higher-order prototype preserves that property, so a stable analog filter always yields a stable digital filter. The converse holds only for poles that map back into the left half-plane, which they do for all finite $z$.

The map from the imaginary axis is where the difficulty lies. Setting $s = j\Omega$ and $z = e^{j\omega}$ in the substitution gives

$$\Omega = \frac{2}{T}\tan\left(\frac{\omega}{2}\right)$$

so the entire analog frequency axis, from zero to infinity, maps onto the unit circle once, from $\omega = 0$ to $\omega = \pi$. That is the map’s virtue: no aliasing, because the whole infinite axis is compressed into a finite arc rather than periodically repeated. It is also the source of the distortion: the compression is nonlinear.

The warping begins at the cutoff

Near $\omega = 0$ the tangent is approximately its argument, so $\Omega \approx \omega/T$ and the mapping is locally linear with the right slope. The correspondence degrades as frequency rises. At $\omega = \pi/2$, corresponding to a quarter of the sample rate, $\Omega = (2/T)\tan(\pi/4) = 2/T$, about 1.27 times the value the linear relation would predict. As $\omega$ approaches $\pi$, the analog frequency approaches infinity.

The practical effect is that a digital filter designed by transforming a prototype with a nominal cutoff of $\Omega_c$ does not have its $-3$ dB point at the frequency that the naive scaling $\omega_c = \Omega_c T$ suggests. If the analog cutoff is placed at $\Omega_c$ and the digital filter’s response is evaluated at $\omega$, the actual cutoff satisfies $\Omega_c = (2/T)\tan(\omega_c/2)$. Solving for the digital cutoff gives $\omega_c = 2\arctan(\Omega_c T/2)$.

For a cutoff at a tenth of the sample rate, $\omega_c = 0.2\pi \approx 0.628$, and the tangent correction is a few percent. At a quarter of the sample rate the correction is 27 percent, and above that it grows without bound. The discrepancy is not an error in the transform; it is the transform working as defined. The fix is to choose the analog prototype’s cutoff so that the desired digital frequency lands where it should, a step called prewarping.

import math

def prewarp(omega_digital):
    """Analog frequency that maps to the given digital frequency under s = tan(w/2)."""
    return math.tan(omega_digital / 2)

def digital_cutoff(Omega_analog):
    return 2 * math.atan(Omega_analog)

w = 0.2 * math.pi                 # desired digital cutoff, radians per sample
W = prewarp(w)                   # design the analog prototype at this cutoff
print(digital_cutoff(W) - w)     # ~0, the correction is exact

Prewarping is applied to the critical frequency, and only one frequency can be matched exactly. If the specification has both a passband edge and a stopband edge, the usual practice is to prewarp one of them, typically the passband edge, and accept the resulting shift at the other. The shift compresses the analog transition band, because the tangent is steeper at higher frequency, so a digital transition of fixed width corresponds to a narrower analog transition when it sits higher in frequency. The filter order derived from the analog specification is therefore conservative if the prewarp is done at the lower edge and slightly optimistic if done at the upper.

Order, ripple, and where the ripples land

The classical prototypes differ in how they distribute approximation error across the band, and the bilinear transform preserves that distribution in the warped coordinate. A Butterworth prototype has a magnitude response that is maximally flat at DC, monotonically decreasing, with the $-3$ dB point at the cutoff for every order. Its transition is the gentlest of the three families for a given order, and its phase response is the most nearly linear in the passband. A Chebyshev type I prototype concentrates its error as equal ripples in the passband and is monotonic in the stopband; for the same order and the same transition width it achieves a much lower stopband level than Butterworth. Chebyshev type II inverts that: flat passband, equal ripples in the stopband. Elliptic prototypes ripple in both bands and produce the sharpest transition for a given order, at the cost of the poorest phase response and the greatest sensitivity of the response to coefficient rounding.

For a given set of specifications, the required order is a formula, not a search. For Butterworth, the order needed to move from a passband attenuation $A_p$ at $\Omega_p$ to a stopband attenuation $A_s$ at $\Omega_s$ is

$$n \ge \frac{\log_{10}\left(\frac{10^{A_s/10}-1}{10^{A_p/10}-1}\right)}{2\log_{10}(\Omega_s/\Omega_p)}$$

with the analog frequencies prewarped first. A specification calling for 0.5 dB at the passband edge and 60 dB at twice that edge needs an order of 12 for Butterworth, 7 for Chebyshev with 0.5 dB ripple, and 5 for elliptic. The Chebyshev advantage over Butterworth grows with the ratio between the stopband and passband edges: at a ratio of two it saves five orders, and at a ratio of ten it saves one, because the hyperbolic cosine in its order formula grows logarithmically in the attenuation and only as the inverse cosine in the frequency ratio. Elliptic filters beat both at every ratio because their ripple is spent in both bands rather than one. Each step down in order reduces the arithmetic but moves the poles closer to the unit circle, which increases the sensitivity to coefficient quantization.

Why high-order sections are not implemented as one polynomial

A tenth-order filter expressed as a single transfer function with an order-10 numerator and denominator has coefficients spanning a wide dynamic range, and the roots of a polynomial are ill-conditioned with respect to its coefficients. Small rounding in the denominator coefficients moves the poles measurably, and a pole that crosses the unit circle turns a stable filter into an oscillator.

The standard implementation is a cascade of second-order sections, each a biquad with its own pole pair and zeros. Each section’s coefficients are computed from the corresponding complex-conjugate pole pair, so the conditioning problem is confined to one pair at a time. The sections are ordered to control the internal dynamic range: putting the sections with the smallest gain first and the largest last keeps intermediate values from over- or underflowing, and pairing each pole pair with the nearest zero pair minimizes the peak gain of each section.

import numpy as np
from scipy import signal

# Fourth-order Butterworth, cutoff at a fifth of the sample rate
sos = signal.butter(4, 0.2, btype='low', output='sos')
w, h = signal.sosfreqz(sos, worN=2048)
print(abs(h[0]))                       # DC gain, exactly 1 for Butterworth
print(20 * np.log10(abs(h).min()))     # deepest stopband point

The sos form stores six coefficients per section in a fixed layout, with the first three for the numerator and the last three for the denominator. A transposed direct form II structure is the usual per-section realization, because it needs only two state variables per section rather than four and has better roundoff behaviour when the coefficients are quantized. The direct form I structure, which keeps separate delay lines for input and output, is more forgiving of overflow but needs twice the state.

The limit at high frequencies

The transform is only usable when the digital cutoff is well below the Nyquist frequency. As $\omega_c$ approaches $\pi$, the prewarped analog frequency grows without limit and the required analog prototype order grows with it, so a specification asking for a sharp transition near half the sample rate yields a filter of impractical order with poles very close to the circle. In that regime the appropriate response is not a higher-order IIR design but a change of architecture: an oversampled front end, an FIR design with linear phase, or a multirate structure that splits the band before filtering.

The transform also produces a zero at $z = -1$ for every pole at infinity in the analog prototype. A lowpass prototype has all its zeros at infinity, so the digital filter has an equal number of zeros at $z = -1$, which forces the response to zero at exactly half the sample rate. That is a desirable property for a lowpass filter and an inconvenient one for a bandpass centered near Nyquist, where the enforced zero can cut into the passband. Multipoint transformations, which map a digital bandpass onto the analog lowpass before the bilinear step, move the zeros to the right places but add their own warping step.

For applications that need exactly linear phase, the IIR route is closed regardless of order, since a rational transfer function with poles inside the circle cannot have exactly linear phase except in trivial cases. The forward-backward two-pass technique described in Zero-Phase Filtering by Forward-Backward Passes recovers zero phase from any stable IIR filter by squaring its magnitude response, which is often an acceptable substitute for offline analysis. When the phase must be exactly linear in a single pass, the alternative is an FIR design whose symmetry guarantees it, and the length that costs is set by the window and transition rules in FIR Design by Windowing and the Gibbs Ripple.