Major cleanup and improvements of filter design code.

This commit is contained in:
Anne de Jong 2020-01-09 11:06:57 +01:00
parent 54173b6ecc
commit 00bd30eb2a
6 changed files with 472 additions and 213 deletions

View File

@ -1,5 +1,4 @@
from .soundpressureweighting import * from .soundpressureweighting import *
from .filterbank_design import OctaveBankDesigner, ThirdOctaveBankDesigner from .filterbank_design import *
from .filterbank_standard_limits import octave_band_limits, third_octave_band_limits
from .fir_design import * from .fir_design import *

View File

@ -4,28 +4,71 @@
Author: J.A. de Jong - ASCEE Author: J.A. de Jong - ASCEE
Description: FIR filter design for octave bands from 16Hz to 16 kHz for a Description: FIR filter design for octave bands from 16Hz to 16 kHz for a
sampling frequency of 48 kHz, FIR filter design for one-third octave bands sampling frequency of 48 kHz, filter design for one-third octave bands.
Resulting filters are supposed to be standard compliant.
See test/octave_fir_test.py for a testing See test/octave_fir_test.py for a testing
""" """
from .fir_design import bandpass_fir_design, freqResponse as frsp from .fir_design import bandpass_fir_design, freqResponse as firFreqResponse
import numpy as np import numpy as np
# For designing second-order sections
from scipy.signal import butter
__all__ = ['OctaveBankDesigner', 'ThirdOctaveBankDesigner'] __all__ = ['OctaveBankDesigner', 'ThirdOctaveBankDesigner']
class FilterBankDesigner: class FilterBankDesigner:
""" """A class responsible for designing FIR filters."""
A class responsible for designing FIR filters
""" def __init__(self, fs):
G = 10**(3/10) """Initialize a filter bank designer.
fr = 1000.
L = 256 # Filter order Args:
fs = 48000. # Sampling frequency fs: Sampling frequency [Hz]
"""
# Default FIR filter length
firFilterLength = 256 # Filter order
self.fs = fs
# Constant G, according to standard
self.G = 10**(3/10)
# Reference frequency for all filter banks
self.fr = 1000.
def testStandardCompliance(self, x, freq, h_dB, filter_class=0):
"""Test whether the filter with given frequency response is compliant
with the standard.
Args:
x: Band designator
freq: Array of frequencies to test for. Note: *Should be fine
enough to follow response!*
h_dB: Filter frequency response in *deciBell*
filter_class: Filter class to test for
Returns:
True if filter is norm-compliant, False if not
"""
# Skip zero-frequency
if np.isclose(freq[0], 0):
freq = freq[1:]
h_dB = h_dB[1:]
freqlim, llim, ulim = self.band_limits(x, filter_class)
# 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])
return bool(np.all(llim_full <= h_dB) and
np.all(ulim_full >= h_dB))
def fm(self, x): def fm(self, x):
""" """Returns the exact midband frequency of the bandpass filter.
Returns the exact midband frequency of the bandpass filter
Args: Args:
x: Midband designator x: Midband designator
@ -34,8 +77,7 @@ class FilterBankDesigner:
return self.G**(x/self.b)*self.fr return self.G**(x/self.b)*self.fr
def fl(self, x): def fl(self, x):
""" """Returns the exact cut-on frequency of the bandpass filter.
Returns the exact cut-on frequency of the bandpass filter
Args: Args:
x: Midband designator x: Midband designator
@ -43,50 +85,75 @@ class FilterBankDesigner:
return self.fm(x)*self.G**(-1/(2*self.b)) return self.fm(x)*self.G**(-1/(2*self.b))
def fu(self, x): def fu(self, x):
""" """Returns the exact cut-off frequency of the bandpass filter.
Returns the exact cut-off frequency of the bandpass filter
Args: Args:
x: Midband designator x: Midband designator
""" """
return self.fm(x)*self.G**(1/(2*self.b)) return self.fm(x)*self.G**(1/(2*self.b))
def createFirFilter(self, fs, x): def createFirFilter(self, x):
""" """Create a FIR filter for band designator b and sampling frequency fs.
Create a FIR filter for band designator b and sampling frequency fs. firdecimation should be obtained from firdecimation() method.
Decimation should be obtained from decimation() method.
Returns: Returns:
filter: 1D ndarray with FIR filter coefficients filter: 1D ndarray with FIR filter coefficients
""" """
assert np.isclose(fs, self.fs), "Invalid sampling frequency" assert np.isclose(fs, self.fs), "Invalid sampling frequency"
fd = fs / np.prod(self.decimation(x)) fd = fs / np.prod(self.firDecimation(x))
# For designing the filter, the lower and upper frequencies need to be # For designing the filter, the lower and upper frequencies need to be
# slightly adjusted to fall within the limits for a class 1 filter. # slightly adjusted to fall within the limits for a class 1 filter.
fl = self.fl(x)*self.fac_l(x) fl = self.fl(x)*self.firFac_l(x)
fu = self.fu(x)*self.fac_u(x) fu = self.fu(x)*self.firFac_u(x)
return bandpass_fir_design(self.L, fd, fl, fu) return bandpass_fir_design(self.firFilterLength, fd, fl, fu)
def freqResponse(self, fs, x, freq): def createSOSFilter(self, x: int):
""" """Create a Second Order Section filter (cascaded BiQuad's) for the
Compute the frequency response for a certain filter given sample rate and band designator.
Args: Args:
fs: Sampling frequency [Hz] x: Band designator
x: Midband designator
""" """
fir = self.createFilter(fs, x)
SOS_ORDER = 5
fs = self.fs
fl = self.fl(x)*self.sosFac_l(x)
fu = self.fu(x)*self.sosFac_u(x)
fnyq = fs/2
# Normalized upper and lower frequencies of the bandpass
fl_n = fl/fnyq
fu_n = fu/fnyq
return butter(SOS_ORDER, [fl_n, fu_n], output='sos', btype='band')
def firFreqResponse(self, x, freq):
"""Compute the frequency response for a certain filter.
Args:
x: Midband designator
freq: Array of frequencies to evaluate on
Returns:
h: Linear filter transfer function [-]
"""
fir = self.createFirFilter(fs, x)
# Decimated sampling frequency [Hz] # Decimated sampling frequency [Hz]
fd = fs / np.prod(self.decimation(x)) fd = fs / np.prod(self.firdecimation(x))
return frsp(fd, freq, fir) return firFreqResponse(fd, freq, fir)
def nominal_txt_tox(self, nom_txt): def nominal_txt_tox(self, nom_txt: str):
""" """Returns the x-value corresponding to a certain nominal txt: '1k' ->
Returns the x-value corresponding to a certain nominal txt: '1k' -> 0 0.
Args:
nom_txt: Text-representation of midband frequency
""" """
for x in self.xs: for x in self.xs:
if self.nominal_txt(x) == nom_txt: if self.nominal_txt(x) == nom_txt:
@ -96,12 +163,10 @@ class FilterBankDesigner:
class OctaveBankDesigner(FilterBankDesigner): class OctaveBankDesigner(FilterBankDesigner):
""" """Octave band filter designer."""
Octave band filter designer
"""
def __init__(self): def __init__(self, fs):
pass super().__init__(fs)
@property @property
def b(self): def b(self):
@ -110,10 +175,58 @@ class OctaveBankDesigner(FilterBankDesigner):
@property @property
def xs(self): def xs(self):
"""All possible band designators for an octave band filter."""
return list(range(-6, 5)) return list(range(-6, 5))
def band_limits(self, x, filter_class=0):
"""Returns the octave band filter limits for filter designator x.
Args:
x: Filter offset power from the reference frequency of 1000 Hz.
filter_class: Either 0 or 1, defines the tolerances on the frequency
response
Returns:
freq, llim, ulim: Tuple of Numpy arrays containing the frequencies of
the corner points of the filter frequency response limits, lower limits
in *deciBell*, upper limits in *deciBell*, respectively.
"""
b = 1
# Exact midband frequency
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_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
mininf = -1e300
if filter_class == 1:
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_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]
if filter_class == 0:
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)
return freqs, lower_limits, upper_limits
def nominal_txt(self, x): def nominal_txt(self, x):
# Text corresponding to the nominal frequency """Returns textual repressentation of corresponding to the nominal
frequency."""
nominals = {4: '16k', nominals = {4: '16k',
3: '8k', 3: '8k',
2: '4k', 2: '4k',
@ -125,12 +238,13 @@ class OctaveBankDesigner(FilterBankDesigner):
-4: '63', -4: '63',
-5: '31.5', -5: '31.5',
-6: '16'} -6: '16'}
assert len(nominals) == len(self.xs)
return nominals[x] return nominals[x]
def fac_l(self, x): def firFac_l(self, x):
""" """Factor with which to multiply the cut-on frequency of the FIR
Factor with which to multiply the cut-on frequency of the FIR filter filter."""
""" assert int(self.fs) == 48000, 'Fir coefs are only valid for 48kHz fs'
if x == 4: if x == 4:
return .995 return .995
elif x in (3, 1): elif x in (3, 1):
@ -140,10 +254,10 @@ class OctaveBankDesigner(FilterBankDesigner):
else: else:
return .96 return .96
def fac_u(self, x): def firFac_u(self, x):
""" """Factor with which to multiply the cut-off frequency of the FIR
Factor with which to multiply the cut-off frequency of the FIR filter filter."""
""" assert int(self.fs) == 48000, 'Fir coefs are only valid for 48kHz fs'
if x == 4: if x == 4:
return 1.004 return 1.004
elif x in (3, 1): elif x in (3, 1):
@ -153,10 +267,9 @@ class OctaveBankDesigner(FilterBankDesigner):
else: else:
return 1.02 return 1.02
def decimation(self, x): def firDecimation(self, x):
""" """Required firdecimation for each filter."""
Required decimation for each filter assert int(self.fs) == 48000, 'Fir coefs are only valid for 48kHz fs'
"""
if x > 1: if x > 1:
return [1] return [1]
elif x > -2: elif x > -2:
@ -167,24 +280,52 @@ class OctaveBankDesigner(FilterBankDesigner):
return [4, 4, 4] return [4, 4, 4]
elif x == -6: elif x == -6:
return [4, 4, 4, 4] return [4, 4, 4, 4]
assert False, 'Overlooked decimation' assert False, 'Overlooked firdecimation'
def sosFac_l(self, x):
"""Left side percentage of change in cut-on frequency for designing the
filter, for OCTAVE band filter.
Args:
x: Filter band designator
"""
# Idea: correct for frequency warping:
if int(self.fs) in [48000, 96000]:
return 1.0
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.
Args:
x: Filter band designator
"""
if int(self.fs) in [48000, 96000]:
return 1.0
else:
raise ValueError('Unimplemented sampling frequency for SOS'
'filter design')
class ThirdOctaveBankDesigner(FilterBankDesigner): class ThirdOctaveBankDesigner(FilterBankDesigner):
def __init__(self): def __init__(self, fs):
super().__init__(fs)
self.xs = list(range(-16, 14)) self.xs = list(range(-16, 14))
# Text corresponding to the nominal frequency # Text corresponding to the nominal frequency
self._nominal_txt = ['25', '31.5', '40', self._nominal_txt = ['25', '31.5', '40',
'50', '63', '80', '50', '63', '80',
'100', '125', '160', '100', '125', '160',
'200', '250', '315', '200', '250', '315',
'400', '500', '630', '400', '500', '630',
'800', '1k', '1.25k', '800', '1k', '1.25k',
'1.6k', '2k', '2.5k', '1.6k', '2k', '2.5k',
'3.15k', '4k', '5k', '3.15k', '4k', '5k',
'6.3k', '8k', '10k', '6.3k', '8k', '10k',
'12.5k', '16k', '20k'] '12.5k', '16k', '20k']
assert len(self.xs) == len(self._nominal_txt) assert len(self.xs) == len(self._nominal_txt)
@ -199,8 +340,61 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
index = x - self.xs[0] index = x - self.xs[0]
return self._nominal_txt[index] return self._nominal_txt[index]
@staticmethod def band_limits(self, x, filter_class=0):
def decimation(x): """Returns the third octave band filter limits for filter designator x.
Args:
x: Filter offset power from the reference frequency of 1000 Hz.
filter_class: Either 0 or 1, defines the tolerances on the frequency
response
Returns:
freq, llim, ulim: Tuple of Numpy arrays containing the frequencies of
the corner points of the filter frequency response limits, lower limits
in *deciBell*, upper limits in *deciBell*, respectively.
"""
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_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
mininf = -1e300
if filter_class == 1:
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]
else:
raise ValueError('Filter class should either be 0 or 1')
upper_limits_neg = upper_limits_pos[:]
upper_limits_neg.reverse()
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]
elif filter_class == 0:
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)
return freqs, lower_limits, upper_limits
def firDecimation(self, x):
assert int(self.fs) == 48000, 'Fir coefs are only valid for 48kHz fs'
if x > 5: if x > 5:
return [1] return [1]
elif x > -1: elif x > -1:
@ -211,10 +405,9 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
return [4, 4, 4] return [4, 4, 4]
elif x > -17: elif x > -17:
return [4, 4, 4, 4] return [4, 4, 4, 4]
assert False, 'Bug: overlooked decimation' assert False, 'Bug: overlooked firdecimation'
@staticmethod def firFac_l(self, x):
def fac_l(x):
if x in (-13, -7, -1, 5, 11, 12, 13): if x in (-13, -7, -1, 5, 11, 12, 13):
return .995 return .995
elif x in (-12, -6, 0, 6): elif x in (-12, -6, 0, 6):
@ -222,8 +415,7 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
else: else:
return .99 return .99
@staticmethod def firFac_u(self, x):
def fac_u(x):
if x in (-14, -13, -8, -7, -1, -2, 3, 4, 5, 10, 11, 12): if x in (-14, -13, -8, -7, -1, -2, 3, 4, 5, 10, 11, 12):
return 1.005 return 1.005
elif x in (12, 13): elif x in (12, 13):
@ -232,3 +424,22 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
return 1.015 return 1.015
else: else:
return 1.01 return 1.01
def sosFac_l(self, x):
"""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
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."""
if np.isclose(self.fs, 48000):
return 1
else:
raise ValueError('Unimplemented sampling frequency for SOS'
'filter design')

View File

@ -1,99 +0,0 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""!
Author: J.A. de Jong - ASCEE
Description: Limit lines for class 1 octave band filter limits according to
the ICS 17.140.50 standard.
"""
__all__ = ['G', 'fr', 'third_octave_band_limits', 'octave_band_limits']
import numpy as np
# Reference frequency
fr = 1000.
G = 10**(3/10)
def third_octave_band_limits(x):
"""
Returns the class 1 third octave band filter limits for filter designator
x.
Args:
x: Filter offset power from the reference frequency of 1000 Hz.
Returns:
freq, ulim, llim: Tuple of Numpy arrays containing the frequencyies,
upper and lower limits of the filter.
"""
b = 3
fm = G**(x/b)*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_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
mininf = -1e300
upper_limits_pos = [.3]*5 + [-2, -17.5, -42, -61, -70, -70]
upper_limits_neg = upper_limits_pos[:]
upper_limits_neg.reverse()
upper_limits = np.array(upper_limits_neg[:-1] + upper_limits_pos)
lower_limits_pos = [-.3, -.4, -.6, -1.3, -5, -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)
return freqs, upper_limits, lower_limits
def octave_band_limits(x):
b = 1
# Exact midband frequency
fm = G**(x/b)*fr
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
mininf = -1e300
lower_limits_pos = [-0.3, -0.4, -0.6, -1.3, -5.0, -5.0] + 4*[mininf]
lower_limits_neg = lower_limits_pos[:]
lower_limits_neg.reverse()
lower_limits = np.asarray(lower_limits_neg[:-1] + lower_limits_pos)
upper_limits_pos = [0.3]*5 + [-2, -17.5, -42, -61, -70]
upper_limits_neg = upper_limits_pos[:]
upper_limits_neg.reverse()
upper_limits = np.asarray(upper_limits_neg[:-1] + upper_limits_pos)
freqs = fm*G**np.asarray(G_power_values)
return freqs, upper_limits, lower_limits
if __name__ == '__main__':
from asceefigs.plot import close, Figure
close('all')
freqs, upper_limits, lower_limits = octave_band_limits(0)
f = Figure()
f.semilogx(freqs, lower_limits)
f.semilogx(freqs, upper_limits)
f.ylim(-80, 1)

View File

@ -106,8 +106,7 @@ def C_fir_design():
def show_Afir(): def show_Afir():
from asceefigs.plot import close, Figure from asceefig.plot import Figure
close('all')
fs = 48000. fs = 48000.
freq_design = np.linspace(0, 17e3, 3000) freq_design = np.linspace(0, 17e3, 3000)
@ -133,8 +132,7 @@ def show_Afir():
def show_Cfir(): def show_Cfir():
from asceefigs.plot import close, Figure from asceefig.plot import Figure
close('all')
fs = 48000. fs = 48000.
freq_design = np.linspace(0, 17e3, 3000) freq_design = np.linspace(0, 17e3, 3000)

View File

@ -38,8 +38,6 @@ class lasp_shelve:
if lasp_shelve.refcount == 0: if lasp_shelve.refcount == 0:
lasp_shelve.shelve.close() lasp_shelve.shelve.close()
lasp_shelve.shelve = None lasp_shelve.shelve = None
# Reference sound pressure level # Reference sound pressure level
P_REF = 2e-5 P_REF = 2e-5
@ -77,7 +75,6 @@ class Window:
def getCurrent(cb): def getCurrent(cb):
return Window.types[cb.currentIndex()] return Window.types[cb.currentIndex()]
class TimeWeighting: class TimeWeighting:
none = (None, 'Raw (no time weighting)') none = (None, 'Raw (no time weighting)')
uufast = (1e-4, '0.1 ms') uufast = (1e-4, '0.1 ms')
@ -105,7 +102,6 @@ class TimeWeighting:
def getCurrent(cb): def getCurrent(cb):
return TimeWeighting.types[cb.currentIndex()] return TimeWeighting.types[cb.currentIndex()]
class FreqWeighting: class FreqWeighting:
""" """
Frequency weighting types Frequency weighting types
@ -132,7 +128,6 @@ class FreqWeighting:
def getCurrent(cb): def getCurrent(cb):
return FreqWeighting.types[cb.currentIndex()] return FreqWeighting.types[cb.currentIndex()]
def getTime(fs, N, start=0): def getTime(fs, N, start=0):
""" """
Return a time array for given number of points and sampling frequency. Return a time array for given number of points and sampling frequency.
@ -145,7 +140,6 @@ def getTime(fs, N, start=0):
assert N > 0 and fs > 0 assert N > 0 and fs > 0
return np.linspace(start, start + N/fs, N, endpoint=False) return np.linspace(start, start + N/fs, N, endpoint=False)
def getFreq(fs, nfft): def getFreq(fs, nfft):
""" """
return an array of frequencies for single-sided spectra return an array of frequencies for single-sided spectra

View File

@ -3,22 +3,68 @@
"""! """!
Author: J.A. de Jong - ASCEE Author: J.A. de Jong - ASCEE
Provides the FIR implementation of the octave filter bank Provides the implementations of (fractional) octave filter banks
""" """
__all__ = ['FirOctaveFilterBank', 'FirThirdOctaveFilterBank'] __all__ = ['FirOctaveFilterBank', 'FirThirdOctaveFilterBank',
'OverallFilterBank', 'SosOctaveFilterBank',
'SosThirdOctaveFilterBank']
from .filter.filterbank_design import OctaveBankDesigner, ThirdOctaveBankDesigner from .filter.filterbank_design import (OctaveBankDesigner,
from .wrappers import Decimator, FilterBank as pyxFilterBank ThirdOctaveBankDesigner)
from .wrappers import (Decimator, FilterBank as pyxFilterBank,
SosFilterBank as pyxSosFilterBank)
import numpy as np import numpy as np
class OverallFilterBank:
"""
Dummy type filter bank. Does nothing special, only returns output in a
sensible way
"""
def __init__(self, fs):
"""
Initialize overall filter bank
"""
self.fs = fs
self.N = 0
self.xs = [0]
def filter_(self, data):
"""
Filter input data
"""
assert data.ndim == 2
assert data.shape[1] == 1, "invalid number of channels, should be 1"
if data.shape[0] == 0:
return {}
# Output given as a dictionary with x as the key
output = {}
tstart = self.N / self.fs
Ncur = data.shape[0]
tend = tstart + Ncur / self.fs
t = np.linspace(tstart, tend, Ncur, endpoint=False)
self.N += Ncur
output['Overall'] = {'t': t, 'data': data, 'x': 0}
return output
def decimation(self, x):
return [1]
class FirFilterBank: class FirFilterBank:
""" """
Single channel octave filter bank implementation Single channel (fractional) octave filter bank implementation, based on FIR
filters and sample rate decimation.
""" """
def __init__(self, fs): def __init__(self, fs, xmin, xmax):
""" """
Initialize a OctaveFilterBank object. Initialize a OctaveFilterBank object.
@ -29,7 +75,10 @@ class FirFilterBank:
assert np.isclose(fs, 48000), "Only sampling frequency" \ assert np.isclose(fs, 48000), "Only sampling frequency" \
" available is 48 kHz" " available is 48 kHz"
maxdecimation = self.decimation(self.xs[0]) self.fs = fs
self.xs = list(range(xmin, xmax + 1))
maxdecimation = self.designer.firDecimation(self.xs[0])
self.decimators = [] self.decimators = []
for dec in maxdecimation: for dec in maxdecimation:
self.decimators.append(Decimator(1, dec)) self.decimators.append(Decimator(1, dec))
@ -43,7 +92,7 @@ class FirFilterBank:
self.filterbanks = [] self.filterbanks = []
# Sort the x values in categories according to the required decimation # Sort the x values in categories according to the required decimation
for x in self.xs: for x in self.xs:
dec = self.decimation(x) dec = self.designer.firDecimation(x)
if len(dec) == 1 and dec[0] == 1: if len(dec) == 1 and dec[0] == 1:
xs_d1.append(x) xs_d1.append(x)
elif len(dec) == 1 and dec[0] == 4: elif len(dec) == 1 and dec[0] == 4:
@ -60,24 +109,25 @@ class FirFilterBank:
xs_all = [xs_d1, xs_d4, xs_d16, xs_d64, xs_d256] xs_all = [xs_d1, xs_d4, xs_d16, xs_d64, xs_d256]
for xs in xs_all: for xs in xs_all:
nominals_txt = [] nominals_txt = []
firs = np.empty((self.L, len(xs)), order='F') if len(xs) > 0:
for i, x in enumerate(xs): firs = np.empty((self.designer.firFilterLength, len(xs)), order='F')
# These are the filters that do not require lasp_decimation for i, x in enumerate(xs):
# prior to filtering # These are the filters that do not require lasp_decimation
nominals_txt.append(self.nominal_txt(x)) # prior to filtering
firs[:, i] = self.createFirFilter(fs, x) nominals_txt.append(self.designer.nominal_txt(x))
filterbank = {'fb': pyxFilterBank(firs, 1024), firs[:, i] = self.designer.createFirFilter(fs, x)
'xs': xs, filterbank = {'fb': pyxFilterBank(firs, 1024),
'nominals': nominals_txt} 'xs': xs,
self.filterbanks.append(filterbank) 'nominals': nominals_txt}
self.filterbanks.append(filterbank)
# Sample input counter. # Sample input counter.
self.N = 0 self.N = 0
# Filter output counters
self.dec = [1, 4, 16, 64, 256] self.dec = [1, 4, 16, 64, 256]
# These intial delays are found experimentally using a toneburst # Filter output counters
# These intial delays are found 'experimentally' using a toneburst
# response. # response.
self.Nf = [915, 806, 780, 582, 338] self.Nf = [915, 806, 780, 582, 338]
@ -105,8 +155,9 @@ class FirFilterBank:
t = np.linspace(tstart, tend, Nf, endpoint=False) t = np.linspace(tstart, tend, Nf, endpoint=False)
self.Nf[dec_stage] += Nf self.Nf[dec_stage] += Nf
for i, nom in enumerate(self.filterbanks[dec_stage]['nominals']): for i, nom_txt in enumerate(self.filterbanks[dec_stage]['nominals']):
output[nom] = {'t': t, 'data': filtered[:, [i]]} x = self.designer.nominal_txt_tox(nom_txt)
output[nom_txt] = {'t': t, 'data': filtered[:, [i]], 'x': x}
return output return output
@ -122,36 +173,141 @@ class FirFilterBank:
# Output given as a dictionary with x as the key # Output given as a dictionary with x as the key
output = {} output = {}
output_unsorted = {}
self.N += data.shape[0] self.N += data.shape[0]
output = {**output, **self.filterd(0, data)} output_unsorted = {**output_unsorted, **self.filterd(0, data)}
for i in range(len(self.decimators)): for i in range(len(self.decimators)):
dec_stage = i+1 dec_stage = i+1
if data.shape[0] > 0: if data.shape[0] > 0:
# Apply a decimation stage # Apply a decimation stage
data = self.decimators[i].decimate(data) data = self.decimators[i].decimate(data)
output = {**output, **self.filterd(dec_stage, data)} output_unsorted = {**output_unsorted,
**self.filterd(dec_stage, data)}
# Create sorted output
for x in self.xs:
nom_txt = self.designer.nominal_txt(x)
output[nom_txt] = output_unsorted[nom_txt]
return output return output
def decimation(self, x):
return self.designer.firDecimation(x)
class FirOctaveFilterBank(FirFilterBank, OctaveBankDesigner):
class FirOctaveFilterBank(FirFilterBank):
""" """
Filter bank which uses FIR filtering for each octave frequency band Filter bank which uses FIR filtering for each octave frequency band
""" """
def __init__(self, fs): def __init__(self, fs, xmin, xmax):
OctaveBankDesigner.__init__(self) self.designer = OctaveBankDesigner()
FirFilterBank.__init__(self, fs) FirFilterBank.__init__(self, fs, xmin, xmax)
class FirThirdOctaveFilterBank(FirFilterBank, ThirdOctaveBankDesigner): class FirThirdOctaveFilterBank(FirFilterBank):
""" """
Filter bank which uses FIR filtering for each one-third octave frequency Filter bank which uses FIR filtering for each one-third octave frequency
band. band.
""" """
def __init__(self, fs): def __init__(self, fs, xmin, xmax):
ThirdOctaveBankDesigner.__init__(self) self.designer = ThirdOctaveBankDesigner(fs)
FirFilterBank.__init__(self, fs) FirFilterBank.__init__(self, fs, xmin, xmax)
class SosFilterBank:
def __init__(self, fs, xmin, xmax):
"""
Initialize a second order sections filterbank
Args:
fs: Sampling frequency [Hz]
xmin: Minimum value for the bands
xmax: Maximum value for the bands
"""
self.fs = fs
self.xs = list(range(xmin, xmax + 1))
nfilt = len(self.xs)
self._fb = pyxSosFilterBank(nfilt, 5)
for i, x in enumerate(self.xs):
sos = self.designer.createSOSFilter(x)
self._fb.setFilter(i, sos)
self.xmin = xmin
self.xmax = xmax
self.N = 0
def filter_(self, data):
"""
Filter input data
"""
assert data.ndim == 2
assert data.shape[1] == 1, "invalid number of channels, should be 1"
if data.shape[0] == 0:
return {}
filtered_data = self._fb.filter_(data)
# Output given as a dictionary with nom_txt as the key
output = {}
tstart = self.N / self.fs
Ncur = data.shape[0]
tend = tstart + Ncur / self.fs
t = np.linspace(tstart, tend, Ncur, endpoint=False)
self.N += Ncur
for i, x in enumerate(self.xs):
# '31.5' to '16k'
nom_txt = self.designer.nominal_txt(x)
output[nom_txt] = {'t': t, 'data': filtered_data[:, [i]], 'x': x}
return output
def decimation(self, x):
return [1]
class SosThirdOctaveFilterBank(SosFilterBank):
"""
Filter bank which uses FIR filtering for each one-third octave frequency
band.
"""
def __init__(self, fs, xmin, xmax):
"""
Initialize a second order sections filterbank.
Args:
fs: Sampling frequency [Hz]
xmin: Minimum value for the bands
xmax: Maximum value for the bands
"""
self.designer = ThirdOctaveBankDesigner(fs)
SosFilterBank.__init__(self, fs, xmin, xmax)
class SosOctaveFilterBank(SosFilterBank):
"""
Filter bank which uses FIR filtering for each one-third octave frequency
band.
"""
def __init__(self, fs, xmin, xmax):
"""
Initialize a second order sections filterbank.
Args:
fs: Sampling frequency [Hz]
xmin: Minimum value for the bands
xmax: Maximum value for the bands
"""
self.designer = OctaveBankDesigner(fs)
SosFilterBank.__init__(self, fs, xmin, xmax)