Z-Transform and Digital Filter Design: Poles, Zeros, and Stability with Python

Introduction

When designing and analyzing digital filters, working directly with difference equations quickly becomes cumbersome. Just as the Laplace transform converts differential equations into algebraic equations in the continuous-time domain, the Z-transform converts difference equations into algebraic equations and makes the theoretical treatment of digital filters dramatically simpler.

With the Z-transform, the following can be understood systematically:

• Deriving the transfer function \(H(z)\) of a digital filter
• Intuitively grasping filter characteristics from the placement of poles and zeros in the complex plane
• Expressing system stability conditions mathematically
• Computing frequency response by evaluating \(H(z)\) on the unit circle

This article progresses from the definition and basic properties of the Z-transform through inverse transforms, pole-zero plots, transfer functions and stability, and the connection to frequency response. Python implementations using scipy are provided throughout.

Reading about Butterworth filter design and FIR vs IIR filter comparison beforehand will deepen your understanding of this article.

Z-Transform Definition and Region of Convergence

Bilateral Z-Transform

The bilateral Z-transform of a discrete-time signal \(x[n]\) (\(n \in \mathbb{Z}\) ) is defined as:

where \(z\) is a complex variable. Writing \(z = re^{j\omega}\) (\(r > 0\) ), equation \((1)\) expands to:

This has the same form as the Discrete-Time Fourier Transform (DTFT) of the signal \(x[n] r^{-n}\) . The Z-transform can therefore be viewed as a generalization of the DTFT to the entire complex plane.

Unilateral Z-Transform

When dealing with causal systems (\(x[n] = 0\) for \(n < 0\) ), the unilateral Z-transform is used:

Causal systems are standard in digital filter design; unless stated otherwise, we assume the unilateral Z-transform.

Region of Convergence (ROC)

The infinite series in equation \((1)\) does not necessarily converge for all values of \(z\) . The set of \(z\) values for which the Z-transform converges is called the Region of Convergence (ROC).

The ROC takes the form of an annular region \(r_1 < |z| < r_2\) centered at the origin. Key properties include:

• Finite-length signals: ROC is the entire complex plane except possibly \(z = 0\) or \(z = \infty\)
• Right-sided signals (non-zero only for \(n \geq 0\) ): ROC is \(|z| > r_1\) (outside some radius)
• Left-sided signals (non-zero only for \(n \leq 0\) ): ROC is \(|z| < r_2\) (inside some radius)
• The ROC never contains poles

As an example, compute the Z-transform of \(x[n] = a^n u[n]\) (\(u[n]\) is the unit step function):

This series converges when \(|az^{-1}| < 1\) , i.e., \(|z| > |a|\) . The ROC is the exterior region \(|z| > |a|\) .

Key Properties of the Z-Transform

Linearity

The ROC contains at least the intersection \(\text{ROC}_1 \cap \text{ROC}_2\) .

Time Shift

Proof: Substitute into the definition and set \(m = n - k\) :

This property is essential for converting difference equations into algebraic equations in the Z-domain. The factor \(z^{-1}\) acts as a “one-sample delay” operator.

Convolution Theorem

Convolution in the time domain becomes multiplication in the Z-domain. This means that filtering (convolution of input signal with impulse response) is expressed as multiplication of transfer functions.

Time Reversal

Z-Domain Differentiation (Multiplication Property)

Initial Value and Final Value Theorems

For causal signals, the following theorems hold.

Initial Value Theorem:

Final Value Theorem:

The final value theorem is useful for finding the steady-state output of a system using only Z-domain algebra, without simulating in the time domain. It applies only when all poles of \((z-1)X(z)\) lie strictly inside the unit circle.

Common Transform Pairs

Inverse Z-Transform

Partial Fraction Expansion Method

The most widely used method. Express \(X(z)\) as a rational function and decompose it into terms that match known transform pairs.

The standard procedure is to expand \(X(z)/z\) (or \(H(z)/z\) ) into partial fractions.

Example: Compute the inverse Z-transform of the following.

Expand \(X(z)/z\) in partial fractions:

Therefore:

Since the ROC is \(|z| > 0.8\) (right-sided signal):

Power Series Expansion (Long Division)

Expand \(X(z)\) as a power series in \(z^{-1}\) :

Performing long division of the numerator by the denominator yields the coefficients directly. This is convenient for FIR filters and numerical verification.

Example: Expand \(X(z) = \dfrac{1 + z^{-1}}{1 - 0.5z^{-1}}\) .

Long division gives \(x[0]=1,\; x[1]=1.5,\; x[2]=0.75,\; x[3]=0.375, \ldots\) , which yields \(x[n] = \delta[n] + 1.5(0.5)^{n-1}u[n-1]\) .

Residue Method (Cauchy Integral Formula)

where \(C\) is a counterclockwise closed contour inside the ROC and \(z_k\) are poles enclosed by \(C\) . This is important for theoretical analysis, but in practice the partial fraction method is more straightforward.

Poles, Zeros, and the Z-Plane

Definitions

For a rational function \(X(z) = N(z)/D(z)\) where \(N\) and \(D\) are polynomials:

• Values \(z_0\) where \(N(z_0) = 0\) : zeros
• Values \(z_p\) where \(D(z_p) = 0\) : poles

A rational function with numerator degree \(M\) and denominator degree \(N\) has \(M\) zeros and \(N\) poles (counting zeros/poles at infinity).

Pole-Zero Plot

Plotting poles as × and zeros as ○ on the complex plane (z-plane) produces a pole-zero plot. This is the single most important tool for visually grasping the characteristics of a digital filter.

The magnitude of the filter frequency response corresponds to \(|H(z)|\) as \(z\) sweeps along the unit circle \(z = e^{j\omega}\) :

• Zero near the unit circle: gain decreases (notch) near that angular frequency
• Pole near the unit circle: gain increases (peak) near that angular frequency
• Pole on or outside the unit circle: system is unstable

Stability Condition

An LTI (Linear Time-Invariant) discrete-time system is BIBO stable if and only if:

This has an intuitive explanation. The time-domain component corresponding to pole \(z_k = r_k e^{j\theta_k}\) is proportional to \(r_k^n e^{j\theta_k n}\) . If \(r_k < 1\) , this decays as \(n \to \infty\) ; if \(r_k > 1\) , it diverges.

Transfer Function and Difference Equations

Definition of the Transfer Function

The transfer function \(H(z)\) of a linear time-invariant system is defined as the ratio of the Z-transform of the output \(Y(z)\) to the Z-transform of the input \(X(z)\) :

By the convolution theorem, \(Y(z) = H(z) X(z)\) , so \(H(z)\) is the Z-transform of the impulse response \(h[n]\) .

From Difference Equation to Transfer Function

A general \(N\) -th order difference equation is:

Taking the Z-transform of both sides and applying the time-shift property \(\mathcal{Z}\{y[n-k]\} = z^{-k} Y(z)\) :

Rearranging:

This is the general form of the IIR filter transfer function.

FIR and IIR Filters in the Z-Domain

FIR filters (Finite Impulse Response) have no feedback (all \(a_k = 0\) ):

This is a polynomial with no poles (except at the origin). FIR filters are designed using zeros only.

IIR filters (Infinite Impulse Response) have feedback (\(a_k \neq 0\) ):

IIR filters have both poles and zeros. They can realize steep frequency responses with fewer coefficients, but require careful stability management. See FIR and IIR Filter Comparison for details.

Pole-Zero Form

Expressing the transfer function as a product of poles and zeros:

where \(q_k\) are zeros and \(p_k\) are poles. This representation corresponds directly to the pole-zero plot and offers an intuitive understanding of filter characteristics.

Stability Analysis

Stable, Unstable, and Marginally Stable Systems

Second-Order Sections and Numerical Stability

Implementing a high-order IIR filter directly can cause poles to migrate outside the unit circle due to finite-precision arithmetic errors. To prevent this, IIR filters are standardly decomposed into and implemented as Second-Order Sections (SOS):

scipy.signal.butter(..., output='sos') returns coefficients in this form.

Frequency Response from Poles and Zeros

Evaluating on the Unit Circle

The frequency response of a digital filter is obtained by evaluating the transfer function \(H(z)\) on the unit circle \(z = e^{j\omega}\) (\(\omega \in [-\pi, \pi]\) ):

This coincides with the DTFT. The Fast Fourier Transform (FFT) computes discrete samples of this.

Geometric Interpretation from Poles and Zeros

Substituting \(z = e^{j\omega}\) into equation \((22)\) :

Each factor can be interpreted as a “vector” in the complex plane:

• \(|e^{j\omega} - q_k|\) : distance from the point \(e^{j\omega}\) on the unit circle to zero \(q_k\)
• \(|e^{j\omega} - p_k|\) : distance from the point \(e^{j\omega}\) on the unit circle to pole \(p_k\)

Therefore:

From this geometric interpretation:

• A pole close to the unit circle at some angular frequency → denominator is small → gain increases (peak)
• A zero close to the unit circle at some angular frequency → numerator is small → gain decreases (notch)
• A zero on the unit circle at some angular frequency → \(|H| = 0\) (complete rejection)

For the theoretical treatment of phase characteristics, see also Phase and Group Delay Analysis.

Python Implementation

Importing Libraries

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, freqz, tf2zpk, zpk2tf, sosfilt

Pole-Zero Plot Function (zplane)

def zplane(b, a, ax=None, title="Pole-Zero Plot"):
    """
    Draw a pole-zero plot for a digital filter.

    Parameters
    ----------
    b : array_like
        Numerator coefficients (b[0] + b[1]*z^-1 + ...)
    a : array_like
        Denominator coefficients (a[0] + a[1]*z^-1 + ...)
    ax : matplotlib Axes, optional
        Target Axes. Creates a new figure if None.
    title : str
        Plot title

    Returns
    -------
    zeros : ndarray
        Array of zeros
    poles : ndarray
        Array of poles
    """
    if ax is None:
        fig, ax = plt.subplots(figsize=(6, 6))

    # Compute poles and zeros
    zeros, poles, gain = tf2zpk(b, a)

    # Draw unit circle
    theta = np.linspace(0, 2 * np.pi, 512)
    ax.plot(np.cos(theta), np.sin(theta), "k--", linewidth=0.8, alpha=0.5, label="Unit Circle")

    # Real and imaginary axes
    ax.axhline(0, color="k", linewidth=0.5, alpha=0.4)
    ax.axvline(0, color="k", linewidth=0.5, alpha=0.4)

    # Zeros (circles)
    ax.scatter(zeros.real, zeros.imag, s=80, marker="o", facecolors="none",
               edgecolors="blue", linewidths=1.5, zorder=5, label=f"Zeros ({len(zeros)})")

    # Poles (crosses)
    ax.scatter(poles.real, poles.imag, s=80, marker="x",
               color="red", linewidths=1.5, zorder=5, label=f"Poles ({len(poles)})")

    # Axis settings
    ax.set_xlabel("Real Part")
    ax.set_ylabel("Imaginary Part")
    ax.set_title(title)
    ax.set_aspect("equal")
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

    # Auto-set limits to accommodate unit circle
    lim = max(1.5, np.max(np.abs(poles)) * 1.2, np.max(np.abs(zeros)) * 1.2)
    ax.set_xlim(-lim, lim)
    ax.set_ylim(-lim, lim)

    return zeros, poles

Comparing Stable and Unstable Systems

# --- Stable system (all poles inside unit circle) ---
b_stable = np.array([1.0])
a_stable = np.array([1.0, -0.5, 0.3])   # poles: |p| < 1

# --- Unstable system (poles outside unit circle) ---
b_unstable = np.array([1.0])
a_unstable = np.array([1.0, -1.5, 0.9]) # poles: |p| >= 1

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
z_s, p_s = zplane(b_stable, a_stable, ax=axes[0],
                  title="Stable System\n(All Poles Inside Unit Circle)")
z_u, p_u = zplane(b_unstable, a_unstable, ax=axes[1],
                  title="Unstable System\n(Poles Outside Unit Circle)")

print(f"Stable system poles:   {np.roots(a_stable)}")
print(f"  Magnitudes: {np.abs(np.roots(a_stable))}")
print(f"Unstable system poles: {np.roots(a_unstable)}")
print(f"  Magnitudes: {np.abs(np.roots(a_unstable))}")

plt.tight_layout()
plt.show()

Sample output:

Stable system poles:   [0.25+0.479j  0.25-0.479j]
  Magnitudes: [0.541  0.541]
Unstable system poles: [0.75+0.661j  0.75-0.661j]
  Magnitudes: [1.0  1.0]  ← on unit circle (marginally stable)

Butterworth Filter Pole-Zero Analysis

# Butterworth low-pass filter (cutoff 100 Hz, sampling 1000 Hz)
fs = 1000   # Hz
fc = 100    # Hz
orders = [2, 4, 6]

fig, axes = plt.subplots(len(orders), 2, figsize=(14, 4 * len(orders)))

for i, N in enumerate(orders):
    # Design in SOS form
    sos = butter(N, fc, fs=fs, output='sos')
    # Convert to transfer function form
    b, a = butter(N, fc, fs=fs, output='ba')

    # Left column: pole-zero plot
    zplane(b, a, ax=axes[i, 0], title=f"Butterworth N={N} Pole-Zero Plot")

    # Right column: frequency response
    w, h = freqz(b, a, worN=4096, fs=fs)
    axes[i, 1].plot(w, 20 * np.log10(np.abs(h) + 1e-12), "b", linewidth=1.5)
    axes[i, 1].axvline(fc, color="gray", linestyle="--", alpha=0.7, label=f"fc={fc}Hz")
    axes[i, 1].axhline(-3, color="red", linestyle=":", alpha=0.7, label="-3 dB")
    axes[i, 1].set_xlim(0, fs / 2)
    axes[i, 1].set_ylim(-80, 5)
    axes[i, 1].set_xlabel("Frequency [Hz]")
    axes[i, 1].set_ylabel("Magnitude [dB]")
    axes[i, 1].set_title(f"Butterworth N={N} Frequency Response")
    axes[i, 1].legend()
    axes[i, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Frequency Response from Poles and Zeros (Geometric Visualization)

def freq_response_from_zpk(zeros, poles, gain, omega_range=None):
    """
    Compute magnitude and phase of the frequency response
    geometrically from poles and zeros.

    Parameters
    ----------
    zeros : array_like
        Array of zeros (complex)
    poles : array_like
        Array of poles (complex)
    gain : float
        Gain constant
    omega_range : array_like, optional
        Angular frequencies to evaluate (rad/sample)

    Returns
    -------
    omega : ndarray
        Angular frequencies
    H_mag : ndarray
        Magnitude response
    H_phase : ndarray
        Phase response [rad]
    """
    if omega_range is None:
        omega_range = np.linspace(0, np.pi, 1024)

    H_mag = np.ones(len(omega_range)) * np.abs(gain)
    H_phase = np.zeros(len(omega_range))

    for idx in range(len(omega_range)):
        z = np.exp(1j * omega_range[idx])

        # Product of distances from zeros
        for q in zeros:
            H_mag[idx] *= np.abs(z - q)
            H_phase[idx] += np.angle(z - q)

        # Product of distances from poles
        for p in poles:
            H_mag[idx] /= np.abs(z - p)
            H_phase[idx] -= np.angle(z - p)

    return omega_range, H_mag, H_phase


# Example: 2nd-order Butterworth LPF via geometric computation
b, a = butter(2, 0.2)   # digital cutoff at 0.2π (normalized)
zeros, poles, gain = tf2zpk(b, a)

omega, H_mag, H_phase = freq_response_from_zpk(zeros, poles, gain)

# Compare with scipy.signal.freqz
w, h = freqz(b, a, worN=1024)

fig, axes = plt.subplots(2, 1, figsize=(10, 7))

axes[0].plot(omega / np.pi, 20 * np.log10(H_mag + 1e-12), "b--",
             linewidth=2, label="Geometric (from z-plane)")
axes[0].plot(w / np.pi, 20 * np.log10(np.abs(h) + 1e-12), "r",
             linewidth=1, alpha=0.7, label="scipy.signal.freqz")
axes[0].set_xlabel("Normalized Frequency [× π rad/sample]")
axes[0].set_ylabel("Magnitude [dB]")
axes[0].set_title("Frequency Response: Geometric vs freqz")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(-60, 5)

axes[1].plot(omega / np.pi, np.degrees(np.unwrap(H_phase)), "b--",
             linewidth=2, label="Geometric")
axes[1].plot(w / np.pi, np.degrees(np.unwrap(np.angle(h))), "r",
             linewidth=1, alpha=0.7, label="scipy.signal.freqz")
axes[1].set_xlabel("Normalized Frequency [× π rad/sample]")
axes[1].set_ylabel("Phase [degrees]")
axes[1].set_title("Phase Response")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Inverse Z-Transform via Partial Fractions (Numerical Verification)

from scipy.signal import residuez

# H(z) = (1 + 0.5z^-1) / (1 - 0.9z^-1 + 0.2z^-2)
b = np.array([1.0, 0.5])
a = np.array([1.0, -0.9, 0.2])

# residuez returns H(z) = sum(r_k / (1 - p_k z^-1)) + k
r, p, k = residuez(b, a)

print("Partial fraction expansion: H(z) = sum(r_k / (1 - p_k * z^-1)) + k(z)")
for i, (ri, pi) in enumerate(zip(r, p)):
    print(f"  r[{i}] = {ri:.4f},  p[{i}] = {pi:.4f}  (|p| = {abs(pi):.4f})")
print(f"  Direct term k = {k}")

# Inverse Z-transform (analytical impulse response)
N = 20
n = np.arange(N)
h_analytical = np.zeros(N, dtype=complex)
for ri, pi in zip(r, p):
    h_analytical += ri * pi**n
h_analytical = h_analytical.real

# Compare with numerical computation via scipy.signal.lfilter
from scipy.signal import lfilter, unit_impulse
impulse = unit_impulse(N)
h_numerical = lfilter(b, a, impulse)

fig, ax = plt.subplots(figsize=(10, 4))
ax.stem(n, h_numerical, linefmt="b-", markerfmt="bo", basefmt="k-",
        label="lfilter (numerical)")
ax.stem(n, h_analytical, linefmt="r--", markerfmt="rx", basefmt="k-",
        label="Partial fraction (analytical)")
ax.set_xlabel("Sample n")
ax.set_ylabel("Amplitude")
ax.set_title("Impulse Response: Partial Fraction vs Direct Computation")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print(f"\nMax error: {np.max(np.abs(h_numerical - h_analytical)):.2e}")

Interconversion: Transfer Function, SOS, and ZPK Forms

from scipy.signal import tf2sos, sos2tf, zpk2sos

# 4th-order Butterworth in each representation
b, a = butter(4, 0.3)

# BA → ZPK
zeros, poles, gain = tf2zpk(b, a)

# ZPK → BA (round-trip)
b_r, a_r = zpk2tf(zeros, poles, gain)

# BA → SOS
sos = tf2sos(b, a)

# SOS → BA
b_sos, a_sos = sos2tf(sos)

print("=== 4th-Order Butterworth Filter: Representation Formats ===")
print(f"\nBA form (b): {b}")
print(f"BA form (a): {a}")
print(f"\nPoles: {poles}")
print(f"Zeros: {zeros}")
print(f"Gain:  {gain:.6f}")
print(f"\nSOS form:\n{sos}")

# Verify stability (all pole magnitudes < 1)
print(f"\nAll pole magnitudes (must all be < 1 for stability):")
for i, pol in enumerate(poles):
    print(f"  p[{i}] = {pol:.4f},  |p| = {abs(pol):.4f} {'OK' if abs(pol) < 1 else 'UNSTABLE'}")

Summary

This article provided a systematic treatment of the Z-transform and its applications.

The Z-transform is the cornerstone of digital signal processing theory. Once you can quickly assess filter characteristics just by examining the pole-zero placement, your understanding of filter design will deepen considerably.

Related Articles

• Fast Fourier Transform (FFT): How It Works and Python Implementation — Foundation of frequency analysis. Restricting the Z-transform to the unit circle gives the DTFT; its discrete samples give the DFT/FFT.
• Butterworth Filter Design: Principles and Python Implementation — A leading example of IIR filter design. The pole-zero plot and transfer function knowledge from this article applies directly.
• FIR and IIR Filter Comparison — Explains the Z-domain difference between FIR (all-zero filter) and IIR (has poles).
• Phase and Group Delay Analysis — Derives group delay from the phase of \(H(e^{j\omega})\) and analyzes the time-domain distortion of filters.

References

• Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
• Proakis, J. G., & Manolakis, D. G. (2006). Digital Signal Processing: Principles, Algorithms, and Applications (4th ed.). Prentice Hall.
• Antoniou, A. (2005). Digital Signal Processing: Signals, Systems, and Filters. McGraw-Hill.
• SciPy scipy.signal documentation