A New Method for Eliminating Bilinear Transform Stretching Artefacts

Table Of Contents

Digitalizing Filters: Addressing the Stretching Artifact

Analog filters are frequently digitized using the Bilinear Transform (BLT). However, this process stretches the frequency response in the high-frequency range, particularly near the Nyquist limit. Since audio filters are ideally independent of the sampling rate, this artifact is often undesirable, and we aim to minimize it for a response closer to the analog version.

Several solutions to this problem have been proposed, including:

• Cutoff Frequency Pre-Warping: This is the most common approach, as detailed in the RBJ Cookbook¹. It involves pre-warping the cutoff and bandwidth parameters before applying the BLT. However, this method leads to undefined behavior when filters are stacked² and still leaves some distortion in the magnitude response.

• Oversampling: Increasing the sampling rate can reduce the stretching artifact but introduces other issues, such as increased CPU usage, higher latency, and potentially new artifacts.

• Impulse Invariance³: This method causes aliasing and is not suitable for all filter types.

• Matched Z-Transform⁴: This technique is unpredictable and suffers from discontinuity when poles or zeros fall outside the Nyquist frequency.

• Filter-Specific Tweaks: Authors like S. J. Orfanidis⁵, M. Masberg⁶, and M. Vicanek⁷ have developed methods that focus on specific filter types, fine-tuning coefficients to minimize distortion. However, these approaches generally lack generality, apply only to certain filter designs, and may sacrifice desirable mathematical properties.

Proposed Solution for Reducing Stretching Artifacts: The Magnitude-Matching Transform

This proposed solution significantly reduces stretching artifacts due to the classic Bilinear Transformation (BLT) while maintaining low CPU resource requirements. The approach involves modifying the coefficients of the analog filter immediately prior to applying the BLT, using the transformations outlined in the following subsection with the coefficient \(\alpha = 0.15\).

First-Order Filters

For a first-order transfer function defined as \(H_1(s) = \frac{b_0 s + b_1}{a_0 s + a_1},\) the transformation of the coefficients can be expressed as follows:

$$ \boxed{ \begin{aligned} b_0 &\leftarrow \sqrt{b_0^{2}+\alpha b_1^{2}} \newline b_1 &\leftarrow b_1 \newline a_0 &\leftarrow \sqrt{a_0^{2}+\alpha a_1^{2}} \newline a_1 &\leftarrow a_1 \end{aligned}} $$

Second-Order Filters

For a second-order transfer function defined as \(H_2(s) = \frac{b_0 s^2 + b_1 s + b_2}{a_0 s^2 + a_1 s + a_2}\) the coefficient transformation is given by:

$$ \boxed{ \begin{aligned} b_0 &\leftarrow \sqrt{\alpha b_1^{2}+\left(\alpha b_2-b_0\right)^{2}} \newline b_1 &\leftarrow \sqrt{2b_2\left(\alpha b_2-b_0+\sqrt{\alpha b_1^{2}+\left(\alpha b_2-b_0\right)^{2}}\right)+b_1^{2}} \newline b_2 &\leftarrow b_2 \newline a_0 &\leftarrow \sqrt{\alpha a_1^{2}+\left(\alpha a_2-a_0\right)^{2}} \newline a_1 &\leftarrow \sqrt{2a_2\left(\alpha a_2-a_0+\sqrt{\alpha a_1^{2}+\left(\alpha a_2-a_0\right)^{2}}\right)+a_1^{2}} \newline a_2 &\leftarrow a_2 \end{aligned}} $$

Expected Results

The following examples illustrate how the magnitude response of filters utilizing our method closely resembles that of the original analog filter, especially when compared to the traditional Bilinear Transformation (BLT). This comparison highlights the efficacy of our approach. By using our method at a sampling rate of 44100 Hz, the minor residual artifacts typically fall within an inaudible range.

 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

import numpy as np
from numpy.fft import fft, ifft, rfft, irfft, fftshift, ifftshift
from scipy import signal
import matplotlib.pyplot as plt
from matplotlib.pyplot import plot

def prewarp_second_order(a, c = 0.15):
    a0, a1, a2 = a
    a0_ = np.sqrt(a1 * a1 * c + (a0 - a2 * c)**2)
    A1 = (2 * a2 * a2 * c + a1 * a1 - 2 * a2 * a0) / 2 # B / 2
    A3 = np.sqrt(0j + (a1 * a1 * (a1 * a1 - 4 * a0 * a2))) / 2
    a1_ = np.real((np.sqrt(A1 - A3) + np.sqrt(A1 + A3)))
    return a0_, a1_, a2

def prewarp_first_order(a, c = 0.15):
    return (a[1]**2*c +a[0]**2)**0.5, a[1]

w = np.logspace(np.log2(0.1), np.log2(np.pi), 10000, base = 2)
coeffs, names = list(), list()
names.append("Resonant Low Pass")
coeffs.append(([0, 0, 1], [1, 0.2, 1]))
names.append("Resonant High Pass")
coeffs.append(([1, 0, 0], [1, 0.2, 1]))
names.append("Peak Filter")
coeffs.append(([1, 1, 1], [1, 0.2, 1]))
for (b, a), name in zip(coeffs, names):
    _, h = signal.freqz(*signal.bilinear(b, a), w)
    plt.plot(np.log2(w[h!=0]), 20*np.log10(np.abs(h[h!=0])), label="Bilinear Transform")
    _, h = signal.freqz(*signal.bilinear(prewarp_second_order(b), prewarp_second_order(a)), w)
    plt.plot(np.log2(w), 20*np.log10(np.abs(h)), label="Magnitude Matching Transform")
    _, h = signal.freqs(b, a, w)
    plt.plot(np.log2(w), 20*np.log10(np.abs(h)), label="Analog (ideal)")
    plt.axvline(x=np.log2(np.pi), ls='--', linewidth=0.7, c="#BBC3CF")
    plt.ylim(-30, 30)
    plt.title("Magnitude of a " + name)
    plt.xlabel('Frequency (octaves)')
    plt.ylabel('Magnitude (dB)')
    plt.legend()
    plt.show()

Comparative Analysis of Filter Methods

The table below compares our proposed approach with alternative methods for reducing stretching artifacts:

Our method is applicable to any filter decomposable into second-order sections (SOS) and offers an excellent magnitude response that closely mirrors the analog version, representing a significant improvement over other methods.

How does it work ?

Working with magnitude

The idea behind this method is to work with the magnitude response instead of the full transfer function, allowing us greater flexibility to distort the magnitude response without concerning ourselves with the phase information.

It is possible to transform a transfer function \(H(s)\) into its magnitude response \(G(\omega)\) and to perform the inverse transformation from \(G(\omega)\) back to \(H(s)\). \(G(\sqrt{x})^2\) must be a rational function for this inverse transformation to work. As a result, we can pre-distort the frequency axis using a Möbius transformation on \(G(\sqrt{x})^2\) such that it compensates for the artifact introduced during the digitalization process. The Möbius transformation offers greater degrees of freedom than conventional cutoff frequency pre-warping, resulting in a more accurate magnitude response compared to the classic pre-warping method.

Mathematical Representation of Frequency Distortion

As we have seen, the bilinear transform, which is used to convert analog filters into digital form, introduces unwanted distortion in the frequency axis. The relationship between the analog and distorted digital frequency axes can be expressed as follows:

$$\omega_d = f_{\text{warp}}(\omega_a)$$ $$f_{\text{warp}}(\omega) = 2\left(\arctan\left(\frac{\omega}{2}\right)\right)$$

Correcting Frequency Distortion Artifacts

In an ideal scenario, we would directly apply \(f_{\text{\text{warp}}}(\omega)\) to pre-distort the frequency axis. However, for the reverse transformation to work, \(G(\sqrt{x})^2\) must remain rational, which requires an approximation.

We approximate \(f_{\text{\text{warp}}}(\omega)\) with the function:

$$\boxed{\tilde{f}_{\text{warp}}(\omega) = \frac{\omega}{\sqrt{\alpha\omega^{2}+1}}, \alpha = 0.15}$$

Application to a first order filter

The transfert function of a fos filter is given by: $$H(s) = \frac{b_0s+b_1}{a_0s+a_1}$$

Its magnitude response is expressed as: $$G(w) = |H(j\omega)| = \sqrt{\frac{b_0^{2}\omega^{2}+b_1^{2}}{a_0^{2}\omega^{2}+a_1^{2}}}$$

Its warped magnitude is therefor: $$G_{\text{warped}}\left(w\right)=G\left(\tilde{f}_{\text{warp}}\left(w\right)\right)=\sqrt{\frac{\left(b_0^{2}+\alpha b_1^{2}\right)w^{2}+b_1^{2}}{\left(a_0^{2}+\alpha a_1^{2}\right)w^{2}+a_1^{2}}}$$

By identifying the terms, we obtain: $$\boxed{H_{\text{warped}}(s) = \frac{\sqrt{b_0^{2}+\alpha b_1^{2}}s+b_1}{\sqrt{a_0^{2}+\alpha a_1^{2}}s+a_1}}$$

Application to a second order filter

The transfert function of a biquad filter is given by: $$H(s) = \frac{b_0s²+b_1s+b_2}{a_0s²+a_1s+a_2}$$

Its magnitude response is expressed as: $$G_{\text{warped}}\left(\omega\right) = |H(j\omega)| = \sqrt{\frac{b_0^{2}w^{4}+\left(b_1^{2}-2b_0a_2\right)w^{2}+b_2^{2}}{a_0^{2}w^{4}+\left(a_1^{2}-2a_0a_2\right)w^{2}+a_2^{2}}}$$

Its warped magnitude is therefor: $$G_{\text{warped}}\left(\omega\right)=G\left(\tilde{f}_{\text{warp}}\left(w\right)\right)=\sqrt{\frac{\left(\alpha b_1^{2}+\left(\alpha b_2-b_0\right)^{2}\right)\omega^{4}+\left(2b_2\left(\alpha b_2-b_0\right)+b_1^{2}\right)\omega^{2}+b_2^{2}}{\left(\alpha a_1^{2}+\left(\alpha a_2-a_0\right)^{2}\right)\omega^{4}+\left(2a_2\left(\alpha a_2-a_0\right)+a_1^{2}\right)\omega^{2}+a_2^{2}}}$$

Let be \(c_{k}\) and \(d_{k}\) the coefficients of the new transfert function such its magnitude is : $$G_{\text{warped}}\left(\omega\right)=|H_{\text{warped}}(j\omega)|=\sqrt{\frac{d_{0}^{2}w^{4}+\left(d_{1}^{2}-2d_{0}c_{2}\right)w^{2}+d_{2}^{2}}{c_{0}^{2}w^{4}+\left(c_{1}^{2}-2c_{0}c_{2}\right)w^{2}+c_{2}^{2}}}$$

Then by identification we get : $$ \boxed{ \begin{aligned} H_{\text{warped}}(s) &= \frac{d_0s²+d_1s+d_2}{c_0s²+c_1s+c_2} \newline d_{0} &= \sqrt{\alpha b_1^{2}+\left(\alpha b_2-b_0\right)^{2}} \newline d_{1} &= \sqrt{2b_2\left(\alpha b_2-b_0+d_{0}\right)+b_1^{2}} \newline d_{2} &= b_2 \newline c_{0} &= \sqrt{\alpha a_1^{2}+\left(\alpha a_2-a_0\right)^{2}} \newline c_{1} &= \sqrt{2a_2\left(\alpha a_2-a_0+c_{0}\right)+a_1^{2}} \newline c_{2} &= a_2 \end{aligned}} $$

Poles zeros interpretation

The Magnitude-Matching transform remaps each pole and zero individually according to:

$$z_k \leftarrow \frac{z_k}{\sqrt{\alpha z_k^2 + 1}}, \quad \alpha z_k^2 + 1 \neq 0$$
$$p_k \leftarrow \frac{p_k}{\sqrt{\alpha p_k^2 + 1}}, \quad \alpha p_k^2 + 1 \neq 0$$

In the case where the number of poles was not matching the number of zeros, we will add those missing poles/zeros in the output filter, for which the value are going to be: $$z, p = \frac{1}{\sqrt{\alpha}}$$

This can be understood as missing poles and zeros being interpreted as poles and zeros near infinity. Therefore, we just need to compute limits to resolve issues related to singularities.

At the opposite, in the case where \(\alpha z_k^2+1 = 0\) or \(\alpha p_k^2+1 = 0\), we will remove the given pole/zero from our new filter.

The global gain \(k\) relative to a \(z\), \(p\), \(k\) decomposition is not affected by the warping.

The code below demonstrates how we can remap poles and zeros in a high-order filter.

 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

def prewarp_poles_zeros(z):
    return z / np.sqrt(0j + 1 + 0.15 * z**2)

w = np.logspace(np.log2(0.1), np.log2(np.pi), 10000, base = 2)
# Taking cheby2 to get len(p)==len(z)
z, p, k = signal.cheby2(8, 40, 1, 'low', analog=True, output = "zpk")

b, a = signal.bilinear(*signal.zpk2tf(z, p, k))
_, h = signal.freqz(b, a, w)
plt.plot(np.log2(w[h!=0]), 20*np.log10(np.abs(h[h!=0])), label="Bilinear Transform")

b, a = signal.bilinear(*signal.zpk2tf(prewarp_poles_zeros(z), prewarp_poles_zeros(p), k))
_, h = signal.freqz(b, a, w)
plt.plot(np.log2(w[h!=0]), 20*np.log10(np.abs(h[h!=0])), label="Magnitude-Matching Transform")

b, a = signal.zpk2tf(z, p, k)
_, h = signal.freqs(b, a, w)
plt.plot(np.log2(w[h!=0]), 20*np.log10(np.abs(h[h!=0])), label="Analog (ideal)")

plt.title("Magnitude of a 4-order Chebyshev Low Pass filter")
plt.xlabel('Frequency (octaves)')
plt.ylabel('Magnitude (dB)')

plt.ylim(-100, 10)
plt.legend()
plt.show()

Increasing the Approximation Order

To achieve a more accurate warping approximation, we can increase the order of the rational function \(\tilde{f}_{\text{\text{warp}}}(\sqrt{x})^2\) by a factor of \(n\). This enhancement provides a closer match to the desired warping function, resulting in improved digital filter performance with respect to the analog frequency response. However, increasing the order also raises the order of the resulting digital filter, which will now require \(n\) times more CPU resources.

For most practical applications, the first-order approximation has shown excellent performance, making it sufficient in most cases. However, higher-order approximations may be useful for highly demanding applications where minimal warping error is critical. In such cases, optimal coefficients for the rational function can be calculated using linear regression techniques to closely approximate the target function. Once the desired magnitude response is achieved, the same process applied in the first-order approximation can be used to determine the corresponding transfer function.

Conclusion

This work introduces a flexible and efficient approach to minimizing the warping artifacts associated with the bilinear transform, focusing on the magnitude response and employing a Möbius transformation for frequency adjustment. We have named this digitalization method the Magnitude-Matching Transform, as it focuses exclusively on aligning the magnitude response without regard for phase information. This method stands as a computationally efficient alternative to traditional pre-warping techniques, providing a closer match to the analog frequency response while maintaining low CPU demands. Given its simplicity and versatility, this approach offers audio developers a go-to solution for designing digital filters with analog-consistent characteristics.