1. Why Another Frequency Domain?

The z-transform is the natural domain for discrete-time systems — poles and zeros live on or inside the unit circle, and stability has a clean geometric interpretation. But designing filters is easiest in the continuous-time s-domain, where decades of classical theory (Butterworth, Chebyshev, elliptic) give analytic prototype responses.

The challenge is bridging these two worlds. Impulse-invariant transformation preserves time-domain shape but suffers from aliasing in the frequency domain. The bilinear transformation eliminates aliasing entirely by establishing a one-to-one, conformal mapping between the entire s-plane (or an equivalent pseudo-analog plane) and the z-plane. The price paid is a non-linear frequency compression known as frequency warping — the central theme of this article.

📋 Scope of this article We derive the mapping \(W = \tfrac{2}{T}\tfrac{z-1}{z+1}\) from first principles, characterise the image of every canonical region of the z-plane, expose the warping relation \(\Omega = \tfrac{2}{T}\tan\!\left(\tfrac{\omega}{2}\right)\), and show how pre-warping restores exact frequency correspondence at critical points.

2. The Bilinear Transformation Defined

Let \(z \in \mathbb{C}\) be the complex z-plane variable and \(T > 0\) the sampling period. The bilinear (Tustin) transformation is the Möbius map:

Forward Transform (z → W) $$W = \frac{2}{T} \cdot \frac{z - 1}{z + 1}$$

Its inverse — mapping from the w-plane back to the z-plane — is obtained by solving for \(z\):

Inverse Transform (W → z) $$z = \frac{1 + \tfrac{T}{2}W}{1 - \tfrac{T}{2}W} = \frac{2 + TW}{2 - TW}$$

Both directions are rational functions of degree one (Möbius transformations), so the mapping is conformal (angle-preserving) at every point except the poles \(z = -1\) and \(W = 2/T\).

2.1 Notation Conventions

Throughout this article we write \(z = \sigma_z + j\omega_z\) (real and imaginary parts in the z-plane — not to be confused with the digital frequency \(\omega\)), and \(W = \Sigma + j\Omega\) (real and imaginary parts in the w-plane, directly analogous to the Laplace variable \(s = \sigma + j\omega\)). The digital radian frequency is \(\omega \in [-\pi, \pi)\), occurring on the unit circle \(z = e^{j\omega}\).

3. Geometric Mapping: Where Does Each Region Go?

The power of understanding the bilinear transform lies in seeing exactly how it maps the canonical regions of the z-plane. Write \(z = x + jy\) and compute \(W = \Sigma + j\Omega\) by direct substitution.

3.1 Derivation of Real and Imaginary Parts

Step 1 — Substitute
$$W = \frac{2}{T} \cdot \frac{(x-1) + jy}{(x+1) + jy}$$
Step 2 — Rationalise (multiply numerator & denominator by conjugate of denominator)
$$W = \frac{2}{T} \cdot \frac{[(x-1)+jy]\,[(x+1)-jy]}{(x+1)^2 + y^2}$$
Step 3 — Expand numerator
Numerator \(= (x^2 - 1 + y^2) + j\,2y\)
Step 4 — Separate real and imaginary parts (master formulae)
$$\Sigma = \frac{2}{T} \cdot \frac{x^2 + y^2 - 1}{(x+1)^2 + y^2}, \qquad \Omega = \frac{2}{T} \cdot \frac{2y}{(x+1)^2 + y^2}$$

These two expressions are the master formulae. From them every geometric mapping result follows immediately.

3.2 The Unit Circle Maps to the Imaginary Axis

On the unit circle, \(|z| = 1\), i.e.\ \(x^2 + y^2 = 1\). Substituting:

$$\Sigma\big|_{|z|=1} = \frac{2}{T} \cdot \frac{1 - 1}{(x+1)^2 + y^2} = 0$$

So the entire unit circle maps to the imaginary axis \(\Sigma = 0\). The imaginary axis of the w-plane is the analog frequency axis — this is exactly what we want for filter design.

3.3 Inside the Unit Circle Maps to the Left Half-Plane

For \(|z| < 1\) we have \(x^2 + y^2 < 1\), so the numerator of \(\Sigma\) is negative and the denominator is always positive, giving \(\Sigma < 0\). The open unit disk maps bijectively onto the open left half of the w-plane — the stability region of continuous-time systems.

3.4 Outside the Unit Circle Maps to the Right Half-Plane

Symmetrically, \(|z| > 1 \Rightarrow x^2 + y^2 > 1 \Rightarrow \Sigma > 0\). The unstable z-plane region maps to the unstable right half w-plane.

Figure 1: The bilinear transform maps the open unit disk bijectively onto the open left half w-plane, and the unit circle onto the imaginary axis. The pole of W at z = −1 maps to W = ∞.

Table 1 — Key point and region mappings under W = (2/T)(z−1)/(z+1)
z-Plane Region Condition w-Plane Image Physical Meaning
Unit circle \(|z| = 1\) Imaginary axis (\(\Sigma = 0\)) Analog frequency axis
Open unit disk \(|z| < 1\) Open left half-plane (\(\Sigma < 0\)) Stable continuous-time region
Exterior of unit disk \(|z| > 1\) Open right half-plane (\(\Sigma > 0\)) Unstable continuous-time region
\(z = +1\) (DC) — \(W = 0\) DC ↔ DC preserved
\(z = -1\) (Nyquist) — \(W \to \infty\) Nyquist ↔ infinite analog freq.
\(z = j\) (\(\omega = \pi/2\)) — \(W = j\tfrac{2}{T}\) Quarter-sample freq. ↔ \(\Omega = 2/T\)

4. The Frequency Warping Relation

The most practically important consequence of the bilinear transform is the non-linear relationship between digital frequency \(\omega\) and pseudo-analog frequency \(\Omega\). We derive it by restricting \(z = e^{j\omega}\) on the unit circle.

Substitute \(z = e^{j\omega}\)
$$W = \frac{2}{T} \cdot \frac{e^{j\omega} - 1}{e^{j\omega} + 1}$$
Factor \(e^{j\omega/2}\) from numerator and denominator
$$W = \frac{2}{T} \cdot \frac{e^{j\omega/2}(e^{j\omega/2} - e^{-j\omega/2})}{e^{j\omega/2}(e^{j\omega/2} + e^{-j\omega/2})} = \frac{2}{T} \cdot \frac{2j\sin(\omega/2)}{2\cos(\omega/2)}$$
Result: W is purely imaginary on the unit circle
$$W = j\,\frac{2}{T}\tan\!\left(\frac{\omega}{2}\right) \quad \Longrightarrow \quad \boxed{\,\Omega = \frac{2}{T}\tan\!\left(\frac{\omega}{2}\right)\,}$$
💡 Key Insight — The Warping Equation The mapping between analog frequency \(\Omega\) and digital frequency \(\omega\) is a tangent function — highly non-linear except near the origin. This non-linearity is what compresses the entire positive real axis \((0,+\infty)\) into \((0,\pi)\), eliminating aliasing entirely.

4.1 Inverse Warping

Given a desired digital cut-off \(\omega_c\), the corresponding pre-warped analog frequency is:

$$\Omega_c = \frac{2}{T}\tan\!\left(\frac{\omega_c}{2}\right)$$

This is used to pre-warp the analog prototype specification before applying the inverse bilinear transform. After transformation, the digital filter will have its transition band exactly at \(\omega_c\).

Figure 2: Frequency warping curve Ω = (2/T)tan(ω/2) vs. digital frequency ω. The dashed line is the ideal linear relationship. Near ω = 0 the curves coincide; near ω = π the tangent diverges, compressing the entire positive real axis onto (0, π) — this is how aliasing is eliminated.

4.2 Key Properties of the Warping Function

  • Near DC (\(\omega \approx 0\)): \(\tan(\omega/2) \approx \omega/2\), so \(\Omega \approx \omega / T = \omega \cdot f_s / (2\pi)\). The mapping is approximately linear at low frequencies.
  • At the Nyquist frequency (\(\omega = \pi\)): \(\tan(\pi/2) \to \infty\). The entire semi-infinite positive analog frequency axis \((0, +\infty)\) is compressed into \((0, \pi)\). This is the mechanism by which aliasing is eliminated.
  • Monotonicity: \(\Omega(\omega)\) is strictly increasing on \([0, \pi)\), so magnitude response shapes (monotone, equiripple) are preserved — only frequency axis labels change.
  • Odd symmetry: \(\Omega(-\omega) = -\Omega(\omega)\), consistent with the real-valued filter requirement.

5. Conformality and the Möbius Structure

The bilinear transformation belongs to the family of Möbius (linear fractional) transformations — maps of the form \(f(z) = (az+b)/(cz+d)\) with \(ad - bc \neq 0\). Every Möbius transformation is a conformal automorphism of the Riemann sphere; it maps circles and lines to circles and lines.

For \(W = \tfrac{2}{T}\cdot\tfrac{z-1}{z+1}\) we identify \(a = 2/T,\ b = -2/T,\ c = 1,\ d = 1\), giving determinant \(ad - bc = 4/T \neq 0\). The specific circle mapped to a line here is the unit circle (a circle passing through no singularity of the transform), mapped to the imaginary axis — a line, which in the projective sense is a circle through \(\infty\).

The conformal property has a concrete engineering consequence: it guarantees that the shape of poles and zero loci in the z-plane is preserved under the mapping. An ellipse of poles in the s-plane becomes a (warped) ellipse in the z-plane, ensuring that Butterworth, Chebyshev, and elliptic filter pole patterns map cleanly.

6. Transfer Function Conversion

Given a continuous-time prototype transfer function \(H_a(s)\), the digital IIR filter is obtained by substituting \(s \leftarrow W\):

$$H(z) = H_a(W)\Big|_{W = \tfrac{2}{T}\cdot\tfrac{z-1}{z+1}} = H_a\!\left(\frac{2}{T}\cdot\frac{z-1}{z+1}\right)$$

Because \(W\) is a rational function of \(z\) of degree one, and \(H_a(s)\) is rational in \(s\), \(H(z)\) is rational in \(z\) of the same order as \(H_a\). No order inflation occurs.

6.1 Worked Example: First-Order Butterworth LPF

The normalised (cut-off at \(\Omega_c = 1\) rad/s) first-order Butterworth prototype is:

Analog Prototype$$H_a(s) = \frac{1}{s + 1}$$

Suppose we want a digital LPF with cut-off at \(\omega_c = \pi/4\) rad/sample and sampling period \(T = 1\) s. Step 1 — pre-warp:

$$\Omega_c = \frac{2}{1}\tan\!\left(\frac{\pi/4}{2}\right) = 2\tan\!\left(\frac{\pi}{8}\right) \approx 0.8284 \text{ rad/s}$$

Step 2 — frequency-scale the prototype to cut-off \(\Omega_c\):

$$H_a'(s) = H_a\!\left(\frac{s}{\Omega_c}\right) = \frac{1}{s/\Omega_c + 1} = \frac{\Omega_c}{s + \Omega_c}$$

Step 3 — apply the bilinear substitution \(s = 2(z-1)/(z+1)\) (with \(T=1\)):

$$H(z) = \frac{\Omega_c}{\dfrac{2(z-1)}{z+1} + \Omega_c} = \frac{\Omega_c(z+1)}{(2+\Omega_c)z + (\Omega_c - 2)}$$

Dividing numerator and denominator by \((2 + \Omega_c)\) and writing \(\alpha = \Omega_c/(2+\Omega_c)\):

$$H(z) = \alpha\cdot\frac{1 + z^{-1}}{1 - (1 - 2\alpha)z^{-1}}$$

With \(\Omega_c \approx 0.8284\): \(\alpha \approx 0.2929\), giving pole at \(z \approx 0.4142\) — well inside the unit circle, confirming stability. The digital cut-off is exactly at \(\omega_c = \pi/4\) by construction of the pre-warp step.

7. Pre-Warping: Correcting for Frequency Compression

Frequency warping would corrupt the cut-off and transition band specifications if not corrected. The remedy is pre-warping: before designing the analog prototype, map every critical frequency specification through the warping relation.

  1. Identify all critical digital frequencies \(\omega_1, \omega_2, \ldots\) (pass-band edge, stop-band edge, etc.).
  2. Pre-warp each: \(\Omega_k = \tfrac{2}{T}\tan(\omega_k/2)\).
  3. Design the analog prototype \(H_a(s)\) to meet the specification at \(\Omega_1, \Omega_2, \ldots\)
  4. Apply the bilinear substitution \(s = \tfrac{2}{T}\tfrac{z-1}{z+1}\) to obtain \(H(z)\).
  5. The digital filter \(H(z)\) meets the original specification at \(\omega_1, \omega_2, \ldots\) exactly, because the pre-warping and warping errors cancel.
⚠️ Pre-warping precision Pre-warping guarantees exact frequency correspondence only at the pre-warped critical frequencies. Between critical frequencies, the shape of the response is warped by \(\tan(\omega/2)\). For narrow-band filters (small \(\omega_c\)) this distortion is negligible; for wide-band filters it must be assessed case by case.

8. Comparison with Impulse-Invariant Transformation

Table 2 — Bilinear Transform vs. Impulse-Invariant Method
Property Bilinear Transform Impulse-Invariant
Aliasing None (compression prevents aliasing) Present — spectra overlap at \(\omega = \pm\pi\)
Frequency axis mapping Non-linear (\(\tan\) warping) Linear (\(\Omega = \omega/T\))
Filter order preserved Yes Yes
Impulse response match No Yes (by definition)
Suitable for HPF/BSF Yes No (severe aliasing)
Pre-warping required Yes (for critical frequencies) No
Preferred for All standard IIR filter types LPF approximations, matched-z

9. Python Implementation

SciPy provides bilinear and bilinear_zpk for direct conversion, and iirdesign/butter with the fs argument automates the entire pipeline including pre-warping. 

import numpy as np
from scipy.signal import butter, freqz, bilinear_zpk

# ── Parameters ──────────────────────────────────────────────
fs   = 8000          # sampling rate [Hz]
T    = 1 / fs        # sampling period [s]
fc   = 1000          # desired digital cut-off [Hz]
wc   = 2 * np.pi * fc / fs   # digital cut-off [rad/sample]

# ── Step 1: pre-warp ────────────────────────────────────────
Omega_c = (2 / T) * np.tan(wc / 2)
print(f"Pre-warped analog cut-off: {Omega_c/(2*np.pi):.2f} Hz")

# ── Step 2: design analog Butterworth prototype ──────────────
# scipy.signal.butter with analog=True uses Omega_c directly
from scipy.signal import buttap, zpkbilinear
z_a, p_a, k_a = buttap(4)           # 4th-order normalised prototype

# Scale prototype poles to Omega_c
p_a_scaled = p_a * Omega_c
k_a_scaled = k_a * Omega_c**4       # gain scales as Omega_c^N

# ── Step 3: apply bilinear transform ────────────────────────
z_d, p_d, k_d = bilinear_zpk(z_a, p_a_scaled, k_a_scaled, fs=fs)
print("Digital poles:", np.round(p_d, 4))

# ── Step 4: frequency response ───────────────────────────────
from scipy.signal import zpk2tf
b, a = zpk2tf(z_d, p_d, k_d)
w, H = freqz(b, a, worN=4096, fs=fs)

# ── Verify: magnitude at fc should be ≈ -3 dB ───────────────
idx = np.argmin(np.abs(w - fc))
print(f"|H(e^jω)| at {fc} Hz = {20*np.log10(np.abs(H[idx])):.3f} dB")

The bilinear_zpk function internally computes \(z_d = (1 + z_{a,k}/(2f_s)) / (1 - z_{a,k}/(2f_s))\) for each analog zero and the analogous formula for poles, corresponding exactly to our inverse mapping \(z = (2 + TW)/(2 - TW)\).

10. Limitations and Practical Considerations

  • Phase non-linearity: The bilinear transform does not preserve the group delay of the analog prototype. If linear phase is required, FIR designs are preferred; the bilinear transform is inappropriate.
  • Wide-band stop-bands: The \(\tan\) compression stretches the stop-band in digital frequency space, so a stop-band specified over a wide analog range appears narrower than expected. This is usually beneficial (actual stop-band attenuation improves), but must be verified.
  • High-order filters and numerical conditioning: Converting high-order (\(N \geq 8\)) filters via direct transfer-function polynomial form is numerically ill-conditioned. Always prefer second-order section (SOS) cascades, obtainable via butter(..., output='sos') in SciPy.
  • Sampling rate selection: The approximation \(\Omega \approx \omega/T\) holds for \(\omega \ll \pi\). For accurate broad-band designs, choose \(f_s\) at least 5–10× above the highest critical frequency.

References

  • A. V. Oppenheim & R. W. Schafer — Discrete-Time Signal Processing, 3rd ed., Pearson, 2010. §7.1–7.3.
  • J. G. Proakis & D. G. Manolakis — Digital Signal Processing: Principles, Algorithms, and Applications, 4th ed., Pearson, 2007. §8.3.
  • L. R. Rabiner & B. Gold — Theory and Application of Digital Signal Processing, Prentice Hall, 1975.
  • S. K. Mitra — Digital Signal Processing: A Computer-Based Approach, 4th ed., McGraw-Hill, 2011. §9.5.
  • W. S. Burdic — Radar Signal Analysis, Prentice Hall, 1968. (Tustin's original formulation in the context of control systems.)