Removed biquad and reworked peq2geq.py to use PEQ.

frequency_response.py

@@ -17,15 +17,14 @@
 from PIL import Image
 import re
 import warnings
-import biquad
 from constants import DEFAULT_F_MIN, DEFAULT_F_MAX, DEFAULT_STEP, DEFAULT_MAX_GAIN, DEFAULT_TREBLE_F_LOWER, \
     DEFAULT_TREBLE_F_UPPER, DEFAULT_TREBLE_GAIN_K, DEFAULT_SMOOTHING_WINDOW_SIZE, \
     DEFAULT_SMOOTHING_ITERATIONS, DEFAULT_TREBLE_SMOOTHING_F_LOWER, DEFAULT_TREBLE_SMOOTHING_F_UPPER, \
     DEFAULT_TREBLE_SMOOTHING_WINDOW_SIZE, DEFAULT_TREBLE_SMOOTHING_ITERATIONS, DEFAULT_TILT, DEFAULT_FS, \
     DEFAULT_F_RES, DEFAULT_BASS_BOOST_GAIN, DEFAULT_BASS_BOOST_FC, \
     DEFAULT_BASS_BOOST_Q, DEFAULT_GRAPHIC_EQ_STEP, HARMAN_INEAR_PREFENCE_FREQUENCIES, \
     HARMAN_ONEAR_PREFERENCE_FREQUENCIES, PREAMP_HEADROOM, DEFAULT_MAX_SLOPE, PEQ_CONFIGS
-from peq import PEQ, OptimizationFinished
+from peq import PEQ, Peaking
 
 warnings.filterwarnings("ignore", message="Values in x were outside bounds during a minimize step, clipping to bounds")
 
@@ -598,16 +597,12 @@ def create_target(self,
         Returns:
             Target for equalization
         """
-        bass_boost = biquad.digital_coeffs(
-            self.frequency,
-            DEFAULT_FS,
-            *biquad.low_shelf(bass_boost_fc, bass_boost_q, bass_boost_gain, DEFAULT_FS)
-        )
+        bass_boost = Peaking(self.frequency, DEFAULT_FS, fc=bass_boost_fc, q=bass_boost_q, gain=bass_boost_gain)
         if tilt is not None:
             tilt = self._tilt(tilt=tilt)
         else:
             tilt = np.zeros(len(self.frequency))
-        return bass_boost + tilt
+        return bass_boost.fr + tilt
 
     def compensate(self,
                    compensation,

peq.py

@@ -123,11 +123,10 @@ def sharpness_penalty(self):
 
 
 class Peaking(PEQFilter):
-    def __init__(self, f, fs,
-                 fc=None, optimize_fc=None, min_fc=DEFAULT_PEAKING_FILTER_MIN_FC, max_fc=DEFAULT_PEAKING_FILTER_MAX_FC,
-                 q=None, optimize_q=None, min_q=DEFAULT_PEAKING_FILTER_MIN_Q, max_q=DEFAULT_PEAKING_FILTER_MAX_Q,
-                 gain=None, optimize_gain=None, min_gain=DEFAULT_PEAKING_FILTER_MIN_GAIN,
-                 max_gain=DEFAULT_PEAKING_FILTER_MAX_GAIN):
+    def __init__(self, f, fs, fc=None, optimize_fc=None, min_fc=DEFAULT_PEAKING_FILTER_MIN_FC,
+                 max_fc=DEFAULT_PEAKING_FILTER_MAX_FC, q=None, optimize_q=None, min_q=DEFAULT_PEAKING_FILTER_MIN_Q,
+                 max_q=DEFAULT_PEAKING_FILTER_MAX_Q, gain=None, optimize_gain=None,
+                 min_gain=DEFAULT_PEAKING_FILTER_MIN_GAIN, max_gain=DEFAULT_PEAKING_FILTER_MAX_GAIN):
         super().__init__(f, fs, fc, optimize_fc, min_fc, max_fc, q, optimize_q, min_q, max_q, gain, optimize_gain,
                          min_gain, max_gain)
 
@@ -221,11 +220,10 @@ def sharpness_penalty(self):
 
 
 class ShelfFilter(PEQFilter, ABC):
-    def __init__(self, f, fs,
-                 fc=None, optimize_fc=None, min_fc=DEFAULT_SHELF_FILTER_MIN_FC, max_fc=DEFAULT_SHELF_FILTER_MAX_FC,
-                 q=None, optimize_q=None, min_q=DEFAULT_SHELF_FILTER_MIN_Q, max_q=DEFAULT_SHELF_FILTER_MAX_Q,
-                 gain=None, optimize_gain=None, min_gain=DEFAULT_SHELF_FILTER_MIN_GAIN,
-                 max_gain=DEFAULT_SHELF_FILTER_MAX_GAIN):
+    def __init__(self, f, fs, fc=None, optimize_fc=None, min_fc=DEFAULT_SHELF_FILTER_MIN_FC,
+                 max_fc=DEFAULT_SHELF_FILTER_MAX_FC, q=None, optimize_q=None, min_q=DEFAULT_SHELF_FILTER_MIN_Q,
+                 max_q=DEFAULT_SHELF_FILTER_MAX_Q, gain=None, optimize_gain=None,
+                 min_gain=DEFAULT_SHELF_FILTER_MIN_GAIN, max_gain=DEFAULT_SHELF_FILTER_MAX_GAIN):
         super().__init__(f, fs, fc, optimize_fc, min_fc, max_fc, q, optimize_q, min_q, max_q, gain, optimize_gain,
                          min_gain, max_gain)
 
@@ -363,16 +361,18 @@ def __init__(self, f, fs, filters=None, target=None,
             for filt in filters:
                 self.add_filter(filt)
         self.target = np.array(target) if target is not None else None
-        self._ix50 = np.sum(self.f < 50)
-        self._ix10k = np.sum(self.f < 10e3)  # Index for ~10 kHz
-        self._ix20k = np.sum(self.f < 20e3)  # Index for ~20 kHz
-        self.history = None
         self._min_f = min_f
         self._max_f = max_f
+        self._min_f_ix = np.argmin(np.abs(self.f - self._min_f))
+        self._max_f_ix = np.argmin(np.abs(self.f - self._max_f))
+        self._ix50 = np.argmin(np.abs(self.f - 50))
+        self._10k_ix = np.argmin(np.abs(self.f - 10000))
+        self._20k_ix = np.argmin(np.abs(self.f - 20000))
         self._max_time = max_time
         self._target_loss = target_loss
         self._min_change_rate = min_change_rate
         self._min_std = min_std
+        self.history = None
 
     @classmethod
     def from_dict(cls, config, fs, target=None):
@@ -414,7 +414,7 @@ def from_dict(cls, config, fs, target=None):
 
     def add_filter(self, filt):
         if filt.fs != self.fs:
-            raise ValueError('Filter sampling rate (fs) must match equalizer sampling rate')
+            raise ValueError(f'Filter sampling rate ({filt.fs}) must match equalizer sampling rate ({self.fs})')
         if not np.array_equal(filt.f, self.f):
             raise ValueError('Filter frequency array (f) must match equalizer frequency array')
         self.filters.append(filt)
@@ -472,17 +472,16 @@ def _optimizer_loss(self, params, parse=True):
         if parse:
             self._parse_optimizer_params(params)
 
-        # Above 10 kHz only the total energy matters so we'll take mean of values between 10 kHz and 20 kHz
+        # Above 10 kHz only the total energy matters so we'll take the average
         fr = self.fr.copy()
         target = self.target.copy()
-        target[self._ix10k:self._ix20k] = np.mean(target[self._ix10k:self._ix20k])
-        fr[self._ix10k:self._ix20k] = np.mean(self.fr[self._ix10k:self._ix20k])
+        target[self._10k_ix:] = np.mean(target[self._10k_ix:])
+        fr[self._10k_ix:] = np.mean(self.fr[self._10k_ix:])
         #target[:self._ix50] = np.mean(target[:self._ix50])  # TODO: Is this good?
         #fr[:self._ix50] = np.mean(fr[:self._ix50])
 
-        # Mean squared error as loss, only up to 20 kHz
-        loss_val = np.mean(np.square(target[:self._ix20k] - fr[:self._ix20k]))
-        # loss_val = np.mean(np.square(target[:self._ix10k] - fr[:self._ix10k]))
+        # Mean squared error as loss, between minimum and maximum frequencies
+        loss_val = np.mean(np.square(target[self._min_f_ix:self._max_f_ix] - fr[self._min_f_ix:self._max_f_ix]))
 
         # Sum penalties from all filters to MSE
         for filt in self.filters:
@@ -617,7 +616,8 @@ def plot(self, fig=None, ax=None):
             ax.grid(True, which='minor')
             ax.xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
             ax.set_xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000])
-        ax.plot(self.f, self.target, color='gray', linewidth=3, label='Target')
+        if self.target is not None:
+            ax.plot(self.f, self.target, color='black', linestyle='--', linewidth=1, label='Target')
         for i, filt in enumerate(self.filters):
             ax.fill_between(filt.f, np.zeros(filt.fr.shape), filt.fr, alpha=0.3, color=f'C{i}')
             ax.plot(filt.f, filt.fr, color=f'C{i}', linewidth=1)

research/peq2geq.py

@@ -1,40 +1,29 @@
 import os
 import sys
+import unittest
 from pathlib import Path
+from argparse import ArgumentParser, SUPPRESS
+import re
+import numpy as np
 ROOT_PATH = Path(__file__).resolve().parent.parent
 if str(ROOT_PATH) not in sys.path:
     sys.path.insert(1, str(ROOT_PATH))
-import numpy as np
-import re
 from frequency_response import FrequencyResponse
-from biquad import peaking, low_shelf, high_shelf, digital_coeffs
-from argparse import ArgumentParser, SUPPRESS
-
-fns = {'PK': peaking, 'LS': low_shelf, 'HS': high_shelf}
-fs = 48000
+from constants import DEFAULT_FS
+from peq import PEQ, Peaking, LowShelf, HighShelf
 
 
-def peq2fr(fc, q, gain, filts):
-    if type(fc) in [float, int]:
-        fc = np.array([fc])
-    if type(q) in [float, int]:
-        q = np.array([q])
-    if type(gain) in [float, int]:
-        gain = np.array([gain])
-    if type(filts) == str:
-        filts = [filts] * len(fc)
-    fr = FrequencyResponse(name='PEG')
-    c = np.zeros(fr.frequency.shape)
-    for i, filt in enumerate(filts):
-        a0, a1, a2, b0, b1, b2 = fns[filt](fc[i], q[i], gain[i], fs=fs)
-        c += digital_coeffs(fr.frequency, fs, a0, a1, a2, b0, b1, b2)
-    fr.raw = c
-    return fr
+classes = {'PK': Peaking, 'LS': LowShelf, 'HS': HighShelf}
 
 
-def peq2geq(fc, q, gain, filts, normalize=False):
-    fr = peq2fr(fc, q, gain, filts)
-    fr.equalization = fr.raw
+def peq2geq(fcs, qs, gains, types, normalize=False):
+    if not (len(fcs) == len(qs) == len(gains) == len(types)):
+        raise ValueError('Different number of Fc, Q, gain and filter types')
+    fr = FrequencyResponse(name='peq2geq')
+    peq = PEQ(fr.frequency, DEFAULT_FS)
+    for fc, q, gain, filter_type in zip(fcs, qs, gains, types):
+        peq.add_filter(classes[filter_type](fr.frequency, DEFAULT_FS, fc=fc, q=q, gain=gain))
+    fr.equalization = peq.fr
     return fr.eqapo_graphic_eq(normalize=normalize)
 
 
@@ -77,22 +66,40 @@ def main():
     parser.add_argument('--normalize', action='store_true', help='Normalize gain?')
     args = parser.parse_args()
     if 'file' in args and args.file:
-        fcs, qs, gains, filts = read_eqapo(args.file)
+        fcs, qs, gains, types = read_eqapo(args.file)
     else:
-        fcs, qs, gains, filts = [], [], [], []
+        fcs, qs, gains, types = [], [], [], []
     if args.fc:
         fcs += args.fc
     if args.q:
         qs += args.q
     if args.gain:
         gains += args.gain
     if args.type:
-        filts += args.type
-    if not (len(fcs) == len(qs) == len(gains) == len(filts)):
-        print('Different number of Fc, Q, gain and filter types')
-        return
-    print(peq2geq(fcs, qs, gains, filts, normalize=args.normalize))
+        types += args.type
+    print(peq2geq(fcs, qs, gains, types, normalize=args.normalize))
 
 
 if __name__ == '__main__':
     main()
+
+
+class TestPeq2Geq(unittest.TestCase):
+    def test(self):
+        f = FrequencyResponse.generate_frequencies()
+        peq = PEQ(f, DEFAULT_FS, filters=[
+            Peaking(f, DEFAULT_FS, fc=500, q=1.41, gain=2),
+            HighShelf(f, DEFAULT_FS, fc=2000, q=0.7, gain=-3)
+        ])
+        geq = peq2geq(
+            [filt.fc for filt in peq.filters],
+            [filt.q for filt in peq.filters],
+            [filt.gain for filt in peq.filters],
+            ['PK', 'HS']
+        )
+        s = geq.split('GraphicEQ: ')[1]
+        pairs = [(float(p.split()[0]), float(p.split()[1])) for p in s.split('; ')]
+        fr_geq = FrequencyResponse(name='geq', frequency=[_f for _f, _g in pairs], raw=[_g for _f, _g in pairs])
+        fr_geq.interpolate()
+        fr_geq.raw -= np.mean(peq.fr - fr_geq.raw)
+        self.assertLess(np.mean(np.abs(peq.fr - fr_geq.raw)), 0.1)