4Author: J.A. de Jong - ASCEE
6Description: FIR filter design for octave bands from 16Hz to 16 kHz for a
7sampling frequency of 48 kHz, filter design for one-third octave bands.
8Resulting filters are supposed to be standard compliant.
10See test/octave_fir_test.py for a testing
17from scipy.signal
import butter
19from .fir_design
import bandpass_fir_design
20from .fir_design
import freqResponse
as firFreqResponse
22__all__ = [
'OctaveBankDesigner',
'ThirdOctaveBankDesigner']
26 """A class responsible for designing FIR filters."""
29 """Initialize a filter bank designer.
32 fs: Sampling frequency [Hz]
45 """Test whether the filter with given frequency response is compliant
50 freq: Array of frequencies to test for. Note: *Should be fine
51 enough to follow response!*
52 h_dB: Filter frequency response in *deciBell*
53 filter_class: Filter class to test for
56 True if filter is norm-compliant, False if not
59 if np.isclose(freq[0], 0):
62 freqlim, llim, ulim = self.
band_limits(x, filter_class)
65 llim_full = np.interp(freq, freqlim, llim, left=-np.inf, right=-np.inf)
66 ulim_full = np.interp(freq,
72 return bool(np.all(llim_full <= h_dB)
and np.all(ulim_full >= h_dB))
75 raise NotImplementedError()
78 """Returns the x-value corresponding to a certain nominal txt: '1k' ->
82 nom_txt: Text-representation of midband frequency
85 if self.nominal_txt(x) == nom_txt:
88 f
'Could not find a nominal frequency corresponding to {nom_txt}. '
89 'Hint: use \'5k\' instead of \'5000\'.')
92 if isinstance(input_, int):
94 elif isinstance(input_, str):
97 elif isinstance(input_, list):
98 if len(input_) == 3
and input_[1]
is None:
102 return np.asarray(list(range(xl, xu + 1)))
107 def getxs(self, nom_txt_start, nom_txt_end):
108 """Returns a list of all filter designators, for given start end end
112 nom_txt_start: Start frequency band, i.e. '31.5'
113 nom_txt_end: End frequency band, i.e. '10k'
120 return list(range(xstart, xend + 1))
123 """Returns the exact midband frequency of the bandpass filter.
126 x: Midband designator
131 return self.
G**(x / self.
b) * self.
fr
134 """Returns the exact cut-on frequency of the bandpass filter.
137 x: Midband designator
140 return self.
fm(x) * self.
G**(-1 / (2 * self.
b))
143 """Returns the exact cut-off frequency of the bandpass filter.
146 x: Midband designator
149 return self.
fm(x) * self.
G**(1 / (2 * self.
b))
152 """Create a FIR filter for band designator b and sampling frequency fs.
153 firdecimation should be obtained from firdecimation() method.
156 filter: 1D ndarray with FIR filter coefficients
158 assert np.isclose(self.
fs, 48000),
"Invalid sampling frequency"
159 fd = self.
fs / np.prod(self.firDecimation(x))
163 fl = self.
fl(x) * self.firFac_l(x)
164 fu = self.
fu(x) * self.firFac_u(x)
169 """Create a Second Order Section filter (cascaded BiQuad's) for the
170 given sample rate and band designator.
178 fl = self.
fl(x) * self.sosFac_l(x)
179 fu = self.
fu(x) * self.sosFac_u(x)
190 return np.tile(np.asarray([0, 0, 0, 1, 0, 0]), (SOS_ORDER, 1))
194 highpass = butter(SOS_ORDER, fl_n, output=
'sos', btype=
'highpass')
195 allpass = np.tile(np.asarray([1, 0, 0, 1, 0, 0]),
196 (SOS_ORDER-highpass.shape[0], 1))
197 return np.vstack((highpass, allpass))
200 return butter(SOS_ORDER, [fl_n, fu_n], output=
'sos', btype=
'band')
203 """Compute the frequency response for a certain filter.
206 x: Midband designator
207 freq: Array of frequencies to evaluate on
210 h: Linear filter transfer function [-]
215 fd = fs / np.prod(self.firdecimation(x))
226 """Create a narrow band spectrum based on a spectrum in (fractional)
227 octave bands. The result is create such that the total energy in each
228 frequency band is constant. The latter can be checked by working it
229 back to (fractional) octave bands, which is doen using the function
230 `getOctaveBandsFromNarrowBand`. Note that the resulting narrow band has
231 units of *power*, not power spectral density. Input should be levels in
235 xl: Band designator of lowest band
236 xu: Band designator of highest band
237 levels_in_bands: levels in dB for each band, should have length
239 npoints: Number of discrete frequency points in linear frequency
241 method: 'flat' to give the same power for all frequencies in a
243 scale: 'lin' for a linear frequency distribution. 'log' for a
244 logspace. Use 'log' with caution and only if you really know what
248 freq, levels_dB. Where levels_dB is an array of narrow band levels
256 freq = np.linspace(fll, fuu, npoints)
258 freq = np.logspace(np.log10(fll), np.log10(fuu), npoints)
260 raise ValueError(f
'Invalid scale parameter: {scale}')
262 levels_narrow = np.empty_like(freq)
265 for i, x
in enumerate(range(xl, xu + 1)):
271 indices_cur = np.where((freq >= fl) & (freq < fu))
273 indices_cur = np.where((freq >= fl) & (freq <= fu))
275 power_cur = 10**(levels_in_bands[i] / 10)
276 power_narrow = power_cur / indices_cur[0].size
277 level_narrow = 10 * np.log10(power_narrow)
278 levels_narrow[indices_cur] = level_narrow
280 return freq, levels_narrow
282 raise ValueError(
'Unimplemented interpolation method')
285 """Put all level results in a certain frequency band (`binning`), by
286 summing up the power.
289 freq: Narrow-band frequency array
290 levels_narrow: Narrow band power levels, should be power, not *PSD*
293 (fm, levels_binned, nom_txt)
313 for x
in range(xl, xu + 1):
317 indices_cur = np.where((freq >= fl) & (freq < fu))
319 indices_cur = np.where((freq >= fl) & (freq <= fu))
321 power_cur = np.sum(10**(levels_narrow[indices_cur] / 10))
322 levels_in_bands.append(10 * np.log10(power_cur))
324 nom_txt.append(self.nominal_txt(x))
325 freq_in_bands.append(self.
fm(x))
327 return np.array(freq_in_bands), np.array(levels_in_bands), nom_txt
331 """Octave band filter designer."""
343 """All possible band designators for an octave band filter."""
344 return list(range(-7, 5))
347 """Returns the octave band filter limits for filter designator x.
350 x: Filter offset power from the reference frequency of 1000 Hz.
351 filter_class: Either 0 or 1, defines the tolerances on the
355 freq, llim, ulim: Tuple of Numpy arrays containing the frequencies
356 of the corner points of the filter frequency response limits, lower
357 limits in *deciBell*, upper limits in *deciBell*, respectively.
362 fm = self.
G**(x / self.
bb) * self.
fr
364 G_power_values_pos = [0, 1 / 8, 1 / 4, 3 / 8, 1 / 2, 1 / 2, 1, 2, 3, 4]
365 G_power_values_neg = [-i
for i
in G_power_values_pos]
366 G_power_values_neg.reverse()
367 G_power_values = G_power_values_neg[:-1] + G_power_values_pos
371 if filter_class == 1:
372 lower_limits_pos = [-0.3, -0.4, -0.6, -1.3, -5.0, -5.0
374 elif filter_class == 0:
375 lower_limits_pos = [-0.15, -0.2, -0.4, -1.1, -4.5, -4.5
377 lower_limits_neg = lower_limits_pos[:]
378 lower_limits_neg.reverse()
379 lower_limits = np.asarray(lower_limits_neg[:-1] + lower_limits_pos)
381 if filter_class == 1:
382 upper_limits_pos = [0.3] * 5 + [-2, -17.5, -42, -61, -70]
383 if filter_class == 0:
384 upper_limits_pos = [0.15] * 5 + [-2.3, -18, -42.5, -62, -75]
385 upper_limits_neg = upper_limits_pos[:]
386 upper_limits_neg.reverse()
387 upper_limits = np.asarray(upper_limits_neg[:-1] + upper_limits_pos)
389 freqs = fm * self.
G**np.asarray(G_power_values)
391 return freqs, lower_limits, upper_limits
394 """Returns textual repressentation of corresponding to the nominal
410 assert len(nominals) == len(self.
xs)
414 """Factor with which to multiply the cut-on frequency of the FIR
416 assert int(self.
fsfs) == 48000,
'Fir coefs are only valid for 48kHz fs'
421 elif x
in (-6, -4, -2, 2, 0):
427 """Factor with which to multiply the cut-off frequency of the FIR
429 assert int(self.
fsfs) == 48000,
'Fir coefs are only valid for 48kHz fs'
434 elif x
in (-6, -4, -2, 2, 0):
440 """Required firdecimation for each filter."""
441 assert int(self.
fsfs) == 48000,
'Fir coefs are only valid for 48kHz fs'
452 assert False,
'Overlooked firdecimation'
455 """Left side percentage of change in cut-on frequency for designing the
456 filter, for OCTAVE band filter.
459 x: Filter band designator
462 if int(self.
fsfs)
in [44100, 48000, 96000]:
468 """Right side percentage of change in cut-on frequency for designing
472 x: Filter band designator
474 if int(self.
fsfs)
in [44100, 48000, 96000]:
484 Initialize ThirdOctaveBankDesigner, a filter bank designer for
485 one-third octave bands
488 fs: Sampling frequency in [Hz] - can be None
491 self.
xs = list(range(-19, 14))
494 '25',
'31.5',
'40',
'50',
'63',
'80',
'100',
'125',
'160',
'200',
495 '250',
'315',
'400',
'500',
'630',
'800',
'1k',
'1.25k',
'1.6k',
496 '2k',
'2.5k',
'3.15k',
'4k',
'5k',
'6.3k',
'8k',
'10k',
'12.5k',
510 index = x - self.
xs[0]
512 elif type(x) == list:
513 index_start = x[0] - self.
xs[0]
514 index_stop = x[-1] - self.
xs[0]
518 """Returns the third octave band filter limits for filter designator x.
521 x: Filter offset power from the reference frequency of 1000 Hz.
522 filter_class: Either 0 or 1, defines the tolerances on the
526 freq, llim, ulim: Tuple of Numpy arrays containing the frequencies
527 of the corner points of the filter frequency response limits, lower
528 limits in *deciBell*, upper limits in *deciBell*, respectively.
531 fm = self.
G**(x / self.
bb) * self.
fr
534 1., 1.02667, 1.05575, 1.08746, 1.12202, 1.12202, 1.29437, 1.88173,
535 3.05365, 5.39195, plusinf
539 0.97402, 0.94719, 0.91958, 0.89125, 0.89125, 0.77257, 0.53143,
540 0.32748, 0.18546, 1 / plusinf
542 f_ratio_neg.reverse()
544 f_ratio = f_ratio_neg + f_ratio_pos
548 if filter_class == 1:
549 upper_limits_pos = [.3] * 5 + [-2, -17.5, -42, -61, -70, -70]
550 elif filter_class == 0:
551 upper_limits_pos = [.15] * 5 + [-2.3, -18, -42.5, -62, -75, -75]
553 raise ValueError(
'Filter class should either be 0 or 1')
555 upper_limits_neg = upper_limits_pos[:]
556 upper_limits_neg.reverse()
557 upper_limits = np.array(upper_limits_neg[:-1] + upper_limits_pos)
559 if filter_class == 1:
561 -.3, -.4, -.6, -1.3, -5, -5, mininf, mininf, mininf, mininf,
564 elif filter_class == 0:
566 -.15, -.2, -.4, -1.1, -4.5, -4.5, mininf, mininf, mininf,
570 lower_limits_neg = lower_limits_pos[:]
571 lower_limits_neg.reverse()
572 lower_limits = np.array(lower_limits_neg[:-1] + lower_limits_pos)
574 freqs = fm * np.array(f_ratio)
576 return freqs, lower_limits, upper_limits
579 assert int(self.
fsfs) == 48000,
'Fir coefs are only valid for 48kHz fs'
590 assert False,
'Bug: overlooked firdecimation'
593 if x
in (-13, -7, -1, 5, 11, 12, 13):
595 elif x
in (-12, -6, 0, 6):
601 if x
in (-14, -13, -8, -7, -1, -2, 3, 4, 5, 10, 11, 12):
605 elif x
in (12, -6, 0):
611 """Left side percentage of change in cut-on frequency for designing the
614 if int(self.
fsfs)
in [48000]:
620 """Right side percentage of change in cut-on frequency for designing
622 if int(self.
fsfs)
in [48000]:
A class responsible for designing FIR filters.
fl(self, x)
Returns the exact cut-on frequency of the bandpass filter.
fm(self, x)
Returns the exact midband frequency of the bandpass filter.
firFreqResponse(self, x, freq)
Compute the frequency response for a certain filter.
createFirFilter(self, x)
Create a FIR filter for band designator b and sampling frequency fs.
fu(self, x)
Returns the exact cut-off frequency of the bandpass filter.
getxs(self, nom_txt_start, nom_txt_end)
Returns a list of all filter designators, for given start end end nominal frequencies.
__init__(self, fs)
Initialize a filter bank designer.
nominal_txt_tox(self, str nom_txt)
Returns the x-value corresponding to a certain nominal txt: '1k' -> 0.
testStandardCompliance(self, x, freq, h_dB, filter_class=0)
Test whether the filter with given frequency response is compliant with the standard.
getOctaveBandFromNarrowBand(self, freq, levels_narrow)
Put all level results in a certain frequency band (binning), by summing up the power.
band_limits(self, x, filter_class)
getNarrowBandFromOctaveBand(self, xl, xu, levels_in_bands, npoints=500, method='flat', scale='lin')
Create a narrow band spectrum based on a spectrum in (fractional) octave bands.
createSOSFilter(self, int x)
Create a Second Order Section filter (cascaded BiQuad's) for the given sample rate and band designato...
sanitize_input(self, input_)
Octave band filter designer.
firFac_l(self, x)
Factor with which to multiply the cut-on frequency of the FIR filter.
sosFac_u(self, x)
Right side percentage of change in cut-on frequency for designing the filter.
sosFac_l(self, x)
Left side percentage of change in cut-on frequency for designing the filter, for OCTAVE band filter.
firDecimation(self, x)
Required firdecimation for each filter.
band_limits(self, x, filter_class=0)
Returns the octave band filter limits for filter designator x.
__init__(self, fs)
Initialize a filter bank designer.
xs(self)
All possible band designators for an octave band filter.
firFac_u(self, x)
Factor with which to multiply the cut-off frequency of the FIR filter.
nominal_txt(self, x)
Returns textual repressentation of corresponding to the nominal frequency.
sosFac_u(self, x)
Right side percentage of change in cut-on frequency for designing the filter.
sosFac_l(self, x)
Left side percentage of change in cut-on frequency for designing the filter.
__init__(self, float fs)
Initialize ThirdOctaveBankDesigner, a filter bank designer for one-third octave bands.
band_limits(self, x, filter_class=0)
Returns the third octave band filter limits for filter designator x.