Bugfix: FilterBankDesigner now always returns second order sections with shape (SOS_ORDER, 6), even for the high pass section. Expanded the third octave bank compliance test. Since it passed for all common sample rates, I removed the warning that the filters might be non-compliant for fs unequal to 48 kHz.

src/lasp/filter/filterbank_design.py

@@ -192,7 +192,10 @@ class FilterBankDesigner:
 
         # Handle fu >= fNyquist: return high pass
         if fu_n >= 1:
-            return butter(SOS_ORDER, fl_n, output='sos', btype='highpass')
+            highpass = butter(SOS_ORDER, fl_n, output='sos', btype='highpass')
+            allpass = np.tile(np.asarray([1, 0, 0, 1, 0, 0]),
+                              (SOS_ORDER-highpass.shape[0], 1))
+            return np.vstack((highpass, allpass))  # same shape=(SOS_ORDER, 6)
 
         return butter(SOS_ORDER, [fl_n, fu_n], output='sos', btype='band')
 
@@ -458,8 +461,6 @@ class OctaveBankDesigner(FilterBankDesigner):
         if int(self.fs) in [44100, 48000, 96000]:
             return 1.0
         else:
-            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):
@@ -472,8 +473,6 @@ class OctaveBankDesigner(FilterBankDesigner):
         if int(self.fs) in [44100, 48000, 96000]:
             return 1.0
         else:
-            warnings.warn(f'Sampling frequency {self.fs} Hz might result in '
-                          'incorrect filters. Use 48 kHz instead.')
             return 1.0
 
 
@@ -614,8 +613,6 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
         if int(self.fs) in [48000]:
             return 1.0
         else:
-            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):
@@ -624,6 +621,4 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
         if int(self.fs) in [48000]:
             return 1.0
         else:
-            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

@@ -13,51 +13,65 @@ 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
+def check_one(fs):
-fnyq = fs/2
+    """Check for fs = fs"""
-fll = designer.fl(xmin) / 1.1  # plot limits
+    # Create filter bank
-fuu = designer.fu(xmax) * 1.1
+    nsections = 5  # filter property
-fig = plt.figure()
+    designer = ThirdOctaveBankDesigner(fs)
-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
+    # Test
-    fu = designer.fu(x)  # upper
+    xmin = designer.nominal_txt_tox('25')  # lowest highest band, nominal name
+    xmax = designer.nominal_txt_tox('20k')  # highest band
 
-    fln = fl/fnyq  # lower cutoff, normalized [0, 1]
+    compliant = True  # remains True if all filters comply
-    fun = fu/fnyq  # upper
+    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
 
-    print(designer.nominal_txt(x))
+        fl = designer.fl(x)  # lower cut off
+        fu = designer.fu(x)  # upper
 
-    sos = designer.createSOSFilter(x)
+        fln = fl/fnyq  # lower cutoff, normalized [0, 1]
-    fnorm, h = sosfreqz(sos, worN=2**18, whole=False)
+        fun = fu/fnyq  # upper
-    freq = fnorm*fnyq/np.pi
 
-    h_dB = 20*np.log10(np.abs(h) + np.finfo(float).eps)
+        sos = designer.createSOSFilter(x)
-    if not designer.testStandardCompliance(x, freq, h_dB):
+        fnorm, h = sosfreqz(sos, worN=2**18, whole=False)
-        print('==============================')
+        freq = fnorm*fnyq/np.pi
-        print(designer.nominal_txt(x), ' does not comply')
-        compliant = False
 
-    plt.semilogx(freq, h_dB)
+        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\n')
+            compliant = False
 
-    # Upper / lower limits
+        plt.semilogx(freq, h_dB)
-    freqlim, ulim, llim = designer.band_limits(x)
-    plt.semilogx(freqlim, ulim)
-    plt.semilogx(freqlim, llim)
 
-    plt.ylim(-5, 1)
+        # Upper / lower limits
-    plt.xlim(fll, fuu)
+        freqlim, ulim, llim = designer.band_limits(x)
+        plt.semilogx(freqlim, ulim)
+        plt.semilogx(freqlim, llim)
 
-print(f"\nPassed test : {compliant}")
+        plt.ylim(-5, 1)
+        plt.xlim(fll, fuu)
+
+    print(f"Passed test : {compliant} for fs = {fs} Hz")
+    return compliant
+
+
+if __name__ == "__main__":
+    # Make list of fs to be tested
+    # Taken from https://en.wikipedia.org/wiki/Sampling_(signal_processing)
+    fs_all = [8000, 11025, 16000, 22050, 32000, 37800, 44056, 44100, 47250,
+              48000, 50000, 50400, 64000, 88200, 96000, 176400, 192000, 352800]
+    allcompliant = True
+    for fs in fs_all:
+        compliant = check_one(fs)
+        if not compliant:
+            allcompliant = False
+print(f"\n\nAll passed test : {allcompliant}")