Audio EQ Cookbook without frequency warping

Question

The famous https://www.w3.org/TR/audio-eq-cookbook/ offers a set of [real] biquad filter calculation formulas that generally work fine.

However, when filter's frequency approaches Nyquist frequency the Q (bandwidth) specification of a filter becomes greatly distorted - usually it shrinks a lot (even though the author mentioned that he performed a necessary pre-warping).

I'm in quest for filter formulas that do not have such strong bandwidth distortion. I need peaking/bell, band-pass, low-pass, high-pass, high-shelf, notch filters. I know this can be done as previously I've bought a peaking/bell/band-pass filter formulas with less distortion, but they are still not perfect and I need other filter types.

So, I'm also willing to pay for the solution if the price is right.

Alternatively, if one could point me to an optimization algorithm that works with Z-domain filters that would be great, too. Unfortunately, most usual optimization algorithms do not work well in Z-domain - they can't optimize a set of parameters to match a desired frequency response (probably due to periodic functions used to calculate the frequency response).

Answer 1

here is a quick look at how the 5 degrees of freedom for the parametric EQ can be viewed. it's my take on what Knud Christensen of tc electronic came up with about a decade ago at an AES convention and this patent.

so, forget about the Cookbook (and the issues of Q and bandwidth therein) and consider (in the s-plane) the parametric EQ as the sum of a bandpass filter (with a $Q$ value) in parallel with a wire:

$$ H(s) = (G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + 1 $$

$G_\text{boost} = 10^{\frac{dB}{20}}$ is the gain of the peak (or valley, if $dB<0$). the gain at DC and at Nyquist is 0 dB. that's a 2nd-order IIR and there are 3 independent parameters. 2 more to go. so we next add an overall gain parameter:

$$ H(s) = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + 1 \right) $$

that's 4 knobs to twist. one more knob to add (without raising the filter order) and we will be done with adding more independent parameters.

so what Knud does here is replace that "wire" (that trailing "$1$" in the transfer function) with a prototype shelving filter that must have the same poles, the same $Q$ and $\omega_0$ as the BPF, so that the denominator is the same. the transfer function of that shelf is:

$$ H_\text{shelf}(s) = \frac{R\left(\frac{s}{\omega_0}\right)^2 + \frac{\sqrt{R}}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} $$

where $R \triangleq 10^{\frac{tilt}{20}}$ and $tilt$ is the gain differential of the shelf in dB. this is what offsets the gain at Nyquist to be different than the gain at DC. After bilinear transformation, Nyquist gets boosted by $tilt$ dB and the gain at DC remains unchanged. Like the $dB$ boost parameter, the $tilt$ parameter can be either positive or negative. $G_0 \cdot R$ is the linear gain at Nyquist.

put that all together and you get:

$$ \begin{align} H(s) & = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + H_\text{shelf}(s) \right) \\ & = G_0 \left((G_\text{boost} - 1)\frac{\frac{1}{Q}\frac{s}{\omega_0}}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} + \frac{R\left(\frac{s}{\omega_0}\right)^2 + \frac{\sqrt{R}}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} \right) \\ & = G_0 \frac{R\left(\frac{s}{\omega_0}\right)^2 + (G_\text{boost}+\sqrt{R}-1)\frac{1}{Q}\frac{s}{\omega_0} + 1}{\left(\frac{s}{\omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\omega_0} + 1} \\ \end{align}$$

no matter how you look at that, this has 5 degrees of freedom and those 5 biquad coefficients are fully defined from these 5 parameters. doesn't matter if you map from $s$ to $z$ using the blinear transform or the trapezoidal rule (effectively the same thing) or any other method that does not change the order of the filter. you may have to fudge the definition of $Q$ or bandwidth, you may have to compensate $\omega_0$ and/or $Q$ for frequency warping effects (like you get with the blinear transform), but if you paid big bucks for something that gets you a 2nd-order IIR filter, it doesn't matter if you implement it with any Direct Form or transposed Direct Form or Lattice or Normalized Ladder or Hal Chamberlin's State-Variable or Andrew Simpson's modeling of linear analog with trapezoidal integration, eventually you get to 5 coefficients and they can be mapped to these 5 independent parameters. it's all the same. whether or not you paid money for a license or not. the math is stronger than any claims made by whomever you are licensing from.

Just FYI, i solved where the true peak or valley frequency is when there is a $tilt$ that is not zero. the frequency where the peak or valley has been nudged over by the tilt is:

$$ \omega_\text{peak} \ = \ \omega_0 \ \sqrt{ \frac{Q^2 \left(R - \frac{1}{R}\right)}{G_\text{boost}^2 - R + 2 Q^2 (R - 1)} + \\ \sqrt{\frac{G_\text{boost}^2 - \frac{1}{R} + 2 Q^2 \left(\frac{1}{R} - 1\right)}{G_\text{boost}^2 - R + 2 Q^2 (R - 1)} + \frac{Q^4 \left(R - \frac{1}{R}\right)^2}{\left(G_\text{boost}^2 - R + 2 Q^2 (R - 1)\right)^2} } } $$

you can see that when $tilt = 0$, then $R=1$ and consequently $\omega_\text{peak} = \omega_0$. the peak gain $G_\text{boost}$ might also have to be adjusted a little and that has yet to be worked out. a good first guess would be $G_\text{boost} \leftarrow \frac{G_\text{boost}}{\sqrt{R}}$ or maybe $G_\text{boost} \leftarrow G_\text{boost}-(\sqrt{R}-1)$.

Answer 2

Using optimization methods, we can get a digital filter's frequency response closer to the target analog filter.

In the following experiment, a 6-order bandpass filter is optimized using Adam, an optimization algorithm often used in machine learning. Frequencies above the passband are excluded from the cost function (assigned zero weight). The optimized filter's response becomes higher than target for frequencies very close to Nyquist, but that difference may be offset by the anti-aliasing filter of the signal source (ADC or sample rate converter).

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

import tensorflow as tf

# Number of sections
M = 3

# Sample rate
f_s = 24000

# Passband center frequency
f0 = 9000

# Number of frequencies to compute
N = 2048

section_colors = np.zeros([M, 3])
for k in range(M):
    section_colors[k] = clr.hsv_to_rgb([(k / (M - 1.0)) / 3.0, 0.5, 0.75])

# Get one of BP poles that maps to LP prototype pole.
def lp_to_bp(s, rbw, w0):
    return w0 * (s * rbw / 2 + 1j * np.sqrt(1.0 - np.power(s * rbw / 2, 2)))

# Frequency response

def freq_response(z, b, a):
    p = b[0]
    q = a[0]
    for k in range(1, len(b)):
        p += b[k] * np.power(z, -k)
    for k in range(1, len(a)):
        q += a[k] * np.power(z, -k)
    return p / q

# Absolute value in decibel

def abs_db(h):
    return 20 * np.log10(np.abs(h))

# Poles of analog low-pass prototype

none, S, none = signal.buttap(M)

# Band limits
c = np.power(2, 1 / 12.0)
f_l = f0 / c
f_u = f0 * c

# Analog frequencies in radians
w0 = 2 * np.pi * f0
w_l = 2 * np.pi * f_l
w_u = 2 * np.pi * f_u

# Relative bandwidth
rbw = (w_u - w_l) / w0

jw0 = 2j * np.pi * f0
z0 = np.exp(jw0 / f_s)

# 1. Analog filter parameters

bc, ac = signal.butter(M, [w_l, w_u], btype='bandpass', analog=True)
ww, H_a = signal.freqs(bc, ac, worN=N)
magnH_a = np.abs(H_a)
f = ww / (2 * np.pi)

omega_d = ww / f_s
z = np.exp(1j * ww / f_s)

# 2. Initial filter design

a = np.zeros([M, 3], dtype=np.double)
b = np.zeros([M, 3], dtype=np.double)
hd = np.zeros([M, N], dtype=np.complex)

# Pre-warp the frequencies

w_l_pw = 2 * f_s * np.tan(np.pi * f_l / f_s)
w_u_pw = 2 * f_s * np.tan(np.pi * f_u / f_s)
w_0_pw = np.sqrt(w_l_pw * w_u_pw)

rbw_pw = (w_u_pw - w_l_pw) / w_0_pw

poles_pw = lp_to_bp(S, rbw_pw, w_0_pw)

# Bilinear transform

T = 1.0 / f_s
poles_d = (1.0 + poles_pw * T / 2) / (1.0 - poles_pw * T / 2)

for k in range(M):
    p = poles_d[k]
    b[k], a[k] = signal.zpk2tf([-1, 1], [p, np.conj(p)], 1)

    g0 = freq_response(z0, b[k], a[k])
    g0 = np.abs(g0)
    b[k] /= g0
    none, hd[k] = signal.freqz(b[k], a[k], worN=omega_d)

plt.figure(2)
plt.title("Initial digital filter (bilinear)")

plt.axis([0, f_s / 2, -90, 10])

plt.plot(f, abs_db(H_a), label='Target response', color='gray', linewidth=0.5)

for k in range(M):
    label = "Section %d" % k
    plt.plot(f, abs_db(hd[k]), color=section_colors[k], alpha=0.5, label=label)

# Combined frequency response of initial digital filter

Hd = np.prod(hd, axis=0)
plt.plot(f, abs_db(Hd), 'k', label='Cascaded filter')
plt.legend(loc='upper left')

plt.figure(3)
plt.title("Initial filter - poles and zeros")
plt.axis([-3, 3, -2.25, 2.25])
unitcircle = plt.Circle((0, 0), 1, color='lightgray', fill=False)
ax = plt.gca()
ax.add_artist(unitcircle)

for k in range(M):
    zeros, poles, gain = signal.tf2zpk(b[k], a[k])
    plt.plot(np.real(poles), np.imag(poles), 'x', color=section_colors[k])
    plt.plot(np.real(zeros), np.imag(zeros), 'o', color='none', markeredgecolor=section_colors[k], alpha=0.5)

# Optimizing filter

tH_a = tf.constant(magnH_a, dtype=tf.float32)

# Assign weights

weight = np.zeros(N)
for i in range(N):
    # In the passband or below?
    if (f[i] <= f_u):
        weight[i] = 1.0

tWeight = tf.constant(weight, dtype=tf.float32)
tZ = tf.placeholder(tf.complex64, [1, N])

# Variables to be changed by optimizer
ta = tf.Variable(a)
tb = tf.Variable(b)
ai = a
bi = b

# TF requires matching types for multiplication;
# cast real coefficients to complex
cta = tf.cast(ta, tf.complex64)
ctb = tf.cast(tb, tf.complex64)

xb0 = tf.reshape(ctb[:, 0], [M, 1])
xb1 = tf.reshape(ctb[:, 1], [M, 1])
xb2 = tf.reshape(ctb[:, 2], [M, 1])

xa0 = tf.reshape(cta[:, 0], [M, 1])
xa1 = tf.reshape(cta[:, 1], [M, 1])
xa2 = tf.reshape(cta[:, 2], [M, 1])

# Numerator:   B = b₀z² + b₁z + b₂
tB = tf.matmul(xb0, tf.square(tZ)) + tf.matmul(xb1, tZ) + xb2

# Denominator: A = a₀z² + a₁z + a₂
tA = tf.matmul(xa0, tf.square(tZ)) + tf.matmul(xa1, tZ) + xa2

# Get combined frequency response
tH = tf.reduce_prod(tB / tA, axis=0)

iterations = 2000
learning_rate = 0.0005

# Cost function
cost = tf.reduce_mean(tWeight * tf.squared_difference(tf.abs(tH), tH_a))
optimizer = tf.train.AdamOptimizer(learning_rate).minimize(cost)

zz = np.reshape(z, [1, N])

with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())

    for epoch in range(iterations):
        loss, j = sess.run([optimizer, cost], feed_dict={tZ: zz})
        if (epoch % 100 == 0):
            print("  Cost: ", j)

    b, a = sess.run([tb, ta])

for k in range(M):
    none, hd[k] = signal.freqz(b[k], a[k], worN=omega_d)

plt.figure(4)
plt.title("Optimized digital filter")

plt.axis([0, f_s / 2, -90, 10])

# Draw the band limits
plt.axvline(f_l, color='black', linewidth=0.5, linestyle='--')
plt.axvline(f_u, color='black', linewidth=0.5, linestyle='--')

plt.plot(f, abs_db(H_a), label='Target response', color='gray', linewidth=0.5)

Hd = np.prod(hd, axis=0)
for k in range(M):
    label = "Section %d" % k
    plt.plot(f, abs_db(hd[k]), color=section_colors[k], alpha=0.5, label=label)

magnH_d = np.abs(Hd)
plt.plot(f, abs_db(Hd), 'k', label='Cascaded filter')
plt.legend(loc='upper left')

plt.figure(5)
plt.title("Optimized digital filter - Poles and Zeros")
plt.axis([-3, 3, -2.25, 2.25])
unitcircle = plt.Circle((0, 0), 1, color='lightgray', fill=False)
ax = plt.gca()
ax.add_artist(unitcircle)

for k in range(M):
    zeros, poles, gain = signal.tf2zpk(b[k], a[k])
    plt.plot(np.real(poles), np.imag(poles), 'x', color=section_colors[k])
    plt.plot(np.real(zeros), np.imag(zeros), 'o', color='none', markeredgecolor=section_colors[k], alpha=0.5)

plt.show()

Answer 3

@Jazz, one of the things we learned in electrical engineering is that any order of differential equation can be broken down to a set (or "system") of 1st-order diff eqs. so if trapezoidal integration, with the same "time step" $\Delta t$ is being used consistently for all continuous-time integrals, for a $N$th-order linear ODE, you can bust that up into $N$ first-order differential equations. then consider just one of those 1st-order diff eqs:

again, consider emulating a capacitor. let the sampling period be $T=\frac{1}{f_\text{s}}$ be the same as the "$\Delta t$" used in the trapezoid rule.

$$ i(t) = C \frac{dv}{dt} $$

or

$$ v(t) = \frac{1}{C} \int\limits_{-\infty}^{t} i(u) \ du $$

in the s-domain it's

$$ V(s) = \frac{1}{s} \left( \frac{1}{C} I(s) \right) $$

so trapazoidal integration at discrete times is:

$$\begin{align} v(nT) & = \frac{1}{C} \int\limits_{-\infty}^{nT} i(u) \ du \\ & = \frac{1}{C} \sum\limits_{k=-\infty}^{n} \quad \int\limits_{kT-T}^{kT} i(u) \ du \\ & \approx \frac{1}{C} \sum\limits_{k=-\infty}^{n} \frac{T}{2} \left( i(kT-T) + i(kT) \right) \\ & = v((n-1)T) + \frac{1}{C} \frac{T}{2} \left( i((n-1)T) + i(nT) \right) \\ \end{align}$$

or as discrete-time sample values

$$ v[n] = v[n-1] + \frac{T}{2C} (i[n] + i[n-1]) $$

applying the Z transform

$$ V(z) = z^{-1} V(z) \ + \ \frac{T}{2C} \left(I(z) + z^{-1} I(z) \right) $$

solving for $V$

$$ V(z) = \frac{T}{2} \frac{1 + z^{-1}}{1 - z^{-1}} \left( \frac{1}{C} I(z) \right) $$

looks like we're substituting

$$ \frac{1}{s} \leftarrow \frac{T}{2} \frac{1 + z^{-1}}{1 - z^{-1}} $$

or

$$ s \leftarrow \frac{2}{T} \frac{z-1}{z+1} $$

which is precisely what the bilinear without compensation for frequency warping does.

Answer 4

I've come up with a design for the 10dB peaking EQ. I've chosen 20 filters with center frequencies between 500 Hz and 16 kHz (Fs = 48 kHz). The top plot below is the design according to RBJ's Audio-EQ-Cookbook, which is good but which leads to bandwidth distortion when the center frequencies get closer to Nyquist. The bottom plot is the new design where the filters very closely match the analog prototype filters:

And this is how the new notch filters look like compared to the Cookbook (bandwidth = 4 octaves, highest $f_0=23$ kHz):

The following figure shows a lowpass filter design ($Q=2$, $f_0=16$ kHz, $F_s=48$ kHz). Note that the new design approximates the analog prototype and for this reason it does not perform as a conventional lowpass filter (it has no zero at Nyquist):