Exploring Bandpass Transformation

Table Of Contents

In this article, we present a method for applying the bandpass transform to any second-order filter decomposition. This approach results in a new second-order filter decomposition with twice the initial order. It can effectively convert elliptic and Chebyshev filters into their bandpass versions or transform shelf filters into peaks.

Understanding the Bandpass Transform

The bandpass transform is defined as:

$$s \leftarrow Q\left(s+\frac{1}{s}\right)$$

By applying this transform to a biquad transfer function:

$$H(s) = \frac{b_0 s^2 + b_1 s + b_2}{a_0 s^2 + a_1 s + a_2}$$

we obtain the transformed transfert function:

$$H_{BP}(s) = \frac{b_{0}s^{4} + \frac{b_{1}}{Q}s^{3} + \left(2b_{0} + \frac{b_{2}}{Q^{2}}\right)s^{2} + \frac{b_{1}}{Q}s + b_{0}}{a_{0}s^{4} + \frac{a_{1}}{Q}s^{3} + \left(2a_{0} + \frac{a_{2}}{Q^{2}}\right)s^{2} + \frac{a_{1}}{Q}s + a_{0}}$$

To simplify this expression, we can rewrite it in the canonical form as follows:

$$H_{BP}(s) = \frac{b_{0}}{a_{0}} \left( \frac{s^{4} + \frac{b_{1}}{b_{0}Q}s^{3} + \left(2 + \frac{b_{2}}{b_{0}Q^{2}}\right)s^{2} + \frac{b_{1}}{b_{0}Q}s + 1}{s^{4} + \frac{a_{1}}{a_{0}Q}s^{3} + \left(2 + \frac{a_{2}}{a_{0}Q^{2}}\right)s^{2} + \frac{a_{1}}{a_{0}Q}s + 1} \right), \quad b_0 \neq 0, \ a_0 \neq 0$$

Our goal is to factor \(H_{BP}(s)\) into two second-order sections, expressed as:

$$H_{BP}(s) = H_1(s) H_2(s)$$

where both \(H_1(s)\) and \(H_2(s)\) have real coefficients. This requires solving the roots of the numerator and denominator polynomials of \(H_{BP}(s)\) and pairing any conjugate roots.

Solving the Fourth-Order Polynomials

The numerator and denominator polynomials of \(H_{BP}(s)\) can be represented in the form:

$$P(x) = x^4 + cx^3 + dx^2 + cx + 1$$

where \(c, d \in \mathbb{R}\).

The four roots \(r_k\) of \(P(x)\) are calculated as follows¹:

$$r_1 = \frac{u + v_1 - c}{4} \quad r_2 = \frac{u - v_1 - c}{4}$$ $$r_3 = \frac{-u + v_2 - c}{4} \quad r_4 = \frac{-u - v_2 - c}{4}$$

with the variables defined as:

$$v_1 = \sqrt{2c^2 - 4d - 8 - cu}$$ $$v_2 = \sqrt{2c^2 - 4d - 8 + cu}$$ $$u = \sqrt{c^2 - 4d + 8}$$

Root Pairing

To factor the fourth-order polynomial into two second-order real polynomials, we need to pair conjugate roots when applicable. Pairing can be complex, as some pairs might swap during modulation, potentially causing unsolvable audio artifacts. Nevertheless, this issue is relatively minor compared to managing the numerical stability of a single fourth-order filter, which is more likely to propagate errors due to its long recursive structure.

The appropriate pairing of the roots \(r_1, r_2, r_3, r_4\) depends on the coefficients \(c\) and \(d\) of polynomial \(P(x)\):

• If \(c^{2} - 4d + 8 \geq 0\):

  • If \(2c^{2} - 4d - 8 - cu \geq 0\) and \(2c^{2} - 4d - 8 + cu \geq 0\), then both \(v_1\) and \(v_2\) are real, allowing flexible pairing.
  • Otherwise (either \(v_1\) or \(v_2\) is purely imaginary), pair roots as \((r_1, r_2)\) and \((r_3, r_4)\).

• If \(c^{2} - 4d + 8 < 0\):

  • If \(c \neq 0\), then \(v_1 = \overline{v_2}\), leading to conjugate pairs \((r_1, r_3)\) and \((r_2, r_4)\).
  • If \(c = 0\):
    • If \(2c^{2} - 4d - 8 \geq 0\) (\(v_1\) and \(v_2\) are real), then pair conjugates \((r_1, r_3)\) and \((r_2, r_4)\).
    • Otherwise (\(v_1\) or \(v_2\) is purely imaginary), pair conjugates \((r_1, r_4)\) and \((r_2, r_3)\).

Let \(r_A, r_B, r_C, r_D\) represent the reordered roots, grouped into pairs \((r_A, r_B)\) and \((r_C, r_D)\). With this arrangement, the polynomial \( P(x) \) can now be factored as:

$$P(x) = \left(x^{2} - (r_A + r_B)x + r_Ar_B\right)\left(x^{2} - (r_C + r_D)x + r_Cr_D\right)$$

Applying Factorization to \(H_{BP}\)

Now, we can apply the factorization obtained from \(P(x)\) to both the numerator and denominator of \(H_{BP}(s)\):

$$ \begin{aligned} H_{BP}(s)&=\frac{b_{0}}{a_{0}}\left(\frac{s^{2}-\left(z_{1}+z_{2}\right)s+z_{1}z_{2}}{s^{2}-\left(p_{1}+p_{2}\right)s+p_{1}p_{2}}\right)\left(\frac{s^{2}-\left(z_{3}+z_{4}\right)s+z_{3}z_{4}}{s^{2}-\left(p_{3}+p_{4}\right)s+p_{3}p_{4}}\right) \newline &=\left(\frac{b_{0}s^{2}-b_{0}\left(z_{1}+z_{2}\right)s+b_{0}z_{1}z_{2}}{a_{0}s^{2}-a_{0}\left(p_{1}+p_{2}\right)s+a_{0}p_{1}p_{2}}\right)\left(\frac{s^{2}-\left(z_{3}+z_{4}\right)s+z_{3}z_{4}}{s^{2}-\left(p_{3}+p_{4}\right)s+p_{3}p_{4}}\right) \end{aligned} $$

Special Cases

Case 1: \(b_0 = 0\) or \(a_0 = 0\)

When \(b_0 = 0\), the function simplifies to:

$$H_{BP}(s)=\frac{\left(b_{1}s^{2}+\frac{b_{2}}{Q}s+b_{1}\right)\frac{s}{Q}}{a_{0}s^{4}+\frac{a_{1}}{Q}s^{3}+\left(2a_{0}+\frac{a_{2}}{Q^{2}}\right)s^{2}+\frac{a_{1}}{Q}s+a_{0}}$$

In this case, the factorization of the numerator is already explicit. The same applies for the case \(a_0 = 0\).

Case 2: \(b_1 = b_2 = 0\) or \(a_1 = a_2 = 0\)

For instance, applying the bandpass transform to an analog high-pass filter where \(b_1 = b_2 = 0\) yields:

$$H_{BP}(s)=\frac{b_{0}\left(s^{2}+1\right)^{2}}{a_{0}s^{4}+\frac{a_{1}}{Q}s^{3}+\left(2a_{0}+\frac{a_{2}}{Q^{2}}\right)s^{2}+\frac{a_{1}}{Q}s+a_{0}}$$

Note that this case is already addressed by the root-solving method, however, using this formula provides a faster computation. The same applies for the case \(a_1 = a_2 = 0\).

Implementation

The code below demonstrates the application of the bandpass transform to several common filters.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# Solve x^4 + cx^3 + dx^2 + xc + 1 = (x^2 + a1*x + a2)(x^2 + a3*x + a4) for a1, a2, a3, a4 reals
# Input : c, d
# Output : a1, a2, a3, a4
def factorize_4order_palindromic_polynomial(c, d):
    # Find roots
    uu = c*c - 4*d + 8
    u = np.sqrt(uu + 0j)
    a = 2*c*c - 4*d - 8
    v1 = np.sqrt(a - 2*c*u)
    v2 = np.sqrt(a + 2*c*u)
    r1 = (u + v1 - c) / 4
    r2 = (u - v1 - c) / 4
    r3 = (- u + v2 - c) / 4
    r4 = (- u - v2 - c) / 4
    # Pair roots
    if (uu < 0):
        if c != 0 or a >= 0:
            r1, r2, r3, r4 = r1, r3, r2, r4
        else:
            r1, r2, r3, r4 = r1, r4, r2, r3
    assert(np.isclose(np.imag(r1 + r2), 0))
    assert(np.isclose(np.imag(r3 + r4), 0))
    # Compute coefficients
    a1, a2 = np.real(-r1 - r2), np.real(r1 * r2)
    a3, a4 = np.real(-r3 - r4), np.real(r3 * r4)
    return a1, a2, a3, a4

def bandpass_transform_sos(sos, Q):
    sos_bp = []
    for section in sos:
        b0, b1, b2, a0, a1, a2 = section
        if a0 != 0:
            a = factorize_4order_palindromic_polynomial(a1 / (a0 * Q), 2 + a2 / (a0 * Q * Q))
            a = a0, a0*a[0], a0*a[1], 1, a[2], a[3]
        else:
            a = a1, a2/Q, a1, 0, 1/Q, 0
        if b0 != 0:
            b = factorize_4order_palindromic_polynomial(b1 / (b0 * Q), 2 + b2 / (b0 * Q * Q))
            b = b0, b0*b[0], b0*b[1], 1, b[2], b[3]
        else:
            b = b1, b2/Q, b1, 0, 1/Q, 0
        sos_bp.append([b[0], b[1], b[2], a[0], a[1], a[2]])
        sos_bp.append([b[3], b[4], b[5], a[3], a[4], a[5]])
    return sos_bp

def gain_sos(sos, w):
    g = np.ones(len(w))
    for section in sos:
        g *= abs(signal.freqs(section[:3], section[3:], w)[1])
    return g

plt.figure(figsize=(11, 18))
i = 1
def plot(sos, title):
    global i
    w = np.logspace(-5, 5, 5000, base = 2)

    # Prototype filter
    ax = plt.subplot(12, 3, i)
    plt.plot(np.log2(w), np.log2(gain_sos(sos, w)), color = "white", alpha = 1)
    plt.xlim(-5,5)
    plt.ylim(-5,2)
    if (i <= 3):
        plt.title("Prototype filter")
    plt.setp(ax.get_xticklabels(), visible=False)
    plt.setp(ax.get_yticklabels(), visible=False)
    ax.tick_params(axis='both', which='both', length=0)
    plt.ylabel(title)
    plt.axis('off')
    plt.tight_layout()
    i += 1

    # Band pass filter
    ax = plt.subplot(12, 3, i)
    for j, Q in enumerate(np.logspace(-0.5, 0.5, 21)):
        plt.plot(np.log2(w), np.log2(gain_sos(bandpass_transform_sos(sos, Q), w)), color = "grey", alpha = 0.5, linewidth=0.8)
    plt.plot(np.log2(w), np.log2(gain_sos(bandpass_transform_sos(sos, 1), w)), color = "white", alpha = 1)
    plt.xlim(-5,5)
    plt.ylim(-5,2)
    if (i <= 3):
        plt.title("Band pass filter")
    plt.axis('off')
    plt.tight_layout()
    i += 1

    # Notch filter
    sos[:,:3], sos[:,3:] = sos[:,2::-1], sos[:,:2:-1] # highpass_transform_sos: Invert a/b coeffs
    ax = plt.subplot(12, 3, i)
    for j, Q in enumerate(np.logspace(-0.5, 0.5, 21)):
        plt.plot(np.log2(w), np.log2(gain_sos(bandpass_transform_sos(sos, Q), w)), color = "grey", alpha = 0.5, linewidth=0.8)
    plt.plot(np.log2(w), np.log2(gain_sos(bandpass_transform_sos(sos, 1), w)), color = "white", alpha = 1)
    plt.xlim(-5,5)
    plt.ylim(-5,2)
    if (i <= 3):
        plt.title("Notch filter")
    plt.axis('off')
    plt.tight_layout()
    i += 1
    
plot(np.array([[0, 0, 1, 0, 1, 1]]), "First Order")
plot(np.array([[0, 0, 1, 1, 0.5, 1]]), "Second Order")
sos = signal.butter(16, 1, 'low', output = "sos", analog=True)
sos[:,:3] = 0, 0, 1 # Quick fix of a scipy issue
plot(sos, "Butterworth")
plot(signal.ellip(16, 0.02, 90, 1, 'low', output = "sos", analog=True), "Brickwall")
sos = signal.cheby1(8, 15, 1, 'low', output = "sos", analog=True)
sos[:,0], sos[:,2] = 0, sos[:,0] + sos[:,2] # Quick fix of a scipy issue
plot(sos, "Chebyshev 1")
plot(signal.cheby2(8, 15, 1, 'low', output = "sos", analog=True), "Chebyshev 2")
plot(signal.ellip(4, 6, 22, 1, 'low', output = "sos", analog=True), "Elliptic")
sos = [[0, 0, 1, 1, 1/2, 1]]
sos = bandpass_transform_sos(sos, 1)
plot(np.array(sos), "Resonant Band Pass")
sos = bandpass_transform_sos(sos, 1)
plot(np.array(sos), "Resonant Band Pass X2")
sos = [[1, 0, 0, 1, 1/2, 1]]
sos = bandpass_transform_sos(sos, 1)
plot(np.array(sos), "Resonant Notch")
sos = bandpass_transform_sos(sos, 1)
plot(np.array(sos), "Resonant Notch X2")
plt.show()