From 786afc3fac57dbb41aa9ca7e07eb30d3ffa381dc Mon Sep 17 00:00:00 2001 From: Casper Date: Fri, 6 Jan 2023 09:44:34 +0100 Subject: [PATCH] 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 | 50 +++++++++++----------- test/scipy_sos_thirdoctavebank.py | 63 ++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 24 deletions(-) create mode 100644 test/scipy_sos_thirdoctavebank.py diff --git a/src/lasp/filter/filterbank_design.py b/src/lasp/filter/filterbank_design.py index 5370725..25b0058 100644 --- a/src/lasp/filter/filterbank_design.py +++ b/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 diff --git a/test/scipy_sos_thirdoctavebank.py b/test/scipy_sos_thirdoctavebank.py new file mode 100644 index 0000000..7fc55a6 --- /dev/null +++ b/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}")