Louis Couka

Chebyshev Filters: From Theory to Musical Applications

Published: 24/03/2022 | Author: Louis Couka

Table Of Contents

Introduction

This article introduces an adaptation of the Chebyshev filter for musical applications. The approach focuses on defining a user-friendly resonance parameter, normalizing both the cutoff frequency and output gain, and applying filter transformations to conform to standard filter types: low-pass, high-pass, band-pass, and notch.

Chebyshev Polynomials

To design the Chebyshev filter, we use the Chebyshev polynomial \(T_n\) of order \(n\). Various definitions for \(T_n\) can be found here. For efficient implementations, one might opt for a recursive approach or precompute the polynomials.

For now, we concentrate on the essential requirements, using the following definition for \(T_n\) and its inverse \(T_n^{-1}\):

$$ T_n(x) = \begin{cases} \cos(n \arccos x) & \text{if } |x| \le 1, \newline \cosh(n \operatorname{arcosh} x) & \text{if } x \ge 1 \end{cases} $$

$$ T_n^{-1}(x) = \begin{cases} \cos\left(\frac{\arccos(x)}{n}\right) & \text{if } |x| \le 1, \newline \cosh\left(\frac{\operatorname{arcosh}(x)}{n}\right) & \text{if } x \ge 1 \end{cases} $$

A plot below illustrates the behavior of these functions.

Python
 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

import numpy as np
import matplotlib.pyplot as plt

# Chebyshev polynomial function T_n(x) from https://en.wikipedia.org/wiki/Chebyshev_polynomials
# This implementation is valid for input values in the range [-1, oo]
def T(n, w):  # Equivalent to ((w - sqrt(w^2 - 1))^n + (w + sqrt(w^2 - 1))^n) / 2
    return np.piecewise(
        w,
        [(-1 <= w) & (w <= 1), w > 1],
        [lambda w: np.cos(n * np.arccos(w)), lambda w: np.cosh(n * np.arccosh(w))]
    )

# Main branch of the inverse Chebyshev function T_n^{-1}(y), valid for [-1, oo]
def T_inv(n, y):
    return np.piecewise(
        y,
        [(-1 <= y) & (y <= 1), y > 1],
        [lambda y: np.cos(np.arccos(y) / n), lambda y: np.cosh(np.arccosh(y) / n)]
    )

n = 8
x = np.linspace(-1, 4, 1000)
plt.plot(x, T(n, x), label="\\(T_n(x)\\)")
plt.plot(x, T_inv(n, x), label="\\(T_n^{-1}(x)\\)")
plt.title(f"Chebyshev Polynomial and Its Inverse (n={n})")
plt.xlabel("x")
plt.ylabel("y")
plt.ylim(-2, 2)
plt.xlim(-4, 4)
plt.legend()
plt.show()

Python - Console Output

Magnitude and Gain Normalization

The gain of the Chebyshev filter is traditionally expressed as:

$$G_{\text{down}}(x) = \sqrt{\frac{1}{1 + \epsilon^2 T(x)^2}}$$

However, to ensure that the filter’s oscillations remain positive and the gain at DC (zero frequency) equals 1, we adjust the formula to:

$$\boxed{G_{\text{up}}(x) = \sqrt{\frac{1 + \epsilon^2}{1 + \epsilon^2 T(x)^2}}}$$

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

w = np.logspace(-5, 5, 50000, base=2) # Generate 10 octaves (~20Hz to 20kHz)
eps = 4  # Example value for epsilon
g_up = np.sqrt((1 + eps**2) / (1 + (eps * np.vectorize(T)(n, w))**2))  # Adjusted gain (G_up)
g_down = np.sqrt(1 / (1 + (eps * np.vectorize(T)(n, w))**2))  # Traditional gain (G_down)

plt.plot(np.log2(w), 20 * np.log10(g_up), label="\\(G_{up}(\omega)\\)")
plt.plot(np.log2(w), 20 * np.log10(g_down), label="\\(G_{down}(\omega)\\)")
plt.title("Gain Normalization")
plt.xlabel('Frequency (octaves)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-30, 30)
plt.legend()
plt.show()

Python - Console Output

Here’s an improved version of the text:

Definition of a Quality Factor (Q)

To control the filter’s resonance using a traditional Q factor, we use:

$$\boxed{\varepsilon = \frac{2Q^{2} - 1}{\sqrt{4Q^{2} - 1}}}$$

This equation ensures that the peak magnitude remains consistent regardless of the filter order and matches that of a standard second-order filter. Therefore, for a second-order Chebyshev filter, the Q control will behave the same as in a typical second-order filter.

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

Q = 6  # Quality factor

# Calculate epsilon based on the Q factor
eps = (2 * Q**2 - 1) / (4 * Q**2 - 1)**0.5

# Loop through specified filter orders
for n in [8, 2]:
    # Calculate and plot the gain g
    g = np.sqrt((1 + eps**2) / (1 + (eps * np.vectorize(T)(n, w))**2))
    plt.plot(np.log2(w), 20 * np.log10(g), label="\\(n=" + str(n) + "\\)")

plt.title("Q-Factor Definition (Q=" + str(Q) + ")")
plt.xlabel('Frequency (octaves)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-30, 30)

# Add the Q factor in dB to the y-axis ticks
Q_dB = 20 * np.log10(Q) # Calculate the Q factor in dB
current_ticks = plt.yticks()[0]
new_ticks = np.append(current_ticks, Q_dB)
plt.yticks(new_ticks, labels=[f"{int(tick)}" if tick != Q_dB else r"\\(20\log_{10}(Q)\\)" for tick in new_ticks])

plt.legend()
plt.show()

Python - Console Output

Normalization of the Cutoff Frequency

Several strategies can be employed to normalize the cutoff frequency of the filter.

The first and most intuitive approach is to normalize the frequency so that the asymptotic gain as \(\omega \to +\infty\) matches that of a Butterworth filter. This ensures that the modulation of the Q-factor only alters the shape of the magnitude around \(\omega_c.\) For our filter, this normalization involve taking as its cut frequency:

$$\omega_{\text{norm1}} = \left(\frac{\varepsilon^{2}4^{n-1}}{1 + \varepsilon^{2}}\right)^{-\frac{1}{2n}}$$

However, in practice, this normalization can result in significant jumps in the frequency response when modulating the Q factor.

A second approach normalizes the frequency so that the filter’s gain at (\omega = 1) matches that of a second-order filter. Since the Q factor corresponds to the linear gain at (\omega = 1), it is sensible to maintain this property as a constraint. For this reason, and because this method does not produce the significant jumps seen with the first method, we choose it for normalization. The cutoff frequency is defined as:

$$\boxed{\omega_{\text{norm2}} = T^{-1}\left(\frac{1}{2Q^{2} - 1}\right)}$$

Finally, another closely related approach normalizes the frequency so that the most resonant section of the filter has a cutoff frequency of 1. This method simplifies the calculations (hint: (\sinh/\cosh = \tanh)) and is given by the following relation, where (\tilde{p}_k) denotes the poles of the Chebyshev filter, which we will define later:

$$ \omega_{\text{norm3}} = |\tilde{p}_{0}| $$

Calculation of Poles and Zeros and Normalization of the Cutoff Frequency

The poles and zeros of a Chebyshev filter are given by the following equations^1:

$$ \tilde{p}_{k} = -\sinh(u)\sin(\theta _{k}) + i\cosh(u)\cos(\theta _{k}) $$

where

$$ u = \frac{1}{n}\operatorname{arcsinh}\left(\frac{1}{\varepsilon}\right), \quad \theta_{k} = \frac{\pi(2k + 1)}{2n}, \quad k = [0, \ldots, n - 1] $$

After applying frequency normalization, we can express our poles as:

$$ \boxed{ p_k = \frac{\tilde{p}_{k}}{\omega _{\text{norm2}}} } $$

Python
 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

import numpy as np
import matplotlib.pyplot as plt

w = np.logspace(-5, 5, 50000, base=2) # Generate 10 octaves (~20Hz to 20kHz)
n = 8  # Filter order

# Set up the figure for poles
plt.figure(figsize=(8, 8))
plt.subplot(111, projection='polar')

# Loop over various Q factors
for Q in [0.5**0.5, 0.71, 0.8, 2, 4]:
    # Calculate the poles for a standard Chebyshev filter
    eps = (2 * Q**2 - 1) / (4 * Q**2 - 1)**0.5 # Compute epsilon to keep the right dB peak
    m = np.arange(0, n // 2)  # Half the number of poles
    theta = np.pi / 2 * (2 * m + 1) / n  # Calculate angles for poles
    ar = np.arcsinh(1 / eps) / n  # Compute inverse hyperbolic sine
    p = -np.sinh(ar) * np.sin(theta) + 1j * np.cosh(ar) * np.cos(theta)  # Calculate pole locations

    # Normalize the poles based on the chosen method
    w_norm = T_inv(n, 1 / (2 * Q**2 - 1))  # Normalize to match Butterworth filter gain at w = 1
    p /= w_norm  # Apply normalization
    p = np.append(p, np.conj(p))  # Include complex conjugate poles

    # Plot pole locations in the s-plane
    plt.figure(1)
    title = 'Q = ' + str(Q if Q != 0.5**0.5 else "\\(0.5^{0.5}\\)")
    plt.plot(np.angle(p), np.abs(p), 'x', markersize=10, label=title)

    # Plot corresponding magnitude response
    s = 1j * w
    H = np.ones_like(s)
    for p_ in p:
        H *= 1 / (1 - s / p_)
    plt.figure(2)
    plt.plot(np.log2(w), 20 * np.log10(np.abs(H)), label=title)

# Finalize the pole location plot
plt.figure(1)
plt.title("Poles Location on the s-plane")
plt.legend()

# Finalize the magnitude response plot
plt.figure(2)
plt.title("Corresponding Magnitude Response")
plt.xlabel('Frequency (octaves)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-30, 30)
plt.legend()

plt.show()

Python - Console Output

Python - Console Output

Calculation of Extrema

The frequencies of the local extrema are given by the expression:

$$w_{m} = \frac{1}{w_{2}} \cos\left(\frac{m\pi}{2n}\right), \quad m = \left[1, \ldots, n\right]$$

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

# Plot the overall magnitude response with a dashed line
plt.plot(np.log2(w), 20 * np.log10(abs(H)), "--", alpha=0.2, label="Magnitude")

# Calculate the frequencies of local extrema
w_extremums = np.cos(np.arange(1, n) * np.pi / (2 * n)) / w_norm

# Plot the local extrema on the same graph
s_extremums = 1j * w_extremums
H_extremums = np.ones_like(s_extremums)
for p_ in p:
    H_extremums *= 1 / (1 - s_extremums / p_)
plt.plot(np.log2(w_extremums), 20 * np.log10(abs(H_extremums)), "+", markersize=10, c="C0", label="Local Extremums")
plt.title("Local Extrema")
plt.xlabel('Frequency (octaves)')
plt.ylabel('Magnitude (dB)')
plt.legend()
plt.ylim(-30, 30)
plt.show()

Python - Console Output

Factorization into Second-Order Sections

In digital implementations, complex poles require the use of second-order sections (SOS) to maintain real coefficients and enhance numerical stability. Each section corresponds to a pair of conjugate poles. The overall transfer function can be expressed as the product of \(n/2\) second-order sections:

$$ \boxed{H_{LP}(s) = \prod_{k=0}^{n/2-1} \frac{1}{a_{k,0}s^{2} + a_{k,1}s + a_{k,2}}} $$

where the coefficients \(a_{k,0}\), \(a_{k,1}\), and \(a_{k,2}\) are defined as:

$$ \begin{aligned} a_{k,0} &= 1 \newline a_{k,1} &= -(p_k + \overline{p_k}) &= -2 , \operatorname{Re}(p_k) \newline a_{k,2} &= p_k \cdot \overline{p_k} &= \operatorname{Re}(p_k)^2 + \operatorname{Im}(p_k)^2 \end{aligned} $$

Highpass, Bandpass, and Notch Versions

Highpass

The highpass version is straightforward: we simply replace the numerator with \(s^2\) for each subsection:

$$ \boxed{H_{HP}(s) = \prod_{k=0}^{n/2-1} \frac{s^2}{a_{k,0}s^{2} + a_{k,1}s + a_{k,2}}} $$

Bandpass

The bandpass transformation is defined as:

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

Applying this to \(H_{LP}\) yields:

$$ \boxed{H_{BP}(s) = \prod_{k=0}^{n/2-1} \frac{\frac{a_{k,2}}{Q^{2}}s^{2}}{s^{4} + \frac{a_{k,1}}{Q}s^{3} + \left(\frac{a_{k,2}}{Q^{2}} + 2\right)s^{2} + \frac{a_{k,1}}{Q}s + 1}}. $$

The resulting fourth-order denominator can be factored into two second-order sections using a root solver described here.

Notch

The notch transformation is given by:

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

Applying this to \(H_{LP}\) results in:

$$ \boxed{H_{NT}(s) = \prod_{k=0}^{n/2-1} \frac{\left(s^{2} + 1\right)^{2}}{s^{4} + \frac{a_{k,1}}{Q}s^{3} + \left(\frac{a_{k,2}}{Q^{2}} + 2\right)s^{2} + \frac{a_{k,1}}{Q}s + 1}}. $$

Similarly to the bandpass version, the fourth-order denominator can also be factored into two second-order sections using a root solver described here.

Factorization to Second-Order Sections

In the digital implementation, complex poles necessitate the use of second-order sections (SOS) to handle real coefficients and ensure improved numerical stability. Each section corresponds to a pair of conjugate poles. The complete transfer function can then be expressed as the product of \(n/2\) second-order sections as follow:

$$ \boxed{H_{LP}(s) = \prod_{k=0}^{n/2-1} \frac{1}{a_{k,0}s^{2} + a_{k,1}s + a_{k,2}}} $$

where the coefficients \(a_{k,0}\), \(a_{k,1}\), and \(a_{k,2}\) are given by:

$$ \begin{aligned} a_{k,0} &= 1 \newline a_{k,1} &= -(p_k + \overline{p_k}) &= -2 , \operatorname{Re}(p_k) \newline a_{k,2} &= p_k \cdot \overline{p_k} &= \operatorname{Re}(p_k)^2 + \operatorname{Im}(p_k)^2 \end{aligned} $$

Highpass, Bandpass and Notch version

Highpass

The Highpass version is trivial as we just need to replace the numerator by \(s^2\) for each sub section like

$$ \boxed{H_{HP}(s) = \prod_{k=0}^{n/2-1} \frac{s^2}{a_{k,0}s^{2} + a_{k,1}s + a_{k,2}}} $$

Bandpass

The Band Pass transform is given by:

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

Applying it on \(H_{LP}\) we get:

$$ \boxed{H_{BP}(s) = \prod_{k=0}^{n/2-1} \frac{\frac{a_{k,2}}{Q^{2}}s^{2}}{s^{4}+\frac{a_{k,1}}{Q}s^{3}+\left(\frac{a_{k,2}}{Q^{2}}+2\right)s^{2}+\frac{a_{k,1}}{Q}s+1}} $$

the 4-order denominator may be factorized in two 2-order sections using the root solver from here.

Notch

The Notch transform is given by:

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

Applying it on \(H_{LP}\) we get:

$$ \boxed{H_{NT}(s)=\prod_{k=0}^{n/2-1}\frac{\left(s^{2}+1\right)^{2}}{s^{4}+\frac{a_{k,1}}{Q}s^{3}+\left(\frac{a_{k,2}}{Q^{2}}+2\right)s^{2}+\frac{a_{k,1}}{Q}s+1}} $$

As well as for the bandpass version, the 4-order denominator may be factorized in two 2-order sections using the root solver from here.

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# Bandpass and Notch versions
Q_bw = 0.3
dw = 1/Q_bw
H_bp, H_nt = np.ones_like(s), np.ones_like(s)
for i in range(n//2):
    a0, a1, a2 = 1, -p[i]-np.conj(p[i]), p[i]*np.conj(p[i])
    H_bp *= (a2*dw**2*s**2) / (s**4 + a1*dw*s**3 + (a2*dw**2+2)*s**2 + a1*dw*s + 1) # ie : H_bp *= 1/(a0*s_bp**2/a2 + a1*s_bp/a2 + 1) with s_bp = Q_bw*(s+1/s)
    H_nt *= (s**4 + 2*s**2 + 1) / (s**4 + a1/a2*dw*s**3 + (a0/a2*dw**2+2)*s**2 + a1/a2*dw*s + 1) # ie : H_nt *= 1/(a0*s_nt**2/a2 + a1*s_nt/a2 + 1) with s_nt = 1/s_bp

plt.plot(np.log2(w), 20*np.log10(abs(H_bp)), label = "Chebyshev BP")
plt.plot(np.log2(w), 20*np.log10(abs(H_nt)), label = "Chebyshev NT")
plt.title("Chebyshev Bandpass and Notch Versions")
plt.xlabel('Frequency (octaves)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-30, 30)
plt.legend()
plt.show()

Python - Console Output

Adding Underdamping Capability

Currently, the Chebyshev filter does not support a \(Q\) value less than or equal to \(\frac{\sqrt{2}}{2}\). To enhance our musical Chebyshev filter’s flexibility, we aim to enable underdamping, aligning its behavior with our other filters.

Strategy

Our proposed strategy is to utilize the Butterworth resonant filter formula for \(Q \leq \frac{\sqrt{2}}{2}\), and the existing Chebyshev filter formula for \(Q > \frac{\sqrt{2}}{2}\). As \(Q\) approaches \(\frac{\sqrt{2}}{2}\), the Chebyshev filter behaves increasingly like a Butterworth filter. This approach ensures a smooth transition between the two filter types. Notably, we retain the property that for second-order filters, our musical Chebyshev acts as a simple second-order filter with an appropriate definition of \(Q\).

Smoothing the Transition

Although the transition between the Chebyshev and Butterworth resonant filters is continuous, a noticeable jump can occur in practice at higher orders. Thus, implementing a smooth morphing between the two filter types is advantageous.

To achieve this, we decompose the \(Q\) factor into two parameters: \(Q_{\text{res}}\) and \(Q_{\text{damp}}\), ensuring their application compensates for each other by maintaining a consistent gain at the cutoff frequency. This constraint leads to the relationship:

$$ Q = Q_{\text{res}} \cdot Q_{\text{damp}} $$

Here, \(Q_{\text{res}}\) is the \(Q\) factor used for evaluating the Chebyshev filter, and we want \(Q_{\text{res}} > \frac{\sqrt{2}}{2}\). Conversely, \(Q_{\text{damp}}\) represents the additional resonance applied to the most resonant section, similar to the approach used for Butterworth resonant filters, described here. We desire \(Q_{\text{damp}} < \frac{\sqrt{2}}{2}\).

The following functions \(f_1\) and \(f_2\) are suitable candidates for this transition, selected for their mathematical simplicity and symmetric behavior in log-log scaling, as the end user will control the \(Q\) factor through a logarithmic mapping. The empirical value of \(k\) defines the range of \(Q\) values over which the transition occurs; we set it to 8 for efficient computational performance:

$$ f_{1}(x) = \left( x^{\alpha} + 1 \right)^{\frac{1}{\alpha}} $$

$$ f_{2}(x) = \left( x^{-\alpha} + 1 \right)^{-\frac{1}{\alpha}} = \frac{x}{f_{1}(x)} $$

With the transition pivot centered on \(Q_0 = \frac{\sqrt{2}}{2}\), we obtain:

$$ \boxed{ \begin{aligned} Q_{\text{res}} &= \left( Q^{\alpha} + Q_{0}^{\alpha} \right)^{\frac{1}{\alpha}} \ Q_{\text{damp}} &= \frac{Q}{Q_{\text{res}}} \end{aligned} } $$

Python
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

def f1(x, alpha):
    return (x**alpha + 1)**(1/alpha)

def f2(x, alpha):
    return (x**(-alpha) + 1)**(-1/alpha)

alpha = 8
x = np.logspace(-0.5, 0.5, 400)
plt.plot(x, f1(x, alpha), label=r"\\(f_1(x)=(x^\alpha + 1)^{{1/\alpha}}\\)")
plt.plot(x, f2(x, alpha), label=r"\\(f_2(x)=(x^{-\alpha} + 1)^{{-1/\alpha}}\\)")
plt.plot(x, x, label=r"\\(f_1(x)f_2(x)=x\\)", linestyle="--", c="white")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("x (log scale)")
plt.ylabel("y (log scale)")
plt.title(r"Log-Log Plots of \\(\left(x^{\alpha}+1\right)^{\frac{1}{\alpha}}\\) and \\(\left(x^{-\alpha}+1\right)^{-\frac{1}{\alpha}}\\) for \\(\alpha = " + str(alpha) + "\\)")
plt.legend()
plt.grid(True, which="both", linestyle='--', linewidth=0.5)
plt.show()

Python - Console Output

Implementation

The code below summarizes the calculation of the second-order sections (SOS) for the musical Chebyshev prototype, as well as the computation of its magnitude without directly using its transfer function. It includes the underdamped capability we just explained, with the smooth transition to the butterworth resonant filter that we studied [here](CITER LE BUTTERWORTH RESONANT).

Python
 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

n = 8
w = np.logspace(-5, 5, 50000, base = 2) # 10 octaves ~= 20 to 20k
s = 1j * w
k = 8
for j, Q in enumerate(np.exp(np.linspace(np.log(0.5**0.5/40), np.log(0.5**0.5*40), 10))):
    # We split Q into Q_res for the resonant part and Q_damp for the damping part
    Q_res = (Q**k + (0.5**0.5)**k)**(1/k)
    Q_damp = Q / Q_res
    
    # Compute half of the poles
    eps = (2*Q_res**2 - 1)/(4*Q_res**2 - 1)**0.5
    m = np.arange(0, n//2)
    theta = np.pi/2 * (2*m+1)/n
    ar = np.arcsinh(1/eps)/n
    p = - np.sinh(ar) * np.sin(theta) + 1j * np.cosh(ar) * np.cos(theta)
    w_norm = T_inv(n, 1/(2*Q_res**2 - 1)) # Shift center frequency such that gain match the butterworth one at w = 1
    p /= w_norm
    
    # Do sos and plot response
    H = np.ones_like(s)
    a0, a1, a2 = np.ones(len(p)), -p-np.conj(p), p*np.conj(p)
    a1[0] /= Q_damp
    for i in range(len(p)):
        H *= 1/(a0[i]*s**2/a2[i] + a1[i]*s/a2[i] + 1)
    plt.plot(np.log2(w), 20*np.log10(abs(H)), label = "\\(Q=" + str(round(Q, 2)) + "\\)", c="C" + str(j)) # Chebyshev computed via sos

    # Plot tchebychev magnitude only (without computing each pole), recompute prev values to optim
    g2 = (1 + eps**2) / (1 + (eps**2*T(n, w * w_norm)**2))
    wPeak = np.sqrt(np.real(p[0])**2 + np.imag(p[0])**2)
    QPeak = wPeak / (-2*np.real(p[0]))
    w2 = w*w / (wPeak*wPeak)
    g2 *= 1 + (1 / (QPeak * QPeak) - 2) * w2 + w2 * w2
    g2 /= 1 + (1 / (Q_damp * Q_damp * QPeak * QPeak) - 2) * w2 + w2 * w2
    #plt.plot(np.log2(w), 20*np.log10(g2**0.5), c="C" + str(j)) # Chebyshev computed via magnitude formula

    # Plot butterworth resonant
    p = np.sin(theta) + 1j * np.cos(theta)
    H = np.ones_like(s)
    a0, a1, a2 = np.ones(len(p)), -p-np.conj(p), p*np.conj(p)
    a1[0] /= np.sqrt(2)*Q
    for i in range(len(p)):
        H *= 1/(a0[i]*s**2/a2[i] + a1[i]*s/a2[i] + 1)
    #plt.plot(np.log2(w), 20*np.log10(abs(H)), c="C" + str(j)) # Butterworth computed via sos
    
plt.ylim(-30, 30)
plt.title("Adding underdamping capability")
plt.xlabel('Frequency (octaves)')
plt.ylabel('Magnitude (dB)')
plt.legend()
plt.show()

Python - Console Output

Python
 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

import numpy as np
import matplotlib.pyplot as plt

# Filter parameters
n = 8
w = np.logspace(-5, 5, 50000, base=2)  # Frequency range: 10 octaves (20 Hz to 20 kHz)
s = 1j * w
k = 8  # Exponent for smoothing

# Loop through varying Q values
for j, Q in enumerate(np.exp(np.linspace(np.log(0.5**0.5 / 40), np.log(0.5**0.5 * 40), 10))):
    # Split Q into resonant (Q_res) and damping (Q_damp) parts
    Q_res = (Q**k + (0.5**0.5)**k)**(1/k)
    Q_damp = Q / Q_res

    # Calculate the first half of the poles for the Chebyshev filter
    eps = (2 * Q_res**2 - 1) / (4 * Q_res**2 - 1)**0.5
    m = np.arange(0, n // 2)
    theta = np.pi / 2 * (2 * m + 1) / n
    ar = np.arcsinh(1 / eps) / n
    p = -np.sinh(ar) * np.sin(theta) + 1j * np.cosh(ar) * np.cos(theta)

    # Normalize poles to match Butterworth filter gain at w = 1
    w_norm = T_inv(n, 1 / (2 * Q_res**2 - 1))
    p /= w_norm

    # Compute SOS response and plot its magnitude
    H = np.ones_like(s)
    a0 = np.ones(len(p))
    a1 = -p - np.conj(p)
    a2 = p * np.conj(p)
    a1[0] /= Q_damp  # Adjust first coefficient for damping
    for i in range(len(p)):
        H *= 1 / (a0[i] * s**2 / a2[i] + a1[i] * s / a2[i] + 1)
    plt.plot(np.log2(w), 20 * np.log10(abs(H)), label=f"\\(Q={round(Q, 2)}\\)", c=f"C{j}")

    # # Compute magnitude response using Chebyshev magnitude formula
    # g2 = (1 + eps**2) / (1 + (eps**2 * T(n, w * w_norm)**2))
    # wPeak = np.sqrt(np.real(p[0])**2 + np.imag(p[0])**2)
    # QPeak = wPeak / (-2 * np.real(p[0]))
    # w2 = w * w / (wPeak**2)
    # g2 *= 1 + (1 / (QPeak**2) - 2) * w2 + w2**2
    # g2 /= 1 + (1 / (Q_damp**2 * QPeak**2) - 2) * w2 + w2**2
    # plt.plot(np.log2(w), 20*np.log10(g2**0.5), c="C" + str(j)) # Chebyshev computed via magnitude formula

    # # Plot the Butterworth resonant response
    # p_butter = np.sin(theta) + 1j * np.cos(theta)
    # H_butter = np.ones_like(s)
    # a0_butter = np.ones(len(p_butter))
    # a1_butter = -p_butter - np.conj(p_butter)
    # a2_butter = p_butter * np.conj(p_butter)
    # a1_butter[0] /= np.sqrt(2) * Q # Adjust first coefficient for Butterworth
    # for i in range(len(p_butter)):
    #     H_butter *= 1 / (a0_butter[i] * s**2 / a2_butter[i] + a1_butter[i] * s / a2_butter[i] + 1)
    # plt.plot(np.log2(w), 20*np.log10(abs(H)), c="C" + str(j)) # Butterworth computed via sos

# Plot settings
plt.ylim(-30, 30)
plt.title("Adding Underdamping Capability to Chebyshev Filter")
plt.xlabel('Frequency (octaves)')
plt.ylabel('Magnitude (dB)')
plt.legend()
plt.show()

Python - Console Output

Conclusion

The exploration and implementation of the Musical Chebyshev filter have highlighted its remarkable versatility in sound design, particularly in creating distinctive phaser-style effects. Integrated into UVI Shade as a Multi Resonant filter, it excels in delivering sounds, especially when modulation is applied.

While the filter’s inherent ringing can be pronounced (sometimes overly sharp for certain musical applications), there is an exciting opportunity to develop a new parameter that could mitigate this effect.

Overall, the well-studied Chebyshev filter stands out as a powerful tool for sound designers!