LASP 1.0
Library for Acoustic Signal Processing
Loading...
Searching...
No Matches
filterbank_design.py
Go to the documentation of this file.
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""!
4Author: J.A. de Jong - ASCEE
5
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.
9
10See test/octave_fir_test.py for a testing
11
12"""
13import warnings
14
15import numpy as np
16# For designing second-order sections
17from scipy.signal import butter
18
19from .fir_design import bandpass_fir_design
20from .fir_design import freqResponse as firFreqResponse
21
22__all__ = ['OctaveBankDesigner', 'ThirdOctaveBankDesigner']
23
24
26 """A class responsible for designing FIR filters."""
27
28 def __init__(self, fs):
29 """Initialize a filter bank designer.
30
31 Args:
32 fs: Sampling frequency [Hz]
33 """
34 # Default FIR filter length
35 self.firFilterLength = 256 # Filter order
36 self.fs = fs
37
38 # Constant G, according to standard
39 self.G = 10**(3 / 10)
40
41 # Reference frequency for all filter banks
42 self.fr = 1000.
43
44 def testStandardCompliance(self, x, freq, h_dB, filter_class=0):
45 """Test whether the filter with given frequency response is compliant
46 with the standard.
47
48 Args:
49 x: Band designator
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
54
55 Returns:
56 True if filter is norm-compliant, False if not
57 """
58 # Skip zero-frequency
59 if np.isclose(freq[0], 0):
60 freq = freq[1:]
61 h_dB = h_dB[1:]
62 freqlim, llim, ulim = self.band_limits(x, filter_class)
63
64 # Interpolate limites to frequency array as given
65 llim_full = np.interp(freq, freqlim, llim, left=-np.inf, right=-np.inf)
66 ulim_full = np.interp(freq,
67 freqlim,
68 ulim,
69 left=ulim[0],
70 right=ulim[-1])
71
72 return bool(np.all(llim_full <= h_dB) and np.all(ulim_full >= h_dB))
73
74 def band_limits(self, x, filter_class):
75 raise NotImplementedError()
76
77 def nominal_txt_tox(self, nom_txt: str):
78 """Returns the x-value corresponding to a certain nominal txt: '1k' ->
79 0.
80
81 Args:
82 nom_txt: Text-representation of midband frequency
83 """
84 for x in self.xs:
85 if self.nominal_txt(x) == nom_txt:
86 return x
87 raise ValueError(
88 f'Could not find a nominal frequency corresponding to {nom_txt}. '
89 'Hint: use \'5k\' instead of \'5000\'.')
90
91 def sanitize_input(self, input_):
92 if isinstance(input_, int):
93 return input_
94 elif isinstance(input_, str):
95 return self.nominal_txt_tox(input_)
96
97 elif isinstance(input_, list):
98 if len(input_) == 3 and input_[1] is None:
99 # This is the "code" to create an array
100 xl = self.sanitize_input(input_[0])
101 xu = self.sanitize_input(input_[2])
102 return np.asarray(list(range(xl, xu + 1)))
103 else:
104 x = [self.sanitize_input(xi) for xi in input_]
105 return np.asarray(x)
106
107 def getxs(self, nom_txt_start, nom_txt_end):
108 """Returns a list of all filter designators, for given start end end
109 nominal frequencies.
110
111 Args:
112 nom_txt_start: Start frequency band, i.e. '31.5'
113 nom_txt_end: End frequency band, i.e. '10k'
114
115 Returns:
116 [x0, x1, ..]
117 """
118 xstart = self.nominal_txt_tox(nom_txt_start)
119 xend = self.nominal_txt_tox(nom_txt_end)
120 return list(range(xstart, xend + 1))
121
122 def fm(self, x):
123 """Returns the exact midband frequency of the bandpass filter.
124
125 Args:
126 x: Midband designator
127 """
128 x = self.sanitize_input(x)
129
130 # Exact midband frequency
131 return self.G**(x / self.b) * self.fr
132
133 def fl(self, x):
134 """Returns the exact cut-on frequency of the bandpass filter.
135
136 Args:
137 x: Midband designator
138 """
139 x = self.sanitize_input(x)
140 return self.fm(x) * self.G**(-1 / (2 * self.b))
141
142 def fu(self, x):
143 """Returns the exact cut-off frequency of the bandpass filter.
144
145 Args:
146 x: Midband designator
147 """
148 x = self.sanitize_input(x)
149 return self.fm(x) * self.G**(1 / (2 * self.b))
150
151 def createFirFilter(self, x):
152 """Create a FIR filter for band designator b and sampling frequency fs.
153 firdecimation should be obtained from firdecimation() method.
154
155 Returns:
156 filter: 1D ndarray with FIR filter coefficients
157 """
158 assert np.isclose(self.fs, 48000), "Invalid sampling frequency"
159 fd = self.fs / np.prod(self.firDecimation(x))
160
161 # For designing the filter, the lower and upper frequencies need to be
162 # slightly adjusted to fall within the limits for a class 1 filter.
163 fl = self.fl(x) * self.firFac_l(x)
164 fu = self.fu(x) * self.firFac_u(x)
165
166 return bandpass_fir_design(self.firFilterLength, fd, fl, fu)
167
168 def createSOSFilter(self, x: int):
169 """Create a Second Order Section filter (cascaded BiQuad's) for the
170 given sample rate and band designator.
171
172 Args:
173 x: Band designator
174 """
175 SOS_ORDER = 5
176
177 fs = self.fs
178 fl = self.fl(x) * self.sosFac_l(x)
179 fu = self.fu(x) * self.sosFac_u(x)
180
181 fnyq = fs / 2
182
183 # Normalized upper and lower frequencies of the bandpass
184 fl_n = fl / fnyq
185 x = self.sanitize_input(x)
186 fu_n = fu / fnyq
187
188 # Handle fl & fu >= fNyquist: return 'no pass'
189 if fl_n >= 1:
190 return np.tile(np.asarray([0, 0, 0, 1, 0, 0]), (SOS_ORDER, 1))
191
192 # Handle fu >= fNyquist: return high pass
193 if fu_n >= 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)) # same shape=(SOS_ORDER, 6)
198
199 # Regular bands
200 return butter(SOS_ORDER, [fl_n, fu_n], output='sos', btype='band')
201
202 def firFreqResponse(self, x, freq):
203 """Compute the frequency response for a certain filter.
204
205 Args:
206 x: Midband designator
207 freq: Array of frequencies to evaluate on
208
209 Returns:
210 h: Linear filter transfer function [-]
211 """
212 fir = self.createFirFilter(fs, x)
213
214 # Decimated sampling frequency [Hz]
215 fd = fs / np.prod(self.firdecimation(x))
216
217 return firFreqResponse(fd, freq, fir)
218
220 xl,
221 xu,
222 levels_in_bands,
223 npoints=500,
224 method='flat',
225 scale='lin'):
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
232 **deciBells**.
233
234 Args:
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
238 (xu+1 - xl)
239 npoints: Number of discrete frequency points in linear frequency
240 array.
241 method: 'flat' to give the same power for all frequencies in a
242 certain band.
243 scale: 'lin' for a linear frequency distribution. 'log' for a
244 logspace. Use 'log' with caution and only if you really know what
245 you are doing!'
246
247 Returns:
248 freq, levels_dB. Where levels_dB is an array of narrow band levels
249 """
250 # Lowest frequency of all frequencies
251 fll = self.fl(xl)
252 # Highest frequency of all frequencies
253 fuu = self.fu(xu)
254
255 if scale == 'lin':
256 freq = np.linspace(fll, fuu, npoints)
257 elif scale == 'log':
258 freq = np.logspace(np.log10(fll), np.log10(fuu), npoints)
259 else:
260 raise ValueError(f'Invalid scale parameter: {scale}')
261
262 levels_narrow = np.empty_like(freq)
263
264 if method == 'flat':
265 for i, x in enumerate(range(xl, xu + 1)):
266 fl = self.fl(x)
267 fu = self.fu(x)
268 # Find the indices in the frequency array which correspond to
269 # the frequency band x
270 if x != xu:
271 indices_cur = np.where((freq >= fl) & (freq < fu))
272 else:
273 indices_cur = np.where((freq >= fl) & (freq <= fu))
274
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
279
280 return freq, levels_narrow
281 else:
282 raise ValueError('Unimplemented interpolation method')
283
284 def getOctaveBandFromNarrowBand(self, freq, levels_narrow):
285 """Put all level results in a certain frequency band (`binning`), by
286 summing up the power.
287
288 Args:
289 freq: Narrow-band frequency array
290 levels_narrow: Narrow band power levels, should be power, not *PSD*
291
292 Returns:
293 (fm, levels_binned, nom_txt)
294 """
295 # Find lower frequency xl
296 for x in self.xs:
297 xl = x
298 fl = self.fl(x)
299 if fl >= freq[0]:
300 break
301
302 # Find upper frequency xu
303 for x in self.xs:
304 xu = x
305 fu = self.fu(xu)
306 if fu >= freq[-1]:
307 break
308
309 freq_in_bands = []
310 levels_in_bands = []
311 nom_txt = []
312
313 for x in range(xl, xu + 1):
314 fl = self.fl(x)
315 fu = self.fu(x)
316 if x != xu:
317 indices_cur = np.where((freq >= fl) & (freq < fu))
318 else:
319 indices_cur = np.where((freq >= fl) & (freq <= fu))
320
321 power_cur = np.sum(10**(levels_narrow[indices_cur] / 10))
322 levels_in_bands.append(10 * np.log10(power_cur))
323
324 nom_txt.append(self.nominal_txt(x))
325 freq_in_bands.append(self.fm(x))
326
327 return np.array(freq_in_bands), np.array(levels_in_bands), nom_txt
328
329
331 """Octave band filter designer."""
332
333 def __init__(self, fs):
334 super().__init__(fs)
335
336 @property
337 def b(self):
338 # Band division, 1 for octave bands
339 return 1
340
341 @property
342 def xs(self):
343 """All possible band designators for an octave band filter."""
344 return list(range(-7, 5))
345
346 def band_limits(self, x, filter_class=0):
347 """Returns the octave band filter limits for filter designator x.
348
349 Args:
350 x: Filter offset power from the reference frequency of 1000 Hz.
351 filter_class: Either 0 or 1, defines the tolerances on the
352 frequency response
353
354 Returns:
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.
358 """
359 b = 1
360
361 # Exact midband frequency
362 fm = self.G**(x / self.bb) * self.fr
363
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
368
369 mininf = -1e300
370
371 if filter_class == 1:
372 lower_limits_pos = [-0.3, -0.4, -0.6, -1.3, -5.0, -5.0
373 ] + 4 * [mininf]
374 elif filter_class == 0:
375 lower_limits_pos = [-0.15, -0.2, -0.4, -1.1, -4.5, -4.5
376 ] + 4 * [mininf]
377 lower_limits_neg = lower_limits_pos[:]
378 lower_limits_neg.reverse()
379 lower_limits = np.asarray(lower_limits_neg[:-1] + lower_limits_pos)
380
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)
388
389 freqs = fm * self.G**np.asarray(G_power_values)
390
391 return freqs, lower_limits, upper_limits
392
393 def nominal_txt(self, x):
394 """Returns textual repressentation of corresponding to the nominal
395 frequency."""
396 nominals = {
397 4: '16k',
398 3: '8k',
399 2: '4k',
400 1: '2k',
401 0: '1k',
402 -1: '500',
403 -2: '250',
404 -3: '125',
405 -4: '63',
406 -5: '31.5',
407 -6: '16',
408 -7: '8'
409 }
410 assert len(nominals) == len(self.xs)
411 return nominals[x]
412
413 def firFac_l(self, x):
414 """Factor with which to multiply the cut-on frequency of the FIR
415 filter."""
416 assert int(self.fsfs) == 48000, 'Fir coefs are only valid for 48kHz fs'
417 if x == 4:
418 return .995
419 elif x in (3, 1):
420 return .99
421 elif x in (-6, -4, -2, 2, 0):
422 return .98
423 else:
424 return .96
425
426 def firFac_u(self, x):
427 """Factor with which to multiply the cut-off frequency of the FIR
428 filter."""
429 assert int(self.fsfs) == 48000, 'Fir coefs are only valid for 48kHz fs'
430 if x == 4:
431 return 1.004
432 elif x in (3, 1):
433 return 1.006
434 elif x in (-6, -4, -2, 2, 0):
435 return 1.01
436 else:
437 return 1.02
438
439 def firDecimation(self, x):
440 """Required firdecimation for each filter."""
441 assert int(self.fsfs) == 48000, 'Fir coefs are only valid for 48kHz fs'
442 if x > 1:
443 return [1]
444 elif x > -2:
445 return [4]
446 elif x > -4:
447 return [4, 4]
448 elif x > -6:
449 return [4, 4, 4]
450 elif x == -6:
451 return [4, 4, 4, 4]
452 assert False, 'Overlooked firdecimation'
453
454 def sosFac_l(self, x):
455 """Left side percentage of change in cut-on frequency for designing the
456 filter, for OCTAVE band filter.
457
458 Args:
459 x: Filter band designator
460 """
461 # Idea: correct for frequency warping:
462 if int(self.fsfs) in [44100, 48000, 96000]:
463 return 1.0
464 else:
465 return 1.0
466
467 def sosFac_u(self, x):
468 """Right side percentage of change in cut-on frequency for designing
469 the filter.
470
471 Args:
472 x: Filter band designator
473 """
474 if int(self.fsfs) in [44100, 48000, 96000]:
475 return 1.0
476 else:
477 return 1.0
478
479
481
482 def __init__(self, fs: float):
483 """
484 Initialize ThirdOctaveBankDesigner, a filter bank designer for
485 one-third octave bands
486
487 Args:
488 fs: Sampling frequency in [Hz] - can be None
489 """
490 super().__init__(fs)
491 self.xs = list(range(-19, 14))
492 # Text corresponding to the nominal frequency
493 self._nominal_txt = ['12.5', '16', '20',
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',
497 '16k', '20k']
498
499 assert len(self.xs) == len(self._nominal_txt)
500
501 @property
502 def b(self):
503 # Band division factor, 3 for one-third octave bands
504 return 3
505
506 def nominal_txt(self, x):
507 # Put the nominal frequencies in a dictionary for easy access with
508 # x as the key.
509 if type(x) == int:
510 index = x - self.xs[0]
511 return self._nominal_txt[index]
512 elif type(x) == list:
513 index_start = x[0] - self.xs[0]
514 index_stop = x[-1] - self.xs[0]
515 return self._nominal_txt[index_start:index_stop + 1]
516
517 def band_limits(self, x, filter_class=0):
518 """Returns the third octave band filter limits for filter designator x.
519
520 Args:
521 x: Filter offset power from the reference frequency of 1000 Hz.
522 filter_class: Either 0 or 1, defines the tolerances on the
523 frequency response
524
525 Returns:
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.
529 """
530
531 fm = self.G**(x / self.bb) * self.fr
532 plusinf = 20
533 f_ratio_pos = [
534 1., 1.02667, 1.05575, 1.08746, 1.12202, 1.12202, 1.29437, 1.88173,
535 3.05365, 5.39195, plusinf
536 ]
537
538 f_ratio_neg = [
539 0.97402, 0.94719, 0.91958, 0.89125, 0.89125, 0.77257, 0.53143,
540 0.32748, 0.18546, 1 / plusinf
541 ]
542 f_ratio_neg.reverse()
543
544 f_ratio = f_ratio_neg + f_ratio_pos
545
546 mininf = -1e300
547
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]
552 else:
553 raise ValueError('Filter class should either be 0 or 1')
554
555 upper_limits_neg = upper_limits_pos[:]
556 upper_limits_neg.reverse()
557 upper_limits = np.array(upper_limits_neg[:-1] + upper_limits_pos)
558
559 if filter_class == 1:
560 lower_limits_pos = [
561 -.3, -.4, -.6, -1.3, -5, -5, mininf, mininf, mininf, mininf,
562 mininf
563 ]
564 elif filter_class == 0:
565 lower_limits_pos = [
566 -.15, -.2, -.4, -1.1, -4.5, -4.5, mininf, mininf, mininf,
567 mininf, mininf
568 ]
569
570 lower_limits_neg = lower_limits_pos[:]
571 lower_limits_neg.reverse()
572 lower_limits = np.array(lower_limits_neg[:-1] + lower_limits_pos)
573
574 freqs = fm * np.array(f_ratio)
575
576 return freqs, lower_limits, upper_limits
577
578 def firDecimation(self, x):
579 assert int(self.fsfs) == 48000, 'Fir coefs are only valid for 48kHz fs'
580 if x > 5:
581 return [1]
582 elif x > -1:
583 return [4]
584 elif x > -7:
585 return [4, 4]
586 elif x > -13:
587 return [4, 4, 4]
588 elif x > -17:
589 return [4, 4, 4, 4]
590 assert False, 'Bug: overlooked firdecimation'
591
592 def firFac_l(self, x):
593 if x in (-13, -7, -1, 5, 11, 12, 13):
594 return .995
595 elif x in (-12, -6, 0, 6):
596 return .98
597 else:
598 return .99
599
600 def firFac_u(self, x):
601 if x in (-14, -13, -8, -7, -1, -2, 3, 4, 5, 10, 11, 12):
602 return 1.005
603 elif x in (12, 13):
604 return 1.003
605 elif x in (12, -6, 0):
606 return 1.015
607 else:
608 return 1.01
609
610 def sosFac_l(self, x):
611 """Left side percentage of change in cut-on frequency for designing the
612 filter."""
613 # Idea: correct for frequency warping:
614 if int(self.fsfs) in [48000]:
615 return 1.0
616 else:
617 return 1.0
618
619 def sosFac_u(self, x):
620 """Right side percentage of change in cut-on frequency for designing
621 the filter."""
622 if int(self.fsfs) in [48000]:
623 return 1.0
624 else:
625 return 1.0
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.
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...
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.