Audio EQ Cookbook sem distorção de frequência

O famoso http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt oferece um conjunto de fórmulas [reais] de cálculo de filtro de biquad que geralmente funcionam bem.

No entanto, quando a frequência do filtro se aproxima da frequência de Nyquist, a especificação Q (largura de banda) de um filtro fica muito distorcida - geralmente diminui muito (embora o autor tenha mencionado que ele realizou um pré-empenamento necessário).

Estou em busca de fórmulas de filtro que não tenham uma distorção tão forte na largura de banda. Preciso de filtros de pico / sino, passa-banda, passa-baixo, passa-alto, prateleira alta e entalhe. Sei que isso pode ser feito como anteriormente, comprei uma fórmula de filtro de pico / sino / passa-banda com menos distorção, mas elas ainda não são perfeitas e preciso de outros tipos de filtro.

Então, também estou disposto a pagar pela solução, se o preço for justo.

Como alternativa, se alguém pudesse me indicar um algoritmo de otimização que funcione com filtros do domínio Z, também seria ótimo. Infelizmente, os algoritmos de otimização mais comuns não funcionam bem no domínio Z - eles não podem otimizar um conjunto de parâmetros para corresponder à resposta de frequência desejada (provavelmente devido às funções periódicas usadas para calcular a resposta de frequência).

Aqui está uma rápida olhada em como os 5 graus de liberdade para o EQ paramétrico podem ser vistos. é minha opinião sobre o que Knud Christensen da tc electronic inventou cerca de uma década atrás em uma convenção da AES e nesta patente .

portanto, esqueça o Cookbook (e os problemas de Q e largura de banda) e considere (no plano s) o EQ paramétrico como a soma de um filtro passa-banda (com um valor ) em paralelo com um fio:$Q$

$$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}}$ é o ganho do pico (ou vale, se ). o ganho em DC e em Nyquist é 0 dB. esse é um IIR de segunda ordem e existem 3 parâmetros independentes. Mais 2 para ir. então adicionamos um parâmetro de ganho geral:$dB<0$

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

são 4 botões para girar. mais um botão para adicionar (sem aumentar a ordem dos filtros) e concluiremos a adição de mais parâmetros independentes.

então, o que Knud faz aqui é substituir esse "fio" (aquele que segue " " na função de transferência) por um protótipo de filtro de prateleira que deve ter os mesmos pólos, os mesmos e que o BPF, para que o denominador seja o mesmo. a função de transferência dessa prateleira é:$1$$Q$$\omega_0$

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

onde e é o diferencial de ganho da prateleira em dB. é isso que compensa o ganho em Nyquist para ser diferente do ganho em DC. Após a transformação bilinear, Nyquist é impulsionado pela dB e o ganho no DC permanece inalterado. Como o parâmetro boost, o parâmetro pode ser positivo ou negativo. é o ganho linear em Nyquist.$R \triangleq 10^{\frac{tilt}{20}}$$tilt$$tilt$$dB$$tilt$$G_0 \cdot R$

junte tudo isso e você obtém:

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

não importa como você veja isso, isso tem 5 graus de liberdade e esses 5 coeficientes de biquad são totalmente definidos a partir desses 5 parâmetros. não importa se você mapeia de para usando a transformação blinear ou a regra trapezoidal (efetivamente a mesma coisa) ou qualquer outro método que não altere a ordem do filtro. pode ser necessário alterar a definição de ou largura de banda, compensar e / ou$s$$z$$Q$$\omega_0$$Q$para efeitos de distorção de frequência (como você obtém com a transformação blinear), mas se você pagou muito dinheiro por algo que recebe um filtro IIR de segunda ordem, não importa se você o implementa com qualquer Formulário direto ou Formulário direto transposto ou No modelo Lattice ou Normalized Ladder ou na variável de estado de Hal Chamberlin ou Andrew Simpson, de analógico linear com integração trapezoidal, você obtém 5 coeficientes e eles podem ser mapeados para esses 5 parâmetros independentes. é tudo a mesma coisa. se você pagou ou não por uma licença ou não. a matemática é mais forte do que qualquer reivindicação feita por quem você está licenciando.

Apenas FYI, eu resolvido onde o verdadeiro pico ou vale frequência é quando não é uma que não é zero. a frequência em que o pico ou vale foi empurrado pela inclinação é:$tilt$

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

você pode ver que, quando , então e, consequentemente, . o ganho de pico também pode precisar ser ajustado um pouco e isso ainda precisa ser resolvido. um bom palpite seria ou talvez .$tilt = 0$$R=1$$\omega_\text{peak} = \omega_0$$G_\text{boost}$$G_\text{boost} \leftarrow \frac{G_\text{boost}}{\sqrt{R}}$$G_\text{boost} \leftarrow G_\text{boost}-(\sqrt{R}-1)$

@Jazz, uma das coisas que aprendemos em engenharia elétrica é que qualquer ordem de equação diferencial pode ser dividida em um conjunto (ou "sistema") de equações diferenciais de primeira ordem. por isso, se a integração trapezoidal, com o mesmo "passo de tempo" está sendo usado de forma consistente para todos os integrais de tempo contínuo, para uma º de ordem linear ODE, você pode busto que até em equações diferenciais de primeira ordem. considere apenas um desses eqs diferenciais de primeira ordem:$\Delta t$$N$$N$

novamente, considere emular um capacitor. permita que o período de amostragem seja seja o mesmo que " " usado na regra do trapézio.$T=\frac{1}{f_\text{s}}$$\Delta t$

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

ou

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

no domínio s é

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

então a integração trapazoidal em momentos discretos é:

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

ou como valores de amostra em tempo discreto

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

aplicando a transformação Z

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

resolvendo para$V$

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

parece que estamos substituindo

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

ou

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

que é precisamente o que o bilinear sem compensação pelo desvio de frequência faz.

Usando métodos de otimização, podemos obter a resposta de frequência de um filtro digital mais próxima do filtro analógico de destino.

No experimento a seguir, um filtro passa-banda de 6 ordens é otimizado usando Adam, um algoritmo de otimização frequentemente usado em aprendizado de máquina. As frequências acima da banda passante são excluídas da função de custo (peso zero atribuído). A resposta do filtro otimizado se torna mais alta que o alvo para frequências muito próximas de Nyquist, mas essa diferença pode ser compensada pelo filtro anti-aliasing da fonte de sinal (ADC ou conversor de taxa de amostragem).

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

Eu criei um design para o EQ de 10dB. Eu escolhi 20 filtros com frequências centrais entre 500 Hz e 16 kHz (Fs = 48 kHz). O gráfico superior abaixo é o design de acordo com o Audio-EQ-Cookbook da RBJ , que é bom, mas leva à distorção da largura de banda quando as frequências centrais se aproximam de Nyquist. O gráfico inferior é o novo design em que os filtros correspondem muito aos filtros de protótipo analógico:

E é assim que os novos filtros de entalhe se parecem com o Cookbook (largura de banda = 4 oitavas, maior $f_0=23$ kHz):

A figura a seguir mostra um design de filtro passa-baixo ($Q=2$, $f_0=16$ kHz, $F_s=48$kHz). Observe que o novo design se aproxima do protótipo analógico e, por esse motivo, ele não funciona como um filtro passa-baixo convencional (não tem zero em Nyquist):