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

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 [ ]: