ASCEE/lasp
ASCEE/lasp
5
1
Fork 0
Code Issues 5 Pull Requests Actions Releases 14 Activity

Typo fix. Added 44.1 kHz as possible sampling frequency for computing filters.

Browse Source
This commit is contained in:
Anne de Jong 2022-10-06 21:34:33 +02:00
parent f7a49dc4ff
commit 90a5d87419
2 changed files with 90 additions and 74 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/lasp/dsp/lasp_slm.cpp
View File

@ -67,7 +67,7 @@ std::vector<unique_ptr<Filter>> 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.

144
src/lasp/filter/filterbank_design.py
View File

@ -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,31 +358,32 @@ 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',
nominals = {
4: '16k',
3: '8k',
2: '4k',
1: '2k',
@ -385,7 +393,8 @@ class OctaveBankDesigner(FilterBankDesigner):
-3: '125',
-4: '63',
-5: '31.5',
-6: '16'}
-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."""