Audio EQ Cookbook sin deformación de frecuencia

El famoso http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt ofrece un conjunto de fórmulas de cálculo de filtro biquad [real] que generalmente funcionan bien.

Sin embargo, cuando la frecuencia del filtro se acerca a la frecuencia de Nyquist, la especificación Q (ancho de banda) de un filtro se distorsiona mucho, por lo general, se reduce mucho (aunque el autor mencionó que realizó un pre-warping necesario).

Estoy buscando fórmulas de filtro que no tengan una distorsión de ancho de banda tan fuerte. Necesito filtros de pico / campana, paso de banda, paso bajo, paso alto, estante alto y muesca. Sé que esto se puede hacer ya que anteriormente compré fórmulas de filtro de pico / campana / paso de banda con menos distorsión, pero todavía no son perfectas y necesito otros tipos de filtros.

Entonces, también estoy dispuesto a pagar por la solución si el precio es correcto.

Alternativamente, si uno pudiera señalarme un algoritmo de optimización que funcione con filtros de dominio Z, también sería genial. Desafortunadamente, los algoritmos de optimización más habituales no funcionan bien en el dominio Z: no pueden optimizar un conjunto de parámetros para que coincida con una respuesta de frecuencia deseada (probablemente debido a las funciones periódicas utilizadas para calcular la respuesta de frecuencia).

Aquí hay un vistazo rápido a cómo se pueden ver los 5 grados de libertad para el ecualizador paramétrico. Es mi opinión sobre lo que Knud Christensen de tc electronic ideó hace aproximadamente una década en una convención de AES y esta patente .

así que olvídate del Cookbook (y los problemas de Q y el ancho de banda que contiene) y considera (en el plano s) el EQ paramétrico como la suma de un filtro de paso de banda (con un $Q$ valor) en paralelo con un cable:

$$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}}$ es la ganancia del pico (o valle, si $dB<0$) la ganancia en DC y en Nyquist es de 0 dB. ese es un IIR de segundo orden y hay 3 parámetros independientes. 2 más para ir. así que a continuación agregamos un parámetro de ganancia general:

$$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)$$

eso son 4 perillas para girar. una perilla más para agregar (sin aumentar el orden del filtro) y terminaremos agregando parámetros más independientes.

entonces lo que Knud hace aquí es reemplazar ese "cable" (ese final "$1$"en la función de transferencia) con un prototipo de filtro de estantería que debe tener los mismos polos, los mismos$Q$ y $\omega_0$como el BPF, por lo que el denominador es el mismo. La función de transferencia de ese estante es:

$$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}$$

dónde $R \triangleq 10^{\frac{tilt}{20}}$ y $tilt$es el diferencial de ganancia de la plataforma en dB. Esto es lo que compensa la ganancia en Nyquist para ser diferente de la ganancia en DC. Después de la transformación bilineal, Nyquist se ve impulsado por$tilt$dB y la ganancia en DC permanece sin cambios. Como el$dB$ parámetro de impulso, el $tilt$ El parámetro puede ser positivo o negativo. $G_0 \cdot R$ es la ganancia lineal en Nyquist.

pon todo junto y obtienes:

$$\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 importa cómo lo veas, esto tiene 5 grados de libertad y esos 5 coeficientes biquad están completamente definidos a partir de estos 5 parámetros. no importa si haces un mapa desde$s$ a $z$usando la transformación blinear o la regla trapezoidal (efectivamente la misma cosa) o cualquier otro método que no cambie el orden del filtro. puede que tengas que falsificar la definición de$Q$ o ancho de banda, puede que tenga que compensar $\omega_0$ y / o $Q$para efectos de deformación de frecuencia (como los que se obtienen con la transformación blinear), pero si pagó mucho dinero por algo que le proporciona un filtro IIR de segundo orden, no importa si lo implementa con cualquier Forma Directa o Forma Directa transpuesta o Lattice o Normalized Ladder o la variable de estado de Hal Chamberlin o Andrew Simpson de análogo lineal con integración trapezoidal, eventualmente se obtienen 5 coeficientes y se pueden asignar a estos 5 parámetros independientes. todo es lo mismo. si pagó o no dinero por una licencia o no. la matemática es más sólida que cualquier reclamo hecho por la persona de la que está otorgando la licencia

Para su información, he resuelto, donde el verdadero pico o valle de frecuencia es cuando no es una$tilt$eso no es cero. La frecuencia en la que el pico o el valle ha sido empujado por la inclinación es:

$$\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} } }$$

puedes ver eso cuando $tilt = 0$, entonces $R=1$ y consecuentemente $\omega_\text{peak} = \omega_0$. la ganancia máxima$G_\text{boost}$también podría tener que ajustarse un poco y eso aún no se ha resuelto. una buena primera suposición sería$G_\text{boost} \leftarrow \frac{G_\text{boost}}{\sqrt{R}}$ o tal vez $G_\text{boost} \leftarrow G_\text{boost}-(\sqrt{R}-1)$.

@Jazz, una de las cosas que aprendimos en ingeniería eléctrica es que cualquier orden de ecuación diferencial se puede dividir en un conjunto (o "sistema") de ecuaciones diferenciales de primer orden. así que si la integración trapezoidal, con el mismo "paso de tiempo"$\Delta t$está siendo usado consistentemente para todas las integrales de tiempo continuo, para un$N$ODE lineal de orden th, puede reventar eso en $N$ecuaciones diferenciales de primer orden. entonces considere solo una de esas ecuaciones diferenciales de primer orden:

nuevamente, considere emular un condensador. deje que el período de muestreo sea$T=\frac{1}{f_\text{s}}$ ser lo mismo que el "$\Delta t$"utilizado en la regla trapezoidal.

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

o

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

en el dominio s es

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

entonces la integración trapazoidal en momentos discretos es:

$$\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}$$

o como valores de muestra de tiempo discreto

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

aplicando la transformación Z

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

resolviendo para $V$

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

parece que estamos sustituyendo

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

o

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

que es precisamente lo que hace el bilineal sin compensación por deformación de frecuencia.

Usando métodos de optimización, podemos acercar la respuesta de frecuencia de un filtro digital al filtro analógico objetivo.

En el siguiente experimento, un filtro de paso de banda de 6 órdenes se optimiza usando Adam, un algoritmo de optimización que se usa a menudo en el aprendizaje automático. Las frecuencias por encima de la banda de paso están excluidas de la función de costo (peso cero asignado). La respuesta del filtro optimizado se vuelve más alta que el objetivo para frecuencias muy cercanas a Nyquist, pero esa diferencia puede ser compensada por el filtro anti-aliasing de la fuente de señal (ADC o convertidor de frecuencia de muestreo).

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()

Se me ocurrió un diseño para el ecualizador de pico de 10dB. Elegí 20 filtros con frecuencias centrales entre 500 Hz y 16 kHz (Fs = 48 kHz). La trama superior a continuación es el diseño de acuerdo con Audio-EQ-Cookbook de RBJ , que es bueno pero que conduce a una distorsión del ancho de banda cuando las frecuencias centrales se acercan a Nyquist. La trama inferior es el nuevo diseño donde los filtros coinciden muy estrechamente con los filtros prototipo analógicos:

Y así es como se ven los nuevos filtros de muesca en comparación con el Cookbook (ancho de banda = 4 octavas, el más alto $f_0=23$ kHz):

La siguiente figura muestra un diseño de filtro de paso bajo ($Q=2$, $f_0=16$ kHz $F_s=48$kHz). Tenga en cuenta que el nuevo diseño se aproxima al prototipo analógico y, por esta razón, no funciona como un filtro de paso bajo convencional (no tiene cero en Nyquist):