From 90a5d8741949496432aa6cdc9dc30e025210a768 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Thu, 6 Oct 2022 21:34:33 +0200 Subject: [PATCH] Typo fix. Added 44.1 kHz as possible sampling frequency for computing filters. --- src/lasp/dsp/lasp_slm.cpp | 2 +- src/lasp/filter/filterbank_design.py | 162 +++++++++++++++------------ 2 files changed, 90 insertions(+), 74 deletions(-) diff --git a/src/lasp/dsp/lasp_slm.cpp b/src/lasp/dsp/lasp_slm.cpp index 6ed72a1..853cca5 100644 --- a/src/lasp/dsp/lasp_slm.cpp +++ b/src/lasp/dsp/lasp_slm.cpp @@ -67,7 +67,7 @@ std::vector> createBandPass(const dmat &coefs) { return bf; } us SLM::suggestedDownSamplingFac(const d fs,const d tau) { - if(tau<0) throw rte("Invalid time weightin time constant"); + if(tau<0) throw rte("Invalid time weighting time constant"); if(fs<=0) throw rte("Invalid sampling frequency"); // A reasonable 'framerate' for the sound level meter, based on the // filtering time constant. diff --git a/src/lasp/filter/filterbank_design.py b/src/lasp/filter/filterbank_design.py index 4c09d18..7773917 100644 --- a/src/lasp/filter/filterbank_design.py +++ b/src/lasp/filter/filterbank_design.py @@ -10,12 +10,14 @@ Resulting filters are supposed to be standard compliant. See test/octave_fir_test.py for a testing """ -from .fir_design import bandpass_fir_design, freqResponse as firFreqResponse -import numpy as np +import warnings +import numpy as np # For designing second-order sections from scipy.signal import butter +from .fir_design import bandpass_fir_design +from .fir_design import freqResponse as firFreqResponse __all__ = ['OctaveBankDesigner', 'ThirdOctaveBankDesigner'] @@ -34,7 +36,7 @@ class FilterBankDesigner: self.fs = fs # Constant G, according to standard - self.G = 10**(3/10) + self.G = 10**(3 / 10) # Reference frequency for all filter banks self.fr = 1000. @@ -61,11 +63,13 @@ class FilterBankDesigner: # Interpolate limites to frequency array as given llim_full = np.interp(freq, freqlim, llim, left=-np.inf, right=-np.inf) - ulim_full = np.interp(freq, freqlim, ulim, - left=ulim[0], right=ulim[-1]) + ulim_full = np.interp(freq, + freqlim, + ulim, + left=ulim[0], + right=ulim[-1]) - return bool(np.all(llim_full <= h_dB) and - np.all(ulim_full >= h_dB)) + return bool(np.all(llim_full <= h_dB) and np.all(ulim_full >= h_dB)) def band_limits(self, x, filter_class): raise NotImplementedError() @@ -81,7 +85,8 @@ class FilterBankDesigner: if self.nominal_txt(x) == nom_txt: return x raise ValueError( - f'Could not find a nominal frequency corresponding to {nom_txt}. Hint: use \'5k\' instead of \'5000\'.') + f'Could not find a nominal frequency corresponding to {nom_txt}. Hint: use \'5k\' instead of \'5000\'.' + ) def sanitize_input(self, input_): if isinstance(input_, int): @@ -94,12 +99,11 @@ class FilterBankDesigner: # This is the "code" to create an array xl = self.sanitize_input(input_[0]) xu = self.sanitize_input(input_[2]) - return np.asarray(list(range(xl, xu+1))) + return np.asarray(list(range(xl, xu + 1))) else: x = [self.sanitize_input(xi) for xi in input_] return np.asarray(x) - def getxs(self, nom_txt_start, nom_txt_end): """Returns a list of all filter designators, for given start end end nominal frequencies. @@ -113,7 +117,8 @@ class FilterBankDesigner: """ xstart = self.nominal_txt_tox(nom_txt_start) xend = self.nominal_txt_tox(nom_txt_end) - return list(range(xstart, xend+1)) + return list(range(xstart, xend + 1)) + def fm(self, x): """Returns the exact midband frequency of the bandpass filter. @@ -123,7 +128,7 @@ class FilterBankDesigner: x = self.sanitize_input(x) # Exact midband frequency - return self.G**(x/self.b)*self.fr + return self.G**(x / self.b) * self.fr def fl(self, x): """Returns the exact cut-on frequency of the bandpass filter. @@ -132,7 +137,7 @@ class FilterBankDesigner: x: Midband designator """ x = self.sanitize_input(x) - return self.fm(x)*self.G**(-1/(2*self.b)) + return self.fm(x) * self.G**(-1 / (2 * self.b)) def fu(self, x): """Returns the exact cut-off frequency of the bandpass filter. @@ -141,7 +146,7 @@ class FilterBankDesigner: x: Midband designator """ x = self.sanitize_input(x) - return self.fm(x)*self.G**(1/(2*self.b)) + return self.fm(x) * self.G**(1 / (2 * self.b)) def createFirFilter(self, x): """Create a FIR filter for band designator b and sampling frequency fs. @@ -155,8 +160,8 @@ class FilterBankDesigner: # For designing the filter, the lower and upper frequencies need to be # slightly adjusted to fall within the limits for a class 1 filter. - fl = self.fl(x)*self.firFac_l(x) - fu = self.fu(x)*self.firFac_u(x) + fl = self.fl(x) * self.firFac_l(x) + fu = self.fu(x) * self.firFac_u(x) return bandpass_fir_design(self.firFilterLength, fd, fl, fu) @@ -171,15 +176,15 @@ class FilterBankDesigner: SOS_ORDER = 5 fs = self.fs - fl = self.fl(x)*self.sosFac_l(x) - fu = self.fu(x)*self.sosFac_u(x) + fl = self.fl(x) * self.sosFac_l(x) + fu = self.fu(x) * self.sosFac_u(x) - fnyq = fs/2 + fnyq = fs / 2 # Normalized upper and lower frequencies of the bandpass - fl_n = fl/fnyq + fl_n = fl / fnyq x = self.sanitize_input(x) - fu_n = fu/fnyq + fu_n = fu / fnyq return butter(SOS_ORDER, [fl_n, fu_n], output='sos', btype='band') @@ -200,9 +205,11 @@ class FilterBankDesigner: return firFreqResponse(fd, freq, fir) - - def getNarrowBandFromOctaveBand(self, xl, xu, - levels_in_bands, npoints=500, + def getNarrowBandFromOctaveBand(self, + xl, + xu, + levels_in_bands, + npoints=500, method='flat', scale='lin'): """Create a narrow band spectrum based on a spectrum in (fractional) @@ -256,7 +263,7 @@ class FilterBankDesigner: power_cur = 10**(levels_in_bands[i] / 10) power_narrow = power_cur / indices_cur[0].size - level_narrow = 10*np.log10(power_narrow) + level_narrow = 10 * np.log10(power_narrow) levels_narrow[indices_cur] = level_narrow return freq, levels_narrow @@ -292,7 +299,7 @@ class FilterBankDesigner: levels_in_bands = [] nom_txt = [] - for x in range(xl, xu+1): + for x in range(xl, xu + 1): fl = self.fl(x) fu = self.fu(x) if x != xu: @@ -301,7 +308,7 @@ class FilterBankDesigner: indices_cur = np.where((freq >= fl) & (freq <= fu)) power_cur = np.sum(10**(levels_narrow[indices_cur] / 10)) - levels_in_bands.append(10*np.log10(power_cur)) + levels_in_bands.append(10 * np.log10(power_cur)) nom_txt.append(self.nominal_txt(x)) freq_in_bands.append(self.fm(x)) @@ -341,9 +348,9 @@ class OctaveBankDesigner(FilterBankDesigner): b = 1 # Exact midband frequency - fm = self.G**(x/self.b)*self.fr + fm = self.G**(x / self.b) * self.fr - G_power_values_pos = [0, 1/8, 1/4, 3/8, 1/2, 1/2, 1, 2, 3, 4] + G_power_values_pos = [0, 1 / 8, 1 / 4, 3 / 8, 1 / 2, 1 / 2, 1, 2, 3, 4] G_power_values_neg = [-i for i in G_power_values_pos] G_power_values_neg.reverse() G_power_values = G_power_values_neg[:-1] + G_power_values_pos @@ -351,41 +358,43 @@ class OctaveBankDesigner(FilterBankDesigner): mininf = -1e300 if filter_class == 1: - lower_limits_pos = [-0.3, -0.4, - - 0.6, -1.3, -5.0, -5.0] + 4*[mininf] + lower_limits_pos = [-0.3, -0.4, -0.6, -1.3, -5.0, -5.0 + ] + 4 * [mininf] elif filter_class == 0: - lower_limits_pos = [-0.15, -0.2, - - 0.4, -1.1, -4.5, -4.5] + 4*[mininf] + lower_limits_pos = [-0.15, -0.2, -0.4, -1.1, -4.5, -4.5 + ] + 4 * [mininf] lower_limits_neg = lower_limits_pos[:] lower_limits_neg.reverse() lower_limits = np.asarray(lower_limits_neg[:-1] + lower_limits_pos) if filter_class == 1: - upper_limits_pos = [0.3]*5 + [-2, -17.5, -42, -61, -70] + upper_limits_pos = [0.3] * 5 + [-2, -17.5, -42, -61, -70] if filter_class == 0: - upper_limits_pos = [0.15]*5 + [-2.3, -18, -42.5, -62, -75] + upper_limits_pos = [0.15] * 5 + [-2.3, -18, -42.5, -62, -75] upper_limits_neg = upper_limits_pos[:] upper_limits_neg.reverse() upper_limits = np.asarray(upper_limits_neg[:-1] + upper_limits_pos) - freqs = fm*self.G**np.asarray(G_power_values) + freqs = fm * self.G**np.asarray(G_power_values) return freqs, lower_limits, upper_limits def nominal_txt(self, x): """Returns textual repressentation of corresponding to the nominal frequency.""" - nominals = {4: '16k', - 3: '8k', - 2: '4k', - 1: '2k', - 0: '1k', - -1: '500', - -2: '250', - -3: '125', - -4: '63', - -5: '31.5', - -6: '16'} + nominals = { + 4: '16k', + 3: '8k', + 2: '4k', + 1: '2k', + 0: '1k', + -1: '500', + -2: '250', + -3: '125', + -4: '63', + -5: '31.5', + -6: '16' + } assert len(nominals) == len(self.xs) return nominals[x] @@ -397,7 +406,7 @@ class OctaveBankDesigner(FilterBankDesigner): return .995 elif x in (3, 1): return .99 - elif x in(-6, -4, -2, 2, 0): + elif x in (-6, -4, -2, 2, 0): return .98 else: return .96 @@ -472,16 +481,12 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): super().__init__(fs) self.xs = list(range(-16, 14)) # Text corresponding to the nominal frequency - self._nominal_txt = ['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'] + self._nominal_txt = [ + '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' + ] assert len(self.xs) == len(self._nominal_txt) @@ -499,8 +504,7 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): elif type(x) == list: index_start = x[0] - self.xs[0] index_stop = x[-1] - self.xs[0] - return self._nominal_txt[index_start:index_stop+1] - + return self._nominal_txt[index_start:index_stop + 1] def band_limits(self, x, filter_class=0): """Returns the third octave band filter limits for filter designator x. @@ -516,13 +520,17 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): in *deciBell*, upper limits in *deciBell*, respectively. """ - fm = self.G**(x/self.b)*self.fr + fm = self.G**(x / self.b) * self.fr plusinf = 20 - f_ratio_pos = [1., 1.02667, 1.05575, 1.08746, 1.12202, 1.12202, - 1.29437, 1.88173, 3.05365, 5.39195, plusinf] + f_ratio_pos = [ + 1., 1.02667, 1.05575, 1.08746, 1.12202, 1.12202, 1.29437, 1.88173, + 3.05365, 5.39195, plusinf + ] - f_ratio_neg = [0.97402, 0.94719, 0.91958, 0.89125, 0.89125, - 0.77257, 0.53143, 0.32748, 0.18546, 1/plusinf] + f_ratio_neg = [ + 0.97402, 0.94719, 0.91958, 0.89125, 0.89125, 0.77257, 0.53143, + 0.32748, 0.18546, 1 / plusinf + ] f_ratio_neg.reverse() f_ratio = f_ratio_neg + f_ratio_pos @@ -530,9 +538,9 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): mininf = -1e300 if filter_class == 1: - upper_limits_pos = [.3]*5 + [-2, -17.5, -42, -61, -70, -70] + upper_limits_pos = [.3] * 5 + [-2, -17.5, -42, -61, -70, -70] elif filter_class == 0: - upper_limits_pos = [.15]*5 + [-2.3, -18, -42.5, -62, -75, -75] + upper_limits_pos = [.15] * 5 + [-2.3, -18, -42.5, -62, -75, -75] else: raise ValueError('Filter class should either be 0 or 1') @@ -541,17 +549,21 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): upper_limits = np.array(upper_limits_neg[:-1] + upper_limits_pos) if filter_class == 1: - lower_limits_pos = [-.3, -.4, -.6, -1.3, -5, -5, mininf, mininf, - mininf, mininf, mininf] + lower_limits_pos = [ + -.3, -.4, -.6, -1.3, -5, -5, mininf, mininf, mininf, mininf, + mininf + ] elif filter_class == 0: - lower_limits_pos = [-.15, -.2, -.4, -1.1, -4.5, -4.5, mininf, mininf, - mininf, mininf, mininf] + lower_limits_pos = [ + -.15, -.2, -.4, -1.1, -4.5, -4.5, mininf, mininf, mininf, + mininf, mininf + ] lower_limits_neg = lower_limits_pos[:] lower_limits_neg.reverse() lower_limits = np.array(lower_limits_neg[:-1] + lower_limits_pos) - freqs = fm*np.array(f_ratio) + freqs = fm * np.array(f_ratio) return freqs, lower_limits, upper_limits @@ -592,13 +604,17 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): filter.""" # Idea: correct for frequency warping: if np.isclose(self.fs, 48000): + return 1.00 + elif np.isclose(self.fs, 41000): + warnings.warn( + f'Frequency {self.fs} might not result in correct filters') + return 1.00 elif np.isclose(self.fs, 32768): return 1.00 else: raise ValueError('Unimplemented sampling frequency for SOS' 'filter design') - def sosFac_u(self, x): """Right side percentage of change in cut-on frequency for designing the filter."""