Biquad Filter in python from http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

a python Bristow-Johnson biquad

#-*- coding: utf8 -*-
# implementation of recipe on: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
import numpy as n, pylab as p

# variaveis de parametrizacao
fa=44100.
f=2.

# tipos: LPF, HPF, BPF, BPF2, notch
# APF, peakingEQ, lowShelf, highShelf
tipo='LPF'

dBgain=20. # para frequência central

Q=.01
BW=.1
S=1.

# variaveis intermediarias
A = 10**(dBgain/40.)
w0=2*n.pi*f/fa
cw0=n.cos(w0)
sw0=n.sin(w0)
if Q:
    alpha = sw0/(2*Q)
elif BW:
    alpha = sw0*n.sinh((n.log(2)/2)*BW*w0/sw0)
elif S:
    alpha = (sw0/2) * n.sqrt( (A+1/A)*(1/S - 1) + 2  )
else:
    print("no quality, bandwidth or slope specified.")

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


# Impulse response by hand
# :::
sinal=[b[0]]
sinal+=[b[1]-sinal[-1]*a[1]]
sinal+=[b[2]-sinal[-1]*a[1]-sinal[-2]*a[2]]
for i in xrange(44100): # for 1 second
    sinal.append(-sinal[-1]*a[1]-sinal[-2]*a[2])

# fft and modulus
fft=n.fft.fft(sinal)
m=n.abs(fft)

# pylab plot and show
p.plot(m)
p.show()



# _o_o_ oOo _o_o_ labMacambira.sf.net _o_o_ oOo _o_o_