diff --git a/src/lasp/filter/filterbank_design.py b/src/lasp/filter/filterbank_design.py index 25b0058..362ebbe 100644 --- a/src/lasp/filter/filterbank_design.py +++ b/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 diff --git a/test/scipy_sos_thirdoctavebank.py b/test/scipy_sos_thirdoctavebank.py index 7fc55a6..70413eb 100644 --- a/test/scipy_sos_thirdoctavebank.py +++ b/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 -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 +def check_one(fs): + """Check for fs = fs""" + # Create filter bank + nsections = 5 # filter property + designer = ThirdOctaveBankDesigner(fs) - fl = designer.fl(x) # lower cut off - fu = designer.fu(x) # upper + # Test + 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] - fun = fu/fnyq # upper + 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 - print(designer.nominal_txt(x)) + fl = designer.fl(x) # lower cut off + fu = designer.fu(x) # upper - sos = designer.createSOSFilter(x) - fnorm, h = sosfreqz(sos, worN=2**18, whole=False) - freq = fnorm*fnyq/np.pi + fln = fl/fnyq # lower cutoff, normalized [0, 1] + fun = fu/fnyq # upper - 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 + sos = designer.createSOSFilter(x) + fnorm, h = sosfreqz(sos, worN=2**18, whole=False) + freq = fnorm*fnyq/np.pi - 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 - freqlim, ulim, llim = designer.band_limits(x) - plt.semilogx(freqlim, ulim) - plt.semilogx(freqlim, llim) + plt.semilogx(freq, h_dB) - plt.ylim(-5, 1) - plt.xlim(fll, fuu) + # Upper / lower limits + 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}")