2
$\begingroup$

I am building a second order modified biquad filter based on an implementation by Will Pirkle in "designing audio effect plugins in c++". However when I look at my implementation's frequency response, things look off: e.g., a +6 dB gain results in more boost than +12 dB. I am basically implementing what I see in the block diagram and design equations from the book: enter image description here enter image description here

My filter looks like this:

import numpy as np
FS = 44100
def first_order_low_shelf(sig, fc, gain, coeffs_only=False):
    y = [0,0]
    x = sig
    # design equations
    theta_c = (2.*np.pi*fc)/FS
    mu = 10.**(gain/20.)
    beta = 4./(1.+mu)
    delta = beta*np.tan(theta_c/2.)
    gamma = (1-delta)/(1+delta)
    a0 = (1-gamma)/2
    a1 = a0
    a2 = 0.
    b1 = -1*gamma
    b2 = 0.
    c0 = mu - 1.
    d0 = 1.
    if coeffs_only:
        numerator   = [c0*a0+d0, c0*a1, c0*a2]
        denominator = [1, c0*b1, c0*b2]
        return numerator, denominator
    
    for i in range(2, len(sig)):
        #difference equation 
        yn = c0*(a0*x[i] + a1*x[i-1] + a2*x[i-2] - b1*y[i-1] - b2*y[i-2]) + d0*x[i]
        y.append(yn)
        
    return y

I am looking at the frequency response for a fixed corner frequency, over gains between [-12,12] using scipy.signal.freqz, which uses the transfer function coefficients as input:

from scipy.signal import freqz
fc = 400
for db in [-12,-6,-3,3,6,12]:
    nums, denoms = first_order_low_shelf(None, fc, db, True)
    f, h = freqz(nums, denoms)
    gain = 20 * np.log10(abs(h))
    plt.semilogx(f*FS/(2*np.pi), gain, label=db)
plt.legend()
plt.xlabel('freqency (Hz)')
plt.ylabel('gain (db)')

enter image description here

The numerator and denominator arrays in first_order_low_shelf reflect those of the transfer function but I can't see why I'm getting such bizarre frequency responses.

$\endgroup$
5
  • 2
    $\begingroup$ If I were you, I would get rid of the redundancies. Both $c_0$ and $d_0$ are redundant coefficients. You can set $d_0=0$ and $c_0=1$ and make adjustments to $a_0, a_1, a_2$ and get precisely the same frequency response. What basis to set the remaining 5 coefficients? I would suggest the Audio EQ Cookbook. I dunno what shape of shelving filter that Will's method shoots for and there are different variations. But you don't need 7 coefficients for a 2nd-order IIR filter. There are only 5 degrees of freedom and you can do it all with 5. $\endgroup$ Nov 28, 2023 at 1:03
  • 2
    $\begingroup$ And, for consistency in nomenclature, I would swap the $a_n$ labels with the $b_n$ labels and be consistent with Oppenhiem and Schafer and most other text books. "$a_1, a_2$" go into the denominator. $\endgroup$ Nov 28, 2023 at 1:05
  • 2
    $\begingroup$ woah! You got $b_2=0$ and $a_2=0$? That means your second-order filter is really a first-order filter. This is depicted way too complicated for a 1st-order shelving filter. There are only 3 degrees of freedom. You should have only 3 coefficients. $\endgroup$ Nov 28, 2023 at 3:25
  • $\begingroup$ Looking at these design equations through a more critical lens is helpful. I now see these redundancies; unfortunately, the chapter glosses over these pedagogical choices - I wish I knew why Pirkle decided to present them in such a fashion. $\endgroup$ Nov 28, 2023 at 21:53
  • $\begingroup$ Long ago, sometime in the 1990's, at an AES convention in LA, I had a one-on-one lunch or dinner with Will. It was when he first started at UMiami and Ken Pohlmann was the head of the Music Engineering department. It's really hard to decode someone else's pedagogy that doesn't entirely make sense at first glance. $\endgroup$ Nov 29, 2023 at 2:43

1 Answer 1

3
$\begingroup$

Your filter coefficients and your difference equation are wrong. You assumed that the output signal $y[n]$ is fed back to the input, but what is actually fed back is the signal before the multiplier $c_0$.

The correct coefficients should be

num = [c0*a0+d0, c0*a1+d0*b1, c0*a2+d0*b2];
den = [1, b1, b2];

With these changes, the frequency responses look like they should:

enter image description here

Derivation of the transfer function:

Let $W(z)$ be $\mathcal{Z}$-transform of the signal before the multiplier $c_0$, and let $H(z)=N(z)/D(z)$ be the transfer function of the second-order filter with coefficients $a_i$ and $b_i$. Then

$$W(z)=H(z)X(z)=\frac{N(z)}{D(z)}X(z)$$

and

$$ Y(z) = c_0W(z) + d_0X(z) = c_0\frac{N(z)}{D(z)}X(z) + d_0X(z) $$

from which it follows that

$$\frac{Y(z)}{X(z)}=c_0\frac{N(z)}{D(z)}+d_0=\frac{c_0N(z)+d_0D(z)}{D(z)}$$

$\endgroup$
2
  • $\begingroup$ working backwards from your transfer function, the difference equation would be y(n) = (c0*a0+d0)*x(n) + (c0*a1+d0*b1)*x(n-1) + (c0*a2+d0*b2)*x(n-2) - b1*y(n-1) -b2*y(n-2). This checks out numerically. However, I'm wondering how you were able to derive this from the block diagram? E.g., it looks like d0 just scales x(n), and I don't see where d0 scales the b coefficients, or how b scales the inputs with positive sign. $\endgroup$ Nov 28, 2023 at 22:06
  • 1
    $\begingroup$ @BrianBarry: I've added the derivation of the transfer function to my answer. $\endgroup$
    – Matt L.
    Nov 29, 2023 at 7:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.