low shelving filter frequency response diagnosis

Question

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:

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

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.

Answer 1

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:

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