Phase Spectrum and Group Delay: Theory and Python Implementation

A complete guide to phase spectrum and group delay in digital filters: mathematical definitions, linear phase FIR filters, nonlinear phase distortion in IIR filters, and all-pass phase equalization with Python.

Introduction

When performing spectral analysis with the Fast Fourier Transform (FFT), the result is commonly displayed as a magnitude spectrum showing how much power each frequency component carries. However, the Fourier transform output is complex-valued, and alongside the magnitude it encodes a phase spectrum that is equally important yet frequently overlooked.

Why does phase matter in practice?

Consider applying a Butterworth filter to a speech or biomedical signal. Even when the magnitude response is nearly ideal, the filtered waveform can exhibit visible distortion — transients spread out, sharp edges become blurred, and the temporal relationships between frequency components are altered. This distortion arises entirely from the filter’s nonlinear phase response: different frequency components are delayed by different amounts as they pass through the filter. The quantity that captures this frequency-dependent delay is the group delay.

This article provides a systematic treatment of phase spectrum and group delay: their mathematical definitions, the physical interpretation of group delay as signal travel time, the linear phase property of FIR filters, the nonlinear phase distortion inherent to IIR filters, and the use of all-pass filters for phase equalization — all illustrated with Python implementations using scipy.signal.

Phase Spectrum: Definition and Properties

Complex Frequency Response and Phase

The frequency response of a discrete-time linear filter evaluated on the unit circle is:

\[H(e^{j\omega}) = |H(e^{j\omega})| \cdot e^{j\phi(\omega)} \tag{1}\]

where \(\omega \in [0, \pi]\) is the normalized angular frequency (\(\omega = \pi\) corresponds to the Nyquist frequency), and:

  • \(|H(e^{j\omega})|\) — magnitude response: the gain applied to each frequency component
  • \(\phi(\omega) = \angle H(e^{j\omega})\) — phase response: the phase shift imparted to each frequency component

The phase response is computed as:

\[\phi(\omega) = \arctan\left(\frac{\mathrm{Im}[H(e^{j\omega})]}{\mathrm{Re}[H(e^{j\omega})]}\right) \tag{2}\]

Wrapped vs. Unwrapped Phase

Because \(\arctan\) has range \((-\pi, \pi]\) , the directly computed phase is wrapped: it exhibits artificial \(\pm\pi\) discontinuities even when the true underlying phase varies continuously. For instance, a continuously decreasing phase may appear as \(-3.10 \to -3.14 \to 3.12 \to 3.10\) , with a \(2\pi\) jump at the boundary.

The unwrapped phase removes these discontinuities by adding the appropriate integer multiples of \(2\pi\) :

\[\phi_{\mathrm{unwrap}}(\omega) = \phi(\omega) + 2\pi \, k(\omega), \quad k(\omega) \in \mathbb{Z} \tag{3}\]

In Python, np.unwrap() performs this operation automatically. Unwrapped phase is essential for meaningful group delay computation, since numerical differentiation of the wrapped phase produces large spurious spikes at each \(2\pi\) discontinuity.

Linear Phase and Nonlinear Phase

The most significant distinction in filter phase behavior is between linear phase and nonlinear phase.

A filter has linear phase if:

\[\phi(\omega) = -\omega \cdot \tau_0 \tag{4}\]

where \(\tau_0\) is a constant called the phase delay. In this case, every frequency component is shifted by a phase proportional to its frequency — which is mathematically equivalent to a pure time delay of \(\tau_0\) samples. The output waveform is an exact replica of the input, shifted in time.

A filter has nonlinear phase when \(\phi(\omega)\) deviates from a straight line. In this case, the delay experienced by each frequency component differs, and a waveform composed of multiple frequencies will emerge distorted.

Group Delay: Definition and Physical Interpretation

Mathematical Definition

The group delay is defined as the negative derivative of the phase response with respect to angular frequency:

\[\tau_g(\omega) = -\frac{d\phi(\omega)}{d\omega} \tag{5}\]

Group delay has units of samples (or seconds when multiplied by the sampling period \(T_s = 1/f_s\) ). When \(\tau_g(\omega) = \tau_0\) (constant), the filter has linear phase. When \(\tau_g(\omega)\) varies with \(\omega\) , the filter has nonlinear phase and introduces group delay distortion.

Substituting the linear phase condition \(\phi(\omega) = -\omega \tau_0\) into equation \((5)\) confirms:

\[\tau_g(\omega) = -\frac{d(-\omega \tau_0)}{d\omega} = \tau_0 \tag{6}\]

Physical Interpretation: Envelope Delay

The physical meaning of group delay becomes clear by examining a narrowband signal — a slowly modulated carrier:

\[x(t) = a(t)\cos(\omega_0 t)\]

where \(a(t)\) is the slowly varying envelope and \(\omega_0\) is the carrier frequency. After passing through a filter \(H\) , the output is approximately:

\[y(t) \approx |H(e^{j\omega_0})| \cdot a(t - \tau_g(\omega_0)) \cdot \cos(\omega_0 t - \phi(\omega_0)) \tag{7}\]

The envelope \(a(t)\) is delayed by the group delay \(\tau_g(\omega_0)\) , while the carrier is phase-shifted by \(\phi(\omega_0)\) . Thus, the group delay quantifies how many samples (or seconds) a signal packet centered at frequency \(\omega_0\) is delayed by the filter.

If group delay is frequency-independent (linear phase), all signal packets arrive at the output simultaneously and the waveform is preserved. If group delay varies with frequency (group delay distortion), different frequency bands arrive at different times, causing the output waveform to be stretched, smeared, or distorted relative to the input.

Consequences of Group Delay Distortion

Group delay distortion causes problems in many practical applications:

  • Waveform distortion: A square wave contains harmonics at odd multiples of the fundamental. If high harmonics are delayed more than the fundamental (or vice versa), the reconstructed waveform loses its sharp edges.
  • Biomedical signal processing: In ECG analysis, the shape of the QRS complex is diagnostically significant. Phase distortion can alter the morphology of these transient waveforms, potentially affecting automated detection algorithms.
  • Audio systems: Group delay variations across the audio band cause audible coloration — transients lose crispness, and the perceived spatial localization of sounds may shift.
  • Digital communications: In bandlimited channels, frequency-dependent group delay causes inter-symbol interference (ISI), spreading each transmitted symbol into adjacent time slots.
  • Feedback control: In closed-loop systems, group delay constitutes a frequency-dependent phase lag that erodes phase margin and can destabilize the loop.

Linear Phase FIR Filters

Why Symmetric Coefficients Guarantee Linear Phase

An \(N\) -tap FIR filter with impulse response \(h[n]\) has transfer function:

\[H(z) = \sum_{n=0}^{N-1} h[n] z^{-n} \tag{8}\]

If the coefficients satisfy the even symmetry condition:

\[h[n] = h[N-1-n], \quad n = 0, 1, \ldots, N-1 \tag{9}\]

then the frequency response can be factored as:

\[H(e^{j\omega}) = e^{-j\omega(N-1)/2} \cdot A(\omega) \tag{10}\]

where \(A(\omega)\) is a real-valued amplitude function. This factorization shows that the phase response is:

\[\phi(\omega) = -\frac{N-1}{2}\omega + \begin{cases} 0 & \text{if } A(\omega) > 0 \\ \pi & \text{if } A(\omega) < 0 \end{cases} \tag{11}\]

In the passband where \(A(\omega) > 0\) , the phase is exactly linear: \(\phi(\omega) = -\frac{N-1}{2}\omega\) . The group delay is therefore constant:

\[\tau_g = \frac{N-1}{2} \quad \text{[samples]} \tag{12}\]

A key insight is that the group delay equals half the filter length minus one. A 65-tap FIR filter introduces a fixed delay of 32 samples at all frequencies within the passband — a pure time shift, with no waveform distortion.

The Four Types of Linear Phase FIR Filters

Linear phase FIR filters are classified into four types based on coefficient symmetry and filter length parity:

TypeSymmetryLengthApplicable filter typesConstraint
IEven: \(h[n] = h[N-1-n]\)OddLowpass, highpass, bandpassNone (most general)
IIEven: \(h[n] = h[N-1-n]\)EvenLowpass\(H(\pi) = 0\) (cannot be highpass)
IIIOdd: \(h[n] = -h[N-1-n]\)OddDifferentiator, Hilbert\(H(0) = H(\pi) = 0\)
IVOdd: \(h[n] = -h[N-1-n]\)EvenHighpass, differentiator\(H(0) = 0\)

Type I is the most commonly used for general-purpose filtering. scipy.signal.firwin produces Type I or II filters depending on whether the specified number of taps is odd or even.

Zero-Phase Filtering for Offline Processing

The constant group delay of a linear phase FIR filter is \(\frac{N-1}{2}\) samples — a non-zero delay. In offline processing where causality is not required, this delay can be eliminated entirely using forward-backward filtering:

  1. Filter the signal in the forward direction to obtain \(y_1[n]\)
  2. Reverse \(y_1[n]\)
  3. Filter again in the forward direction, then reverse again

The net result has zero phase shift and a squared magnitude response:

\[H_{\text{zero-phase}}(e^{j\omega}) = |H(e^{j\omega})|^2 \tag{13}\]

scipy.signal.filtfilt implements this procedure. It doubles the effective filter order and cannot be used for real-time (causal) processing, but it completely eliminates phase distortion.

Nonlinear Phase in IIR Filters

Why IIR Filters Cannot Achieve Linear Phase

An IIR filter with rational transfer function:

\[H(z) = \frac{B(z)}{A(z)} = \frac{\sum_{k=0}^{M} b_k z^{-k}}{1 + \sum_{k=1}^{N} a_k z^{-k}} \tag{14}\]

has poles located inside the unit circle. The phase response of \(H(e^{j\omega})\) can be decomposed as a sum of contributions from each zero and pole:

\[\phi(\omega) = \sum_k \angle(e^{j\omega} - z_k) - \sum_k \angle(e^{j\omega} - p_k) \tag{15}\]

where \(z_k\) and \(p_k\) are the zeros and poles respectively. Since the poles are located asymmetrically relative to the unit circle (they are inside it, with no mirror-image poles outside), the phase response is generally a nonlinear function of \(\omega\) . This is a fundamental property of stable IIR filters — it cannot be circumvented by choice of filter coefficients.

Group Delay Comparison of Classical IIR Filters

The following table summarizes the group delay characteristics of classical IIR filter designs:

FilterTransition bandwidthPassband group delayPrimary use case
ButterworthModerateRelatively smoothGeneral-purpose, moderate phase
Chebyshev Type INarrowMore distortedAmplitude-critical applications
Chebyshev Type IINarrowModerateStopband attenuation priority
Elliptic (Cauer)NarrowestMost distortedMinimum-order design
BesselWideMaximally flatWaveform preservation priority

The Bessel filter is designed with the explicit objective of maximizing the flatness of group delay in the passband (Maximally Flat Group Delay criterion). It achieves the best time-domain waveform fidelity among IIR filters, at the cost of a very gradual transition between passband and stopband.

The Fundamental Trade-off

The relationship between magnitude sharpness and group delay flatness represents a fundamental trade-off in IIR filter design:

  • Sharper amplitude transition → more poles clustered near the unit circle in the transition region → larger group delay variation
  • Flatter group delay → poles more evenly distributed → gentler amplitude roll-off

When both sharp magnitude cutoff and flat group delay are required simultaneously, the only general solution is to use an IIR filter for amplitude shaping followed by an all-pass filter for phase equalization.

All-Pass Filters for Phase Equalization

Definition and Transfer Function

An all-pass filter (APF) has unit magnitude response at all frequencies:

\[|H_{AP}(e^{j\omega})| = 1, \quad \forall \omega \tag{16}\]

while its phase response — and therefore its group delay — can be freely shaped by choice of pole/zero locations. The first-order all-pass filter transfer function is:

\[H_{AP}^{(1)}(z) = \frac{z^{-1} - a}{1 - az^{-1}}, \quad |a| < 1 \tag{17}\]

This filter has a zero at \(z = 1/a\) (outside the unit circle) and a pole at \(z = a\) (inside the unit circle), placed in a mirror-image relationship with respect to the unit circle. This mirror-image placement ensures that the magnitude response is identically 1 for all \(\omega\) .

For a complex pole at \(z = re^{j\theta}\) (\(r < 1\) ), the all-pass filter and its complex conjugate are combined into a second-order real-coefficient all-pass filter:

\[H_{AP}(z) = \frac{z^{-1} - a^*}{1 - az^{-1}} \tag{18}\]

The practical second-order implementation with real coefficients is:

\[H_{AP}^{(2)}(z) = \frac{r^2 - 2r\cos\theta \cdot z^{-1} + z^{-2}}{1 - 2r\cos\theta \cdot z^{-1} + r^2 z^{-2}} \tag{19}\]

Phase and Group Delay of the First-Order All-Pass Filter

For the first-order all-pass filter with real coefficient \(a\) (\(0 < a < 1\) ), the phase response is:

\[\phi_{AP}(\omega) = -\omega - 2\arctan\left(\frac{a\sin\omega}{1 - a\cos\omega}\right) \tag{20}\]

The group delay is:

\[\tau_g^{AP}(\omega) = \frac{1 - a^2}{1 - 2a\cos\omega + a^2} \tag{21}\]

This is always positive (the all-pass filter only adds delay, never advances the signal), and exhibits a peak near \(\omega = 0\) . For the second-order filter, the peak is located near \(\omega = \theta\) (the angle of the pole), and its height is controlled by \(r\) (the pole radius). A value of \(r\) close to 1 produces a sharp, tall peak; a smaller \(r\) gives a broad, low peak.

Group Delay Equalization Strategy

Phase equalization using all-pass filters follows this procedure:

  1. Measure the group delay \(\tau_g^{\text{target}}(\omega)\) of the IIR filter to be equalized using scipy.signal.group_delay.
  2. Identify the frequency regions where \(\tau_g^{\text{target}}(\omega)\) is smallest relative to the maximum — these are the regions where additional delay needs to be added.
  3. Design all-pass sections (choosing \(r\) and \(\theta\) for each second-order section) to add delay in these regions, bringing the total group delay closer to a flat response.
  4. Cascade the all-pass sections after the original IIR filter.

Since \(|H_{AP}| = 1\) , the amplitude response of the original filter is preserved exactly. The equalization is imperfect in general (particularly at the passband edge), but significant improvement is achievable with a few all-pass sections.

Python Implementation

Setup: Library Imports

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, sosfilt, freqz, firwin, group_delay, lfilter, iirfilter
from scipy.signal import cheby1, ellip, bessel

Phase Response: Wrapped vs. Unwrapped

# --- Filter design ---
fs = 1000          # sampling frequency [Hz]
fc = 100           # cutoff frequency [Hz]
N_iir = 8          # IIR filter order
N_fir = 64         # FIR filter tap count

# Butterworth IIR filter (ba form)
b_iir, a_iir = butter(N_iir, fc, fs=fs, output='ba')

# Linear-phase FIR filter (Hamming window)
b_fir = firwin(N_fir, fc, fs=fs, window='hamming')
a_fir = np.array([1.0])

# --- Frequency response ---
worN = 4096
w, H_iir = freqz(b_iir, a_iir, worN=worN, fs=fs)
_, H_fir = freqz(b_fir, a_fir, worN=worN, fs=fs)

# --- Phase computation ---
phase_iir_wrapped = np.angle(H_iir)          # wrapped phase
phase_iir = np.unwrap(np.angle(H_iir))       # unwrapped phase
phase_fir = np.unwrap(np.angle(H_fir))       # FIR (linear phase)

# --- Plot ---
fig, axes = plt.subplots(3, 1, figsize=(10, 10))

# Magnitude
axes[0].plot(w, 20 * np.log10(np.abs(H_iir) + 1e-12),
             label=f'IIR Butterworth N={N_iir}')
axes[0].plot(w, 20 * np.log10(np.abs(H_fir) + 1e-12),
             label=f'FIR Hamming {N_fir} taps', linestyle='--')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('Magnitude Response')
axes[0].set_xlim(0, fs / 2)
axes[0].set_ylim(-80, 5)
axes[0].axvline(fc, color='gray', linestyle=':', alpha=0.5, label='Cutoff')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Wrapped vs unwrapped phase (IIR)
axes[1].plot(w, np.degrees(phase_iir_wrapped),
             label='IIR (wrapped)', color='orange')
axes[1].plot(w, np.degrees(phase_iir),
             label='IIR (unwrapped)', color='blue')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Phase [degrees]')
axes[1].set_title('Phase Response: Wrapped vs Unwrapped (IIR Butterworth)')
axes[1].set_xlim(0, fs / 2)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# IIR vs FIR phase
axes[2].plot(w, np.degrees(phase_iir),
             label=f'IIR Butterworth N={N_iir} (nonlinear)')
axes[2].plot(w, np.degrees(phase_fir),
             label=f'FIR Hamming {N_fir} taps (linear)', linestyle='--')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Phase [degrees]')
axes[2].set_title('Phase Response: IIR vs FIR')
axes[2].set_xlim(0, fs / 2)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Group Delay Calculation and Comparison

# --- Group delay via scipy.signal.group_delay ---
worN = 4096
w_gd, gd_iir = group_delay((b_iir, a_iir), w=worN, fs=fs)
w_gd, gd_fir = group_delay((b_fir, a_fir), w=worN, fs=fs)

# --- Finite difference approximation (for illustration) ---
w_arr, H_iir_arr = freqz(b_iir, a_iir, worN=worN * 2, fs=fs)
omega_arr = w_arr * 2 * np.pi / fs   # Hz → rad/sample
phase_arr = np.unwrap(np.angle(H_iir_arr))

# τ_g(ω) ≈ −Δφ/Δω
gd_diff = -np.diff(phase_arr) / np.diff(omega_arr)
w_diff = (w_arr[:-1] + w_arr[1:]) / 2   # midpoints

# --- Plot ---
fig, axes = plt.subplots(2, 1, figsize=(10, 8))

axes[0].plot(w_gd, gd_iir,
             label=f'IIR Butterworth N={N_iir}', color='blue')
axes[0].plot(w_gd, gd_fir,
             label=f'FIR Hamming {N_fir} taps (constant = {(N_fir-1)/2:.1f} samples)',
             linestyle='--', color='red')
axes[0].axhline((N_fir - 1) / 2, color='red', linestyle=':', alpha=0.5)
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Group Delay [samples]')
axes[0].set_title('Group Delay: IIR vs FIR')
axes[0].set_xlim(0, fs / 2)
axes[0].set_ylim(-5, 50)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(w_diff, gd_diff,
             label='IIR (finite difference)', color='purple', alpha=0.7)
axes[1].plot(w_gd, gd_iir,
             label='IIR (scipy.signal.group_delay)', color='blue', linestyle='--')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Group Delay [samples]')
axes[1].set_title('Group Delay: Finite Difference vs scipy.signal.group_delay')
axes[1].set_xlim(0, fs / 2)
axes[1].set_ylim(-5, 50)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

The FIR filter maintains a constant group delay of exactly \((N-1)/2 = 31.5\) samples across all frequencies, while the Butterworth IIR filter shows a pronounced group delay peak near the cutoff frequency, dropping back to a lower value in the deep stopband.

Waveform Distortion Due to Phase Nonlinearity

# --- Test signal: sum of harmonics ---
fs = 1000
t = np.arange(0, 0.5, 1/fs)

# Multicomponent signal (10, 30, 50, 70, 90 Hz)
freqs_sig = [10, 30, 50, 70, 90]
x = sum(np.sin(2 * np.pi * f * t) / f for f in freqs_sig)

# --- Filtering ---
fc_test = 200   # high cutoff — all components pass through

# IIR Butterworth (nonlinear phase)
b_iir_t, a_iir_t = butter(6, fc_test, fs=fs, output='ba')
y_iir = lfilter(b_iir_t, a_iir_t, x)

# FIR linear-phase filter
b_fir_t = firwin(101, fc_test, fs=fs, window='hamming')
y_fir = lfilter(b_fir_t, [1.0], x)

# Compensate FIR group delay shift
fir_delay = 50   # (101 - 1) / 2 = 50 samples
y_fir_aligned = np.concatenate([y_fir[fir_delay:], np.zeros(fir_delay)])

# --- Plot ---
fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)

axes[0].plot(t, x, 'k', linewidth=1.5)
axes[0].set_title('Original Signal (10 + 30 + 50 + 70 + 90 Hz)')
axes[0].set_ylabel('Amplitude')
axes[0].grid(True, alpha=0.3)

axes[1].plot(t, y_iir, color='blue')
axes[1].plot(t, x, 'k--', alpha=0.3, linewidth=1, label='Original')
axes[1].set_title('After IIR Butterworth (nonlinear phase — waveform distorted)')
axes[1].set_ylabel('Amplitude')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].plot(t, y_fir_aligned, color='red')
axes[2].plot(t, x, 'k--', alpha=0.3, linewidth=1, label='Original')
axes[2].set_title('After FIR linear-phase (delay-corrected — shape preserved)')
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Amplitude')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

All-Pass Filter Design and Phase Equalization

def allpass_filter_2nd(r, theta):
    """
    Compute second-order all-pass filter coefficients.

    Transfer function:
        H_AP(z) = (r^2 - 2r*cos(θ)*z^{-1} + z^{-2})
                / (1 - 2r*cos(θ)*z^{-1} + r^2*z^{-2})

    Parameters
    ----------
    r : float
        Pole radius (0 < r < 1). Larger r → sharper group delay peak.
    theta : float
        Pole angle in radians. Controls peak frequency location.

    Returns
    -------
    b, a : ndarray
        Numerator and denominator coefficients (b = numerator, a = denominator).
    """
    c = 2 * r * np.cos(theta)
    r2 = r ** 2
    b = np.array([r2, -c, 1.0])
    a = np.array([1.0, -c, r2])
    return b, a


# --- All-pass filter characterization ---
worN = 4096
fig, axes = plt.subplots(3, 1, figsize=(10, 10))

r_values = [0.5, 0.7, 0.9]
theta_ap = np.pi / 4   # peak group delay near ω = π/4

for r in r_values:
    b_ap, a_ap = allpass_filter_2nd(r, theta_ap)
    w_ap, H_ap = freqz(b_ap, a_ap, worN=worN, fs=fs)
    w_gd_ap, gd_ap = group_delay((b_ap, a_ap), w=worN, fs=fs)

    axes[0].plot(w_ap, 20 * np.log10(np.abs(H_ap) + 1e-12), label=f'r={r}')
    axes[1].plot(w_ap, np.degrees(np.unwrap(np.angle(H_ap))), label=f'r={r}')
    axes[2].plot(w_gd_ap, gd_ap, label=f'r={r}')

axes[0].set_title('All-Pass Filter: Magnitude Response (must be 0 dB)')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_ylim(-3, 3)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].set_title('All-Pass Filter: Phase Response')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Phase [degrees]')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].set_title('All-Pass Filter: Controllable Group Delay Peak')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Group Delay [samples]')
axes[2].set_ylim(0, 30)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()
# --- Phase equalization: Butterworth + all-pass filter ---
fs = 1000
fc = 100
N_butter = 6

b_butter, a_butter = butter(N_butter, fc, fs=fs, output='ba')
w_hz, gd_butter = group_delay((b_butter, a_butter), w=4096, fs=fs)

# Design all-pass sections to flatten the group delay
r_ap = 0.65
theta_ap1 = 2 * np.pi * (fc * 0.9) / fs   # near cutoff
theta_ap2 = 2 * np.pi * (fc * 0.7) / fs   # slightly below cutoff

b_ap1, a_ap1 = allpass_filter_2nd(r_ap, theta_ap1)
b_ap2, a_ap2 = allpass_filter_2nd(r_ap, theta_ap2)

# Cascade: convolve coefficients
b_total = np.convolve(np.convolve(b_butter, b_ap1), b_ap2)
a_total = np.convolve(np.convolve(a_butter, a_ap1), a_ap2)

w_hz, gd_total = group_delay((b_total, a_total), w=4096, fs=fs)
w_hz, gd_ap1 = group_delay((b_ap1, a_ap1), w=4096, fs=fs)
w_hz, gd_ap2 = group_delay((b_ap2, a_ap2), w=4096, fs=fs)

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

# Magnitude response (all-pass sections should not affect amplitude)
w_hz_mag, H_butter = freqz(b_butter, a_butter, worN=4096, fs=fs)
_, H_total = freqz(b_total, a_total, worN=4096, fs=fs)
axes[0].plot(w_hz_mag,
             20 * np.log10(np.abs(H_butter) + 1e-12),
             label='Butterworth only', color='blue')
axes[0].plot(w_hz_mag,
             20 * np.log10(np.abs(H_total) + 1e-12),
             label='Butterworth + APF (equalized)', color='red', linestyle='--')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('Magnitude Response: All-Pass Equalization Preserves Amplitude')
axes[0].set_xlim(0, fs / 2)
axes[0].set_ylim(-80, 5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Group delay comparison
axes[1].plot(w_hz, gd_butter,
             label=f'Butterworth N={N_butter}', color='blue')
axes[1].plot(w_hz, gd_ap1 + gd_ap2,
             label='APF contribution (ap1 + ap2)', color='green', linestyle='--')
axes[1].plot(w_hz, gd_total,
             label='Butterworth + APF (equalized)', color='red', linewidth=2)
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Group Delay [samples]')
axes[1].set_title('Group Delay Equalization with All-Pass Filter')
axes[1].set_xlim(0, fc * 2)
axes[1].set_ylim(0, 30)
axes[1].axvline(fc, color='gray', linestyle=':', alpha=0.5, label='Cutoff')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

The equalized system (Butterworth + all-pass) shows a flatter group delay profile near the cutoff frequency, while the magnitude response remains unchanged — demonstrating the principle that amplitude and phase shaping can be decoupled using all-pass filters.

IIR Filter Group Delay Comparison

fs = 1000
fc = 100
N = 6

filters = {
    'Butterworth': butter(N, fc, fs=fs, output='ba'),
    'Chebyshev I (1 dB)': cheby1(N, 1, fc, fs=fs, output='ba'),
    'Elliptic (1 dB, 60 dB)': ellip(N, 1, 60, fc, fs=fs, output='ba'),
    'Bessel': bessel(N, fc, norm='mag', fs=fs, output='ba'),
}

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

for name, (b, a) in filters.items():
    w_hz, H = freqz(b, a, worN=4096, fs=fs)
    w_gd, gd = group_delay((b, a), w=4096, fs=fs)
    axes[0].plot(w_hz, 20 * np.log10(np.abs(H) + 1e-12), label=name)
    axes[1].plot(w_gd, gd, label=name)

axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title(f'Magnitude Response: IIR Filter Comparison (N={N}, fc={fc} Hz)')
axes[0].set_xlim(0, 300)
axes[0].set_ylim(-80, 5)
axes[0].axvline(fc, color='gray', linestyle=':', alpha=0.5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Group Delay [samples]')
axes[1].set_title('Group Delay: IIR Filter Comparison')
axes[1].set_xlim(0, 300)
axes[1].set_ylim(0, 50)
axes[1].axvline(fc, color='gray', linestyle=':', alpha=0.5)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

The Bessel filter exhibits the flattest group delay across the passband, confirming its Maximally Flat Group Delay property. The elliptic filter achieves the steepest magnitude roll-off but at the cost of the largest group delay distortion — illustrating the fundamental sharpness-versus-phase-linearity trade-off.

Summary

ConceptDefinition / PropertyPractical significance
Phase spectrum\(\phi(\omega) = \angle H(e^{j\omega})\) ; requires unwrapping for continuityDetermines how each frequency component is time-shifted
Group delay\(\tau_g(\omega) = -d\phi/d\omega\) ; units: samplesQuantifies frequency-dependent signal packet delay
Linear phase\(\phi(\omega) = -\omega\tau_0\) ; constant group delay \(\tau_0\)Output is a pure time-shifted copy of the input
FIR linear phaseSymmetric coefficients \(h[n]=h[N-1-n]\) ; \(\tau_g = (N-1)/2\) samplesExact linear phase guaranteed; magnitude and phase independently designable
IIR nonlinear phaseFundamental property of stable IIR filters with poles inside the unit circleTrade-off against amplitude sharpness is unavoidable
All-pass filter\(\|H_{AP}\| = 1\) for all \(\omega\) ; phase freely controllablePhase equalization without changing amplitude
Phase equalizationCascade IIR + all-pass; optimize for flat total group delayReduces waveform distortion in real-time IIR systems

Practical Filter Selection Guide

RequirementRecommended approach
Waveform shape preservation is critical (ECG, impulse response)Linear-phase FIR filter (firwin)
Sharp magnitude cutoff, phase distortion acceptableElliptic IIR filter (ellip)
Flat group delay in passband with IIRBessel filter (bessel)
IIR filter with phase distortion reductionIIR + all-pass equalization
Offline processing, zero phase shift requiredfiltfilt / sosfiltfilt (forward-backward filtering)
Real-time processing with minimum total delayLow-order IIR (Butterworth) + group delay monitoring

References