FilterBankDesigner now handles sampling frequencies other than 48 kHz; added norm compliance test script. If part of the (1/3 or 1/1) octave band lies above the Nyquist frequency, a high pass filter is returned. If the whole band lies above the Nyquist frequency, a no pass filter is returned.

src/lasp/filter/filterbank_design.py

@@ -186,6 +186,14 @@ class FilterBankDesigner:
         x = self.sanitize_input(x)
         fu_n = fu / fnyq
 
+        # Handle fl & fu >= fNyquist: return 'no pass'
+        if fl_n >= 1:
+            return np.tile(np.asarray([0, 0, 0, 1, 0, 0]), (SOS_ORDER, 1))
+
+        # Handle fu >= fNyquist: return high pass
+        if fu_n >= 1:
+            return butter(SOS_ORDER, fl_n, output='sos', btype='highpass')
+
         return butter(SOS_ORDER, [fl_n, fu_n], output='sos', btype='band')
 
     def firFreqResponse(self, x, freq):
@@ -450,8 +458,9 @@ class OctaveBankDesigner(FilterBankDesigner):
         if int(self.fs) in [44100, 48000, 96000]:
             return 1.0
         else:
-            raise ValueError('Unimplemented sampling frequency \'{} Hz\' for '
-                             'SOS filter design. Try 48 kHz.'.format(self.fs))
+            warnings.warn(f'Sampling frequency {self.fs} Hz might result in '
+                          'incorrect filters. Use 48 kHz instead.')
+            return 1.0
 
     def sosFac_u(self, x):
         """Right side percentage of change in cut-on frequency for designing
@@ -463,8 +472,9 @@ class OctaveBankDesigner(FilterBankDesigner):
         if int(self.fs) in [44100, 48000, 96000]:
             return 1.0
         else:
-            raise ValueError('Unimplemented sampling frequency \'{} Hz\' for '
-                             'SOS filter design. Try 48 kHz.'.format(self.fs))
+            warnings.warn(f'Sampling frequency {self.fs} Hz might result in '
+                          'incorrect filters. Use 48 kHz instead.')
+            return 1.0
 
 
 class ThirdOctaveBankDesigner(FilterBankDesigner):
@@ -475,8 +485,7 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
         one-third octave bands
 
         Args:
-            fs: Sampling frequency in [Hz]
-
+            fs: Sampling frequency in [Hz] - can be None
         """
         super().__init__(fs)
         self.xs = list(range(-16, 14))
@@ -485,8 +494,7 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
             '25', '31.5', '40', '50', '63', '80', '100', '125', '160', '200',
             '250', '315', '400', '500', '630', '800', '1k', '1.25k', '1.6k',
             '2k', '2.5k', '3.15k', '4k', '5k', '6.3k', '8k', '10k', '12.5k',
-            '16k', '20k'
-        ]
+            '16k', '20k']
 
         assert len(self.xs) == len(self._nominal_txt)
 
@@ -603,25 +611,19 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
         """Left side percentage of change in cut-on frequency for designing the
         filter."""
         # Idea: correct for frequency warping:
-        if np.isclose(self.fs, 48000):
-            return 1.00
-        elif np.isclose(self.fs, 44100):
-            warnings.warn(f'Sampling frequency {self.fs} Hz might result in '
-                          'incorrect filters')
-            return 1.00
-        elif np.isclose(self.fs, 32768):
-            return 1.00
+        if int(self.fs) in [48000]:
+            return 1.0
         else:
-            raise ValueError('Unimplemented sampling frequency \'{} Hz\' for '
-                             'SOS filter design. Try 48 kHz.'.format(self.fs))
+            warnings.warn(f'Sampling frequency {self.fs} Hz might result in '
+                          'incorrect filters. Use 48 kHz instead.')
+            return 1.0
 
     def sosFac_u(self, x):
         """Right side percentage of change in cut-on frequency for designing
         the filter."""
-        if np.isclose(self.fs, 48000):
-            return 1
-        elif np.isclose(self.fs, 32768):
-            return 1.00
+        if int(self.fs) in [48000]:
+            return 1.0
         else:
-            raise ValueError('Unimplemented sampling frequency \'{} Hz\' for '
-                             'SOS filter design. Try 48 kHz.'.format(self.fs))
+            warnings.warn(f'Sampling frequency {self.fs} Hz might result in '
+                          'incorrect filters. Use 48 kHz instead.')
+            return 1.0

test/scipy_sos_thirdoctavebank.py

@@ -0,0 +1,63 @@
+"""
+This script checks whether the 3rd octave equalizer is norm compliant.
+
+Created on Tue Dec 31 14:43:19 2019
+@author: anne
+"""
+
+from lasp.filter import (ThirdOctaveBankDesigner)
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal import freqz, iirdesign, sosfreqz, butter, cheby2
+
+
+plt.close('all')
+
+# %% Settings of 'filter under test'
+fs = 44.1e3
+nsections = 5  # filter property
+
+designer = ThirdOctaveBankDesigner(fs)
+
+xmin = designer.nominal_txt_tox('25')  # lowest highest band, nominal name
+xmax = designer.nominal_txt_tox('20k')  # highest band
+
+# %% Run check
+compliant = True  # remains True if all filters comply
+fnyq = fs/2
+fll = designer.fl(xmin) / 1.1  # plot limits
+fuu = designer.fu(xmax) * 1.1
+fig = plt.figure()
+plt.xlabel('f (Hz)')
+plt.ylabel('Magnitude (dB)')
+for x in range(xmin, xmax+1):  # loop over bands
+
+    fl = designer.fl(x)  # lower cut off
+    fu = designer.fu(x)  # upper
+
+    fln = fl/fnyq  # lower cutoff, normalized [0, 1]
+    fun = fu/fnyq  # upper
+
+    print(designer.nominal_txt(x))
+
+    sos = designer.createSOSFilter(x)
+    fnorm, h = sosfreqz(sos, worN=2**18, whole=False)
+    freq = fnorm*fnyq/np.pi
+
+    h_dB = 20*np.log10(np.abs(h) + np.finfo(float).eps)
+    if not designer.testStandardCompliance(x, freq, h_dB):
+        print('==============================')
+        print(designer.nominal_txt(x), ' does not comply')
+        compliant = False
+
+    plt.semilogx(freq, h_dB)
+
+    # Upper / lower limits
+    freqlim, ulim, llim = designer.band_limits(x)
+    plt.semilogx(freqlim, ulim)
+    plt.semilogx(freqlim, llim)
+
+    plt.ylim(-5, 1)
+    plt.xlim(fll, fuu)
+
+print(f"\nPassed test : {compliant}")