Skip to content

Instantly share code, notes, and snippets.

@mryssng

mryssng/biquad_filter.py

Last active February 8, 2020 15:33
Show Gist options
  • Save mryssng/8a33f908864d144f4ed5732d65808db6 to your computer and use it in GitHub Desktop.
Save mryssng/8a33f908864d144f4ed5732d65808db6 to your computer and use it in GitHub Desktop.
Download ZIP
Implementation of EQ Biquad Filter for python
Raw
biquad_filter.py
# -*- coding: utf-8 -*-
"""
Created on 2020/01/31
Description:
Show digital filters frequency responce.
Based on "Cookbook formulae for audio EQ biquad filter coefficients"
by Robert Bristow-Johnson.
Author: Ytse Jam
Copyright © 2020 Ytse Jam. Licensed under MIT.
"""
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
def allPass(fc, q):
a = [0.0] * 3
b = [0.0] * 3
w0 = 2.0 * np.pi * fc / 44100
sinw0 = np.sin(w0)
cosw0 = np.cos(w0)
alpha = sinw0 / (2.0 * q)
a[0] = 1.0 + alpha
a[1] = -2.0 * cosw0
a[2] = 1.0 - alpha
b[0] = 1.0 - alpha
b[1] = -2.0 * cosw0
b[2] = 1.0 + alpha
return a, b
def peakingEQ(fc, dBGain, bw):
a = [0.0] * 3
b = [0.0] * 3
w0 = 2.0 * np.pi * fc / 44100
sinw0 = np.sin(w0)
cosw0 = np.cos(w0)
alpha = sinw0 * np.sinh(np.log(2.0) / 2.0 * bw * w0 / sinw0)
if dBGain != 0.0:
A = np.power(10.0, (dBGain / 40.0))
else:
A = 1.0
a[0] = 1.0 + alpha / A
a[1] = -2.0 * cosw0
a[2] = 1.0 - alpha / A
b[0] = 1.0 + alpha * A
b[1] = -2.0 * cosw0
b[2] = 1.0 - alpha * A
return a, b
def lowShelf(fc, dBGain, bw):
a = [0.0] * 3
b = [0.0] * 3
w0 = 2.0 * np.pi * fc / 44100
sinw0 = np.sin(w0)
cosw0 = np.cos(w0)
alpha = sinw0 * np.sinh(np.log(2.0) / 2.0 * bw * w0 / sinw0)
if dBGain != 0.0:
A = np.power(10.0, (dBGain / 40.0))
else:
A = 1.0
a[0] = (A + 1.0) + (A - 1.0) * cosw0 + 2.0 * np.sqrt(A) * alpha
a[1] = -2.0 * ((A - 1.0) + (A + 1.0) * cosw0)
a[2] = (A + 1.0) + (A - 1.0) * cosw0 - 2.0 * np.sqrt(A) * alpha
b[0] = A * ((A + 1.0) - (A - 1.0) * cosw0 + 2.0 * np.sqrt(A) * alpha)
b[1] = 2.0 * A * ((A - 1.0) - (A + 1.0) * cosw0)
b[2] = A * ((A + 1.0) - (A - 1.0) * cosw0 - 2.0 * np.sqrt(A) * alpha)
return a, b
def highShelf(fc, dBGain, bw):
a = [0.0] * 3
b = [0.0] * 3
w0 = 2.0 * np.pi * fc / 44100
sinw0 = np.sin(w0)
cosw0 = np.cos(w0)
alpha = sinw0 * np.sinh(np.log(2.0) / 2.0 * bw * w0 / sinw0)
if dBGain != 0.0:
A = np.power(10.0, (dBGain / 40.0))
else:
A = 1.0
a[0] = (A + 1.0) - (A - 1.0) * cosw0 + 2.0 * np.sqrt(A) * alpha
a[1] = 2.0 * ((A - 1.0) - (A + 1.0) * cosw0)
a[2] = (A + 1.0) - (A - 1.0) * cosw0 - 2.0 * np.sqrt(A) * alpha
b[0] = A * ((A + 1.0) + (A - 1.0) * cosw0 + 2.0 * np.sqrt(A) * alpha)
b[1] = -2.0 * A * ((A - 1.0) + (A + 1.0) * cosw0)
b[2] = A * ((A + 1.0) + (A - 1.0) * cosw0 - 2.0 * np.sqrt(A) * alpha)
return a, b
def plotData(w, h1, h2=0, h3=0, h4=0, h5=0, h6=0, title="", name=""):
plt.clf()
f = w / (2.0 * np.pi) * 44100
plt.subplot(211)
if h1 != 0:
a1 = 20 * np.log10(abs(h1[1]))
plt.plot(f, a1, label=h1[0])
else:
a1 = 0
if h2 != 0:
a2 = 20 * np.log10(abs(h2[1]))
plt.plot(f, a2, label=h2[0])
else:
a2 = 0
if h3 != 0:
a3 = 20 * np.log10(abs(h3[1]))
plt.plot(f, a3, label=h3[0])
else:
a3 = 0
if h4 != 0:
a4 = 20 * np.log10(abs(h4[1]))
plt.plot(f, a4, label=h4[0])
else:
a4 = 0
if h5 != 0:
a5 = 20 * np.log10(abs(h5[1]))
plt.plot(f, a5, label=h5[0])
else:
a5 = 0
if h6 != 0:
a6 = 20 * np.log10(abs(h6[1]))
plt.plot(f, a6, label=h6[0])
else:
a6 = 0
# ymax = max(max(a1), max(a2), max(a3), max(a4), max(a5), max(a6))
# ymin = min(min(a1), min(a2), min(a3), min(a4), min(a5), min(a6))
# ymax = max(ymax, abs(ymin))
plt.legend()
plt.xscale('log')
plt.xticks([10, 100, 1000, 10000, 100000], ["10", "100", "1k", "10k", "100k"])
if title == '':
plt.title('Magnitude')
else:
plt.title(title + ' ' + 'Magnitude')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.xlim(10, 100000)
# plt.ylim(-5, 5)
plt.grid(which='both', axis='both')
plt.subplot(212)
if h1 != 0:
plt.plot(f, np.angle(h1[1]) * 180 / np.pi, label=h1[0])
if h2 != 0:
plt.plot(f, np.angle(h2[1]) * 180 / np.pi, label=h2[0])
if h3 != 0:
plt.plot(f, np.angle(h3[1]) * 180 / np.pi, label=h3[0])
if h4 != 0:
plt.plot(f, np.angle(h4[1]) * 180 / np.pi, label=h4[0])
if h5 != 0:
plt.plot(f, np.angle(h5[1]) * 180 / np.pi, label=h5[0])
if h6 != 0:
plt.plot(f, np.angle(h6[1]) * 180 / np.pi, label=h6[0])
plt.legend()
plt.xscale('log')
plt.xticks([10, 100, 1000, 10000, 100000], ["10", "100", "1k", "10k", "100k"])
plt.yticks([-180, -90, 0, 90, 180], ["-180", "-90", "0", "90", "180"])
if title == '':
plt.title('Phase')
else:
plt.title(title + ' ' + 'Phase')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [deg]')
plt.margins(0, 0.1)
plt.xlim(10, 100000)
plt.ylim(-180, 180)
plt.grid(which='both', axis='both')
plt.tight_layout()
if name == "":
plt.show()
else:
plt.savefig(name)
def main():
# PeakingEQ
a1, b1 = peakingEQ(1000, 5, 0.05) # cutoff = 150Hz, gain = +5, bw = 0.05
a2, b2 = peakingEQ(1000, 5, 0.50) # cutoff = 150Hz, gain = +5, bw = 0.50
a3, b3 = peakingEQ(1000, 5, 1.00) # cutoff = 150Hz, gain = +5, bw = 1.00
a4, b4 = peakingEQ(1000, -5, 0.05) # cutoff = 150Hz, gain = -5, bw = 0.05
a5, b5 = peakingEQ(1000, -5, 0.50) # cutoff = 150Hz, gain = -5, bw = 0.50
a6, b6 = peakingEQ(1000, -5, 1.00) # cutoff = 150Hz, gain = -5, bw = 1.00
w, h1 = signal.freqz(b1, a1, worN=1024*5)
w, h2 = signal.freqz(b2, a2, worN=1024*5)
w, h3 = signal.freqz(b3, a3, worN=1024*5)
w, h4 = signal.freqz(b4, a4, worN=1024*5)
w, h5 = signal.freqz(b5, a5, worN=1024*5)
w, h6 = signal.freqz(b6, a6, worN=1024*5)
plotData(w, \
h1=['g=+5, bw=0.05', h1], \
h2=['g=+5, bw=0.10', h2], \
h3=['g=+5, bw=1.00', h3], \
h4=['g=-5, bw=0.05', h4], \
h5=['g=-5, bw=0.10', h5], \
h6=['g=-5, bw=1.00', h6], \
title='PeakingEQ', \
name='peakingEQ_g+5.png')
# AllPass
a1, b1 = allPass(1000, 1.0) # cutoff = 1000Hz, bw = 1.0
w, h = signal.freqz(b1, a1, worN=1024*5)
plotData(w, h1=['cutoff=1000, bw=1.0', h1], title='AllPass', name='allPass.png')
# LowShelf
a1, b1 = lowShelf(300, 2.0, 0.05) # cutoff = 300Hz, gain = +2, bw = 0.05
a2, b2 = lowShelf(300, 2.0, 0.50) # cutoff = 300Hz, gain = +2, bw = 0.50
a3, b3 = lowShelf(300, 2.0, 1.0) # cutoff = 300Hz, gain = +2, bw = 1.0
a4, b4 = lowShelf(300, 2.0, 2.0) # cutoff = 300Hz, gain = +2, bw = 2.0
a5, b5 = lowShelf(300, 2.0, 5.0) # cutoff = 300Hz, gain = +2, bw = 5.0
a6, b6 = lowShelf(300, 2.0, 10.0) # cutoff = 300Hz, gain = +2, bw = 10.0
w, h1 = signal.freqz(b1, a1, worN=1024*5)
w, h2 = signal.freqz(b2, a2, worN=1024*5)
w, h3 = signal.freqz(b3, a3, worN=1024*5)
w, h4 = signal.freqz(b4, a4, worN=1024*5)
w, h5 = signal.freqz(b5, a5, worN=1024*5)
w, h6 = signal.freqz(b6, a6, worN=1024*5)
plotData(w, \
h1=['g=+2, bw=0.05', h1], \
h2=['g=+2, bw=0.50', h2], \
h3=['g=+2, bw=1.0', h3], \
h4=['g=+2, bw=2.0', h4], \
h5=['g=+2, bw=5.0', h5], \
h6=['g=+2, bw=10.0', h6], \
title='LowShelf', \
name='lowShelf_g+2.png')
# HighShelf
a1, b1 = highShelf(1000, 5.0, 0.05) # cutoff = 1000Hz, gain = +5, bw = 0.05
a2, b2 = highShelf(1000, 5.0, 0.50) # cutoff = 1000Hz, gain = +5, bw = 0.50
a3, b3 = highShelf(1000, 5.0, 1.0) # cutoff = 1000Hz, gain = +5, bw = 1.0
a4, b4 = highShelf(1000, 5.0, 2.0) # cutoff = 1000Hz, gain = +5, bw = 2.0
a5, b5 = highShelf(1000, 5.0, 5.0) # cutoff = 1000Hz, gain = +5, bw = 5.0
a6, b6 = highShelf(1000, 5.0, 10.0) # cutoff = 1000Hz, gain = +5, bw = 10.0
w, h1 = signal.freqz(b1, a1, worN=1024*5)
w, h2 = signal.freqz(b2, a2, worN=1024*5)
w, h3 = signal.freqz(b3, a3, worN=1024*5)
w, h4 = signal.freqz(b4, a4, worN=1024*5)
w, h5 = signal.freqz(b5, a5, worN=1024*5)
w, h6 = signal.freqz(b6, a6, worN=1024*5)
plotData(w, \
h1=['g=+5, bw=0.05', h1], \
h2=['g=+5, bw=0.50', h2], \
h3=['g=+5, bw=1.0', h3], \
h4=['g=+5, bw=2.0', h4], \
h5=['g=+5, bw=5.0', h5], \
h6=['g=+5, bw=10.0', h6], \
title='HighShelf', \
name='highShelf_g+5.png')
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment