How to compute parameter values (frequency, gain, Q-value) to EQ a loudspeaker

[ironhorse128]

I would like to be able to compute parameter values (frequency, gain, Q-value) to EQ a loudspeaker or headphone. I guess it is a more or less complex optimization problem given the frequency response of the loudspeaker/ headphone and a target curve.

I know that there are some programs (e.g. REW) that lets you do that. I would like to read up on the math and algorithms behind these. I did not find suitable resources while my google search, so I thought I might ask ASR.

My final target would be to program a web-tool that lets you upload frequency response to compute a suitable EQ configuration (frequency, gain, Q-value).

So if anyone knows a good start to read up on this topic (books, web-resources, open-source projects, programming libraries,…), please let me know.

 

[NTK]

This is the form of the "bell-shaped" filter transfer function I usually use. (Note that the gain given in the equation is NOT in dB.)

The response of the transfer function Ĥ(s) at a given frequency can be computed by substituting "s" with "j*2π*f", where "j" is square root of -1, "π" is 3.1415926..., and "f" is frequency in Hz. The result is a complex number. Its magnitude is the amplitude, and angle (argument) is the phase. When cascading multiple filters, just multiply the complex responses of each filter to get the combined frequency response.

Here is some python code to calculate and plot the FR of 2 of these filters, and when they combined.

Spoiler: Code

import numpy as np
import matplotlib.pyplot as plt

def peq_tf(gain, bw, f0):
    """
      Return transfer function of a 2nd order bell-shaped filter as the coefficients
      of the numerator and denominator polynomials.
     
      Call signature
        num, den = peq_tf(gain, bw, freq)
      Inputs
        'gain': amplitude at center frequency, 1 = unity gain, >1 = boost, <1 = cut
        'bw': bandwidth in octaves
        'f0': center frequency in Hz
      Outputs
        tuple of ('num', 'den')
        'num': transfer function numerator polynomial coefficients, degree 2 first
        'den': transfer function denominator polynomial coefficients, degree 2 first
    """
    w = 2 * np.pi * f0
    num = np.array([1, 2*np.sqrt(gain)*np.sinh(np.log(2)*bw/2)*w, w**2])
    den = np.array([1, (2/np.sqrt(gain))*np.sinh(np.log(2)*bw/2)*w, w**2])
    return num, den

def peq_fr(tf_coeff, freqs):
    """
      Compute complex frequency response of a 2nd order filter.
     
      Call signature
        fr = peq_fr(tf_coeff, freqs)
      Inputs
        'tf_coeff': tuple of (num, den) specifying the polynomial coefficients of the
                    transfer function, degree 2 first
        'freqs': frequencies at which the responses are computed, in Hz
      Outputs
        The complex frequency responses of the filter at 'freqs'
    """
    j_ws = 2j * np.pi * freqs
    num, den = tf_coeff
    return (num[0]*j_ws**2 + num[1]*j_ws + num[2]) / (den[0]*j_ws**2 + den[1]*j_ws + den[2])

# PEQ filter 1 -- center freq 2000 Hz, bandwidth 1.0, gain 2.0
filt1_params = [2.0, 1.0, 2000]
filt1_tf = peq_tf(*filt1_params)

# PEQ filter 2 -- center freq 200 Hz, bandwidth 3.0, gain 0.5
filt2_params = [0.5, 3.0, 200]
filt2_tf = peq_tf(*filt2_params)

# The frequencies at which we want to compute the filter responses
freqs = np.logspace(np.log10(20), np.log10(20480), 501)

# Compute the complex frequency response of each filter and their combined response
filt1_fr = peq_fr(filt1_tf, freqs)
filt2_fr = peq_fr(filt2_tf, freqs)
filt_combined_fr = filt1_fr * filt2_fr

# Plot the magnitudes of the complex filter responses
fig, axs = plt.subplots(3, 1, figsize=(5, 8), sharex=True)
axs[0].semilogx(freqs, np.abs(filt1_fr), c='C1', label='Filter 1')
axs[1].semilogx(freqs, np.abs(filt2_fr), c='C2', label='Filter 2')
axs[2].semilogx(freqs, np.abs(filt_combined_fr), c='C3', label='Combined')
for ax in axs:
    ax.legend(loc='upper right')
    ax.grid(which='both', axis='both')
    ax.set_ylim(0, 2.5)
fig.tight_layout()
plt.show()

Plot output

I have no real experience on how to find the gains, band-widths (or Q) and center frequencies of a series of these filters to fit the response to a target response. I guess the naïve approach will be to simply feed the problem into a nonlinear regression solver, using some reasonable initial guesses, and with no guarantee of success. There are probably better methods.

 

[ironhorse128]

@NTK thanks a lot. That is a really good starting point .

 

[andyc56]

Someone asked John M. what algorithm he used in REW to do the filter computations. I can't find that post now, but he pointed the user to this article.

 

[NTK]

Sorry to necro this thread. I was approached by a new member on how to compute the frequency response of low and high shelving filters (similar to what I did for peaking filters in post #2).

Below are the transfer functions I modified from the ones provided in the link to "Audio-EQ-Cookbook.txt" given in:

RBJ Audio-EQ-Cookbook — Musicdsp.org documentation

And my Python code to test them out:

Spoiler: Python Code and the plots produced

import numpy as np
import matplotlib.pyplot as plt

def bw2Q(bw):
    """
      Return Q given the bandwidth (bw) in octaves
      Reference: RaneNote 167 (PDF page 5)
                 https://www.ranecommercial.com/legacy/note167.html
    """
    return np.sqrt(2**bw)/(2**bw - 1)

def Q2bw(Q):
    """
      Return the bandwidth (bw) in octaves given Q
      Reference: RaneNote 167 (PDF page 5)
                 https://www.ranecommercial.com/legacy/note167.html
    """
    return (2/np.log(2))*np.log(0.5*(1/Q + np.sqrt(1/Q**2 + 4)))

def low_shelf_tf(gain, bw, f0):
    """
      Return transfer function of a 2nd order low shelving filter as the coefficients
      of the numerator and denominator polynomials.
      Reference: See link to "Audio-EQ-Cookbook.txt" in
                 https://www.musicdsp.org/en/latest/Filters/197-rbj-audio-eq-cookbook.html
     
      Call signature
        num, den = low_shelf_tf(gain, bw, freq)
      Inputs
        'gain': amplitude at DC, 1 = unity gain, >1 = boost, <1 = cut
        'bw': bandwidth in octaves
        'f0': conrner frequency in Hz
      Outputs
        tuple of ('num', 'den')
        'num': transfer function numerator polynomial coefficients, degree 2 first
        'den': transfer function denominator polynomial coefficients, degree 2 first
    """
    Q = bw2Q(bw)
    A = np.sqrt(gain)
    w = 2 * np.pi * f0
    num = A * np.array([1, (np.sqrt(A)/Q)*w, A*w**2])
    den = np.array([A, (np.sqrt(A)/Q)*w, w**2])
    return num, den

def high_shelf_tf(gain, bw, f0):
    """
      Return transfer function of a 2nd order high shelving filter as the coefficients
      of the numerator and denominator polynomials.
      Reference: See link to "Audio-EQ-Cookbook.txt" in
                 https://www.musicdsp.org/en/latest/Filters/197-rbj-audio-eq-cookbook.html
     
      Call signature
        num, den = high_shelf_tf(gain, bw, freq)
      Inputs
        'gain': amplitude at infinite frequency, 1 = unity gain, >1 = boost, <1 = cut
        'bw': bandwidth in octaves
        'f0': conrner frequency in Hz
      Outputs
        tuple of ('num', 'den')
        'num': transfer function numerator polynomial coefficients, degree 2 first
        'den': transfer function denominator polynomial coefficients, degree 2 first
    """
    Q = bw2Q(bw)
    A = np.sqrt(gain)
    w = 2 * np.pi * f0
    num = A * np.array([A, (np.sqrt(A)/Q)*w, w**2])
    den = np.array([1, (np.sqrt(A)/Q)*w, A*w**2])
    return num, den

def peq_fr(tf_coeff, freqs):
    """
      Compute complex frequency response of a 2nd order filter.
     
      Call signature
        fr = peq_fr(tf_coeff, freqs)
      Inputs
        'tf_coeff': tuple of (num, den) specifying the polynomial coefficients of the
                    transfer function, degree 2 first
        'freqs': frequencies at which the responses are computed, in Hz
      Outputs
        The complex frequency responses of the filter at 'freqs'
    """
    j_ws = 2j * np.pi * freqs
    num, den = tf_coeff
    return (num[0]*j_ws**2 + num[1]*j_ws + num[2]) / (den[0]*j_ws**2 + den[1]*j_ws + den[2])

# Filter 1 -- low shelf @ 200 Hz, bandwidth 1.8997 (Q = sqrt(0.5)), gain 2.0
filt1_params = [2.0, 1.8997, 200]
filt1_tf = low_shelf_tf(*filt1_params)

# Filter 2 -- high shelf @ 2000 Hz, bandwidth 3.0 (Q = 0.40406), gain 0.5
filt2_params = [0.5, 3.0, 2000]
filt2_tf = high_shelf_tf(*filt2_params)

# The frequencies at which we want to compute the filter responses
freqs = np.logspace(np.log10(20), np.log10(20480), 501)

# Compute the complex frequency response of each filter and their combined response
filt1_fr = peq_fr(filt1_tf, freqs)
filt2_fr = peq_fr(filt2_tf, freqs)
filt_combined_fr = filt1_fr * filt2_fr

# Plot the magnitudes of the complex filter responses
fig, axs = plt.subplots(3, 1, figsize=(5, 8), sharex=True)
axs[0].semilogx(freqs, np.abs(filt1_fr), c='C1', label='Low Shelf')
axs[1].semilogx(freqs, np.abs(filt2_fr), c='C2', label='High Shelf')
axs[2].semilogx(freqs, np.abs(filt_combined_fr), c='C3', label='Combined')
for ax in axs:
    ax.legend(loc='upper right')
    ax.grid(which='both', axis='both')
    ax.set_ylabel(r'$\left| \frac{V_{out}}{V_{in}} \right|$')
    ax.set_ylim(0, 2.5)
axs[2].set_xlabel('Frequency, Hz')
fig.tight_layout()
plt.show()

Perhaps others may find it useful.