From 3a894879f0906d4f88203cfe282ab56bb62df3b5 Mon Sep 17 00:00:00 2001
From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" <j.a.dejong@ascee.nl>
Date: Tue, 20 Apr 2021 22:25:43 +0200
Subject: [PATCH] Added biquad filter code

---
 lasp/filter/biquad.py | 43 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 43 insertions(+)
 create mode 100644 lasp/filter/biquad.py

diff --git a/lasp/filter/biquad.py b/lasp/filter/biquad.py
new file mode 100644
index 0000000..41e6566
--- /dev/null
+++ b/lasp/filter/biquad.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""!
+Author: J.A. de Jong - ASCEE V.O.F.
+
+Description: Filter design implementation of common biquad filters that are
+often used in parametric equalizers.
+"""
+__all__ = ['peaking', 'biquadTF']
+from scipy.signal import bilinear_zpk, zpk2sos, freqz_zpk, sosfreqz
+from scipy.interpolate import interp1d
+import numpy as np
+
+def peaking(fs, f0, Q, gain):
+    """
+    Design of peaking biquad filter
+
+    Args:
+        fs: Sampling frequency [Hz]
+        f0: Center frequency
+        Q: Quality factor (~ inverse of bandwidth)
+        gain: Increase in level at the center frequency
+    """
+    A = np.sqrt(10**(gain/20))
+    omg0 = 2*np.pi*f0/fs
+    alpha = np.sin(omg0)/Q/2
+    b0 = 1+alpha*A
+    b1 = -2*np.cos(omg0)
+    b2 = 1-alpha*A
+    a0 = 1 + alpha/A
+    a1 = -2*np.cos(omg0)
+    a2 = 1-alpha/A
+
+    return np.array([b0/a0, b1/a0, b2/a0, a0/a0, a1/a0, a2/a0])
+
+def biquadTF(fs, freq, ba):
+    """
+    Computes the transfer function of the biquad.
+
+    """
+    freq2, h = sosfreqz(ba, worN=48000, fs=fs)
+    interpolator = interp1d(freq2, h, kind='quadratic') 
+    return interpolator(freq)