Implementation of The Audio EQ Cookbook.


In [227]:
%pylab inline
from scipy import signal
# parameters
fs=44100.


Populating the interactive namespace from numpy and matplotlib

In [231]:
# design the filter
# ftype: LPF, HPF, BPF, BPF2, notch, APF,
#        peakingEQ, lowShelf, highShelf
def design_filter(ftype, f0, **kwargs):
    # optional params with defaults
    Q = kwargs.get('Q', .5)
    BW = kwargs.get('BW', .1)
    S = kwargs.get('S', 1.) # shelf only
    dBgain = kwargs.get('dBgain', 20.) # peakingEQ + shelf only
    
    A = 10**(dBgain/40.) # peakingEQ + shelf only
    w0 = 2.*pi*f0/fs
    cw0 = cos(w0)
    sw0 = sin(w0)
    if Q:
        alpha = sw0/(2*Q)
    elif BW:
        alpha = sw0*sinh((log(2)/2)*BW*w0/sw0)
    elif S:
        alpha = (sw0/2) * sqrt( (A+1/A)*(1/S - 1) + 2  )
    else:
        print("no quality, bandwidth or slope specified.")

    # filter coeffs
    if ftype=='LPF':
        b=[(1-cw0)/2, 1-cw0, (1-cw0)/2]
        a=[1+alpha, -2*cw0, 1-alpha]
    elif ftype=='HPF':
        b=[(1+cw0)/2,-(1+cw0),(1+cw0)/2]
        a=[1+alpha,-2*cw0,1-alpha]
    elif ftype=='BPF':
        b=[alpha,0,-alpha]
        a=[1+alpha,-2*cw0,1-alpha]
    elif ftype=='BPF2':
        b=[alpha,0,-alpha]
        a=[1+alpha,-2*cw0, 1-alpha]
    elif ftype=='notch':
        b=[1,-2*cw0,1]
        a=[1+alpha, -2*cw0, 1-alpha]
    elif ftype=='APF':
        b=[1-alpha,-2*cw0,1+alpha]
        a=[1+alpha, -2*cw0, 1-alpha]
    elif ftype=='peakingEQ':
        b=[1+alpha*A, -2*cw0, 1-alpha*A]
        a=[1+alpha/A, -2*cw0, 1-alpha/A]
    elif ftype=='lowShelf':
        b=[A*(  A+1 - (A-1)*cw0+2*sqrt(A)*alpha  ),
           2*A*(  A-1 - (A+1)*cw0  ),
           A*(  A+1 - (A-1)*cw0 - 2*sqrt(A)*alpha)]
        a=[A+1+(A-1)*cw0+2*sqrt(A)*alpha,
           -2*(  A-1 + (A+1)*cw0  ),
           A+1 + (A-1)*cw0 - 2*sqrt(A)*alpha]
    elif ftype=='highShelf':
        b=[A*(  A+1 + (A-1)*cw0+2*sqrt(A)*alpha  ),
           -2*A*(  A-1 + (A+1)*cw0  ),
           A*(  A+1 + (A-1)*cw0 - 2*sqrt(A)*alpha)]
        a=[A+1-(A-1)*cw0+2*sqrt(A)*alpha,
           2*(  A-1 - (A+1)*cw0  ),
           A+1 - (A-1)*cw0 - 2*sqrt(A)*alpha]

    return a, b

In [232]:
def get_impulse_response1(a,b):
    w, h = signal.freqz(b, a, 2048)
    f = (fs/2./pi)*w
    return f, h

def get_impulse_response2(a,b):
    l = 4096
    # dirac pulse
    x = [0., 0., 1.] + [0.] * (l - 1)
    # normalize coefficients
    b0 = b[0]/a[0]
    b1 = b[1]/a[0]
    b2 = b[2]/a[0]
    a1 = a[1]/a[0]
    a2 = a[2]/a[0]
    # filter implementation, direct form 1
    y = []
    y1 = 0.
    y2 = 0.
    for n in xrange(2, l + 2):
        y0 =  b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y1 - a2*y2
        y.append(y0)
        y2 = y1
        y1 = y0
        
    h = fft.rfft(y)
    f = fft.rfftfreq(len(y), d=1./fs)
    return f, h

In [233]:
def plot_impulse_response(a,b,label):
    f, h = get_impulse_response2(a,b)
    #plot(f, 20 * np.log10(abs(h)), label=label)
    semilogx(f, abs(h), label=label)
    xlabel('Frequency (Hz)')
    xlim([0,fs/2.])
    ylabel('Amplitude (db)')
    autoscale(enable=True, axis='y')
    grid('on')
    legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

Testing the various filters. The code below plots impulse responses for a few variants.


In [234]:
# low pass: f0
ftype = 'LPF'
a1,b1 = design_filter(ftype, 100.0)
a2,b2 = design_filter(ftype, 1000.0)
a3,b3 = design_filter(ftype, 10000.0)
plot_impulse_response(a1, b1, "f0=100.0")
plot_impulse_response(a2, b2, "f0=1000.0")
plot_impulse_response(a3, b3, "f0=10000.0")
show()



In [235]:
# low pass: Q
ftype = 'LPF'
a1,b1 = design_filter(ftype, 1000.0, Q=.1)
a2,b2 = design_filter(ftype, 1000.0, Q=1.)
a3,b3 = design_filter(ftype, 1000.0, Q=2.)
plot_impulse_response(a1, b1, "f0=1000.0, q=.1")
plot_impulse_response(a2, b2, "f0=1000.0, q=1.")
plot_impulse_response(a3, b3, "f0=1000.0, q=2.")
show()



In [236]:
# high pass
ftype = 'HPF'
a1,b1 = design_filter(ftype, 100.0)
a2,b2 = design_filter(ftype, 1000.0)
a3,b3 = design_filter(ftype, 10000.0)
plot_impulse_response(a1, b1, "f0=100.0")
plot_impulse_response(a2, b2, "f0=1000.0")
plot_impulse_response(a3, b3, "f0=10000.0")
show()



In [237]:
# band pass
ftype = 'BPF'
a1,b1 = design_filter(ftype, 100.0)
a2,b2 = design_filter(ftype, 1000.0)
a3,b3 = design_filter(ftype, 10000.0)
plot_impulse_response(a1, b1, "f0=100.0")
plot_impulse_response(a2, b2, "f0=1000.0")
plot_impulse_response(a3, b3, "f0=10000.0")
show()



In [220]:
# peaking eq: f0
ftype = 'peakingEQ'
a1,b1 = design_filter(ftype, 100.0)
a2,b2 = design_filter(ftype, 1000.0)
a3,b3 = design_filter(ftype, 10000.0)
plot_impulse_response(a1, b1, "f0=100.0")
plot_impulse_response(a2, b2, "f0=1000.0")
plot_impulse_response(a3, b3, "f0=10000.0")
show()



In [238]:
# peaking eq: Q
ftype = 'peakingEQ'
a1,b1 = design_filter(ftype, 100.0, Q=.1)
a2,b2 = design_filter(ftype, 100.0, Q=.5)
a3,b3 = design_filter(ftype, 100.0, Q=1.)
plot_impulse_response(a1, b1, "f0=100.0, q=.1")
plot_impulse_response(a2, b2, "f0=100.0, q=.5")
plot_impulse_response(a3, b3, "f0=100.0, q=1.")
show()



In [239]:
# peaking eq: dbgain
ftype = 'peakingEQ'
a1,b1 = design_filter(ftype, 100.0, Q=.5, dBgain=-10)
a2,b2 = design_filter(ftype, 100.0, Q=.5, dBgain=0)
a3,b3 = design_filter(ftype, 100.0, Q=.5, dBgain=+10)
plot_impulse_response(a1, b1, "f0=100.0, dbgain=-10")
plot_impulse_response(a2, b2, "f0=100.0, dbgain=0")
plot_impulse_response(a3, b3, "f0=100.0, dbgain=+10")
show()



In [ ]: