Update SLM API. Now works properly with different start and stop bands, and only overall band
This commit is contained in:
parent
685329d307
commit
22f1cb59ff
@ -4,206 +4,207 @@
|
||||
#include "lasp_tracer.h"
|
||||
|
||||
typedef struct Slm {
|
||||
Sosfilterbank *prefilter; /// Pre-filter, A, or C. If NULL, not used.
|
||||
Sosfilterbank *bandpass; /// Filterbank. If NULL, not used
|
||||
Sosfilterbank **splowpass; /// Used for time-weighting of the squared signal
|
||||
d ref_level; /// Reference value for computing decibels
|
||||
us downsampling_fac; /// Every x'th sample is returned.
|
||||
us cur_offset; /// Storage for offset point in input arrays
|
||||
vd Leq; /// Storage for the computed equivalent levels so far.
|
||||
Sosfilterbank *prefilter; /// Pre-filter, A, or C. If NULL, not used.
|
||||
Sosfilterbank *bandpass; /// Filterbank. If NULL, not used
|
||||
Sosfilterbank **splowpass; /// Used for time-weighting of the squared signal
|
||||
d ref_level; /// Reference value for computing decibels
|
||||
us downsampling_fac; /// Every x'th sample is returned.
|
||||
us cur_offset; /// Storage for offset point in input arrays
|
||||
vd Leq; /// Storage for the computed equivalent levels so far.
|
||||
|
||||
} Slm;
|
||||
|
||||
Slm *Slm_create(Sosfilterbank *prefilter, Sosfilterbank *bandpass, const d fs,
|
||||
const d tau, const d ref_level, us *downsampling_fac) {
|
||||
fsTRACE(15);
|
||||
assertvalidptr(downsampling_fac);
|
||||
const d tau, const d ref_level, us *downsampling_fac) {
|
||||
fsTRACE(15);
|
||||
assertvalidptr(downsampling_fac);
|
||||
|
||||
Slm *slm = NULL;
|
||||
if (ref_level <= 0) {
|
||||
WARN("Invalid reference level");
|
||||
return NULL;
|
||||
} else if (fs <= 0) {
|
||||
WARN("Invalid sampling frequency");
|
||||
return NULL;
|
||||
}
|
||||
|
||||
slm = (Slm *)a_malloc(sizeof(Slm));
|
||||
slm->ref_level = ref_level;
|
||||
slm->prefilter = prefilter;
|
||||
slm->bandpass = bandpass;
|
||||
|
||||
/// Compute the downsampling factor. This one is chosen based on the
|
||||
/// lowpass filter. Which has a -3 dB point of f = 1/(tau*2*pi). See LASP
|
||||
/// documentation for the computation of its minus 20 dB point. We set the
|
||||
/// reduction in its 'sampling frequency' such that its noise is at a level
|
||||
/// of 20 dB less than its 'signal'.
|
||||
if (tau > 0) {
|
||||
const d fs_slm = 1 / (2 * number_pi * tau) * (1 - 0.01) / 0.01;
|
||||
slm->downsampling_fac = (us)(fs / fs_slm);
|
||||
slm->cur_offset = 0;
|
||||
*downsampling_fac = slm->downsampling_fac;
|
||||
} else {
|
||||
*downsampling_fac = 1;
|
||||
slm->downsampling_fac = 1;
|
||||
}
|
||||
|
||||
/// Create the single pole lowpass
|
||||
us filterbank_size;
|
||||
if (bandpass) {
|
||||
filterbank_size = Sosfilterbank_getFilterbankSize(bandpass);
|
||||
} else {
|
||||
filterbank_size = 1;
|
||||
}
|
||||
|
||||
if (tau > 0) {
|
||||
|
||||
vd lowpass_sos = vd_alloc(6);
|
||||
d b0 = 1.0 / (1 + 2 * tau * fs);
|
||||
*getvdval(&lowpass_sos, 0) = b0;
|
||||
*getvdval(&lowpass_sos, 1) = b0;
|
||||
*getvdval(&lowpass_sos, 2) = 0;
|
||||
*getvdval(&lowpass_sos, 3) = 1;
|
||||
*getvdval(&lowpass_sos, 4) = (1 - 2 * tau * fs) * b0;
|
||||
*getvdval(&lowpass_sos, 5) = 0;
|
||||
|
||||
slm->splowpass = a_malloc(filterbank_size * sizeof(Sosfilterbank *));
|
||||
for (us ch = 0; ch < filterbank_size; ch++) {
|
||||
/// Allocate a filterbank with one channel and one section.
|
||||
slm->splowpass[ch] = Sosfilterbank_create(1, 1);
|
||||
Sosfilterbank_setFilter(slm->splowpass[ch], 0, lowpass_sos);
|
||||
Slm *slm = NULL;
|
||||
if (ref_level <= 0) {
|
||||
WARN("Invalid reference level");
|
||||
return NULL;
|
||||
} else if (fs <= 0) {
|
||||
WARN("Invalid sampling frequency");
|
||||
return NULL;
|
||||
}
|
||||
vd_free(&lowpass_sos);
|
||||
} else {
|
||||
/// No low-pass filtering. Tau set to zero
|
||||
slm->splowpass = NULL;
|
||||
}
|
||||
|
||||
feTRACE(15);
|
||||
return slm;
|
||||
slm = (Slm *)a_malloc(sizeof(Slm));
|
||||
slm->ref_level = ref_level;
|
||||
slm->prefilter = prefilter;
|
||||
slm->bandpass = bandpass;
|
||||
|
||||
/// Compute the downsampling factor. This one is chosen based on the
|
||||
/// lowpass filter. Which has a -3 dB point of f = 1/(tau*2*pi). See LASP
|
||||
/// documentation for the computation of its minus 20 dB point. We set the
|
||||
/// reduction in its 'sampling frequency' such that its noise is at a level
|
||||
/// of 20 dB less than its 'signal'.
|
||||
us ds_fac;
|
||||
if (tau > 0) {
|
||||
const d fs_slm = 1 / (2 * number_pi * tau) * (1 - 0.01) / 0.01;
|
||||
dVARTRACE(15, fs_slm);
|
||||
ds_fac = (us) (fs / fs_slm);
|
||||
// If we get 0, it should be 1
|
||||
if(ds_fac == 0) ds_fac++;
|
||||
} else {
|
||||
ds_fac = 1;
|
||||
}
|
||||
slm->downsampling_fac = ds_fac;
|
||||
*downsampling_fac = ds_fac;
|
||||
slm->cur_offset = 0;
|
||||
|
||||
/// Create the single pole lowpass
|
||||
us filterbank_size;
|
||||
if (bandpass) {
|
||||
filterbank_size = Sosfilterbank_getFilterbankSize(bandpass);
|
||||
} else {
|
||||
filterbank_size = 1;
|
||||
}
|
||||
|
||||
if (tau > 0) {
|
||||
|
||||
vd lowpass_sos = vd_alloc(6);
|
||||
d b0 = 1.0 / (1 + 2 * tau * fs);
|
||||
*getvdval(&lowpass_sos, 0) = b0;
|
||||
*getvdval(&lowpass_sos, 1) = b0;
|
||||
*getvdval(&lowpass_sos, 2) = 0;
|
||||
*getvdval(&lowpass_sos, 3) = 1;
|
||||
*getvdval(&lowpass_sos, 4) = (1 - 2 * tau * fs) * b0;
|
||||
*getvdval(&lowpass_sos, 5) = 0;
|
||||
|
||||
slm->splowpass = a_malloc(filterbank_size * sizeof(Sosfilterbank *));
|
||||
for (us ch = 0; ch < filterbank_size; ch++) {
|
||||
/// Allocate a filterbank with one channel and one section.
|
||||
slm->splowpass[ch] = Sosfilterbank_create(1, 1);
|
||||
Sosfilterbank_setFilter(slm->splowpass[ch], 0, lowpass_sos);
|
||||
}
|
||||
vd_free(&lowpass_sos);
|
||||
} else {
|
||||
/// No low-pass filtering. Tau set to zero
|
||||
slm->splowpass = NULL;
|
||||
}
|
||||
|
||||
feTRACE(15);
|
||||
return slm;
|
||||
}
|
||||
|
||||
dmat Slm_run(Slm *slm, vd *input_data) {
|
||||
fsTRACE(15);
|
||||
assertvalidptr(slm);
|
||||
assert_vx(input_data);
|
||||
fsTRACE(15);
|
||||
assertvalidptr(slm);
|
||||
assert_vx(input_data);
|
||||
|
||||
/// First step: run the input data through the pre-filter
|
||||
vd prefiltered;
|
||||
if (slm->prefilter)
|
||||
prefiltered = Sosfilterbank_filter(slm->prefilter, input_data);
|
||||
else {
|
||||
prefiltered = dmat_foreign(input_data);
|
||||
}
|
||||
dmat bandpassed;
|
||||
if (slm->bandpass) {
|
||||
bandpassed = Sosfilterbank_filter(slm->bandpass, &prefiltered);
|
||||
} else {
|
||||
bandpassed = dmat_foreign(&prefiltered);
|
||||
}
|
||||
us filterbank_size = bandpassed.n_cols;
|
||||
/// First step: run the input data through the pre-filter
|
||||
vd prefiltered;
|
||||
if (slm->prefilter)
|
||||
prefiltered = Sosfilterbank_filter(slm->prefilter, input_data);
|
||||
else {
|
||||
prefiltered = dmat_foreign(input_data);
|
||||
}
|
||||
dmat bandpassed;
|
||||
if (slm->bandpass) {
|
||||
bandpassed = Sosfilterbank_filter(slm->bandpass, &prefiltered);
|
||||
} else {
|
||||
bandpassed = dmat_foreign(&prefiltered);
|
||||
}
|
||||
us filterbank_size = bandpassed.n_cols;
|
||||
|
||||
/// Next step: square all values. We do this in-place. Then we filter for
|
||||
/// each channel.
|
||||
d ref_level = slm->ref_level;
|
||||
d *tmp;
|
||||
/// Next step: square all values. We do this in-place. Then we filter for
|
||||
/// each channel.
|
||||
d ref_level = slm->ref_level;
|
||||
d *tmp;
|
||||
|
||||
/// Pre-calculate the size of the output data
|
||||
us downsampling_fac = slm->downsampling_fac;
|
||||
us samples_bandpassed = bandpassed.n_rows;
|
||||
iVARTRACE(15, samples_bandpassed);
|
||||
us cur_offset = slm->cur_offset;
|
||||
/// Pre-calculate the size of the output data
|
||||
us downsampling_fac = slm->downsampling_fac;
|
||||
us samples_bandpassed = bandpassed.n_rows;
|
||||
iVARTRACE(15, samples_bandpassed);
|
||||
iVARTRACE(15, downsampling_fac);
|
||||
us cur_offset = slm->cur_offset;
|
||||
|
||||
/// Compute the number of samples output
|
||||
int nsamples_output = (samples_bandpassed - cur_offset) / downsampling_fac;
|
||||
while (nsamples_output * downsampling_fac + cur_offset < samples_bandpassed)
|
||||
nsamples_output++;
|
||||
if (nsamples_output < 0)
|
||||
nsamples_output = 0;
|
||||
|
||||
iVARTRACE(15, nsamples_output);
|
||||
iVARTRACE(15, cur_offset);
|
||||
|
||||
dmat levels;
|
||||
if (slm->splowpass) {
|
||||
levels = dmat_alloc(nsamples_output, filterbank_size);
|
||||
} else {
|
||||
levels = dmat_alloc(samples_bandpassed, filterbank_size);
|
||||
}
|
||||
|
||||
for (us ch = 0; ch < bandpassed.n_cols; ch++) {
|
||||
iVARTRACE(15, ch);
|
||||
vd chan = dmat_column(&bandpassed, ch);
|
||||
/// Inplace squaring of the signal
|
||||
for (us sample = 0; sample < bandpassed.n_rows; sample++) {
|
||||
tmp = getdmatval(&bandpassed, sample, ch);
|
||||
*tmp = *tmp * *tmp;
|
||||
/// Compute the number of samples output
|
||||
int nsamples_output = samples_bandpassed;
|
||||
if(downsampling_fac > 1) {
|
||||
nsamples_output = (samples_bandpassed - cur_offset) / downsampling_fac;
|
||||
while (nsamples_output * downsampling_fac + cur_offset < samples_bandpassed)
|
||||
nsamples_output++;
|
||||
if (nsamples_output < 0)
|
||||
nsamples_output = 0;
|
||||
}
|
||||
|
||||
// Now that all data for the channel is squared, we can run it through
|
||||
// the low-pass filter
|
||||
if (slm->splowpass) {
|
||||
cur_offset = slm->cur_offset;
|
||||
iVARTRACE(15, nsamples_output);
|
||||
iVARTRACE(15, cur_offset);
|
||||
dmat levels = dmat_alloc(nsamples_output, filterbank_size);
|
||||
|
||||
/// Apply single-pole lowpass filter for current filterbank channel
|
||||
vd power_filtered = Sosfilterbank_filter(slm->splowpass[ch], &chan);
|
||||
dbgassert(chan.n_rows == power_filtered.n_rows, "BUG");
|
||||
for (us ch = 0; ch < bandpassed.n_cols; ch++) {
|
||||
iVARTRACE(15, ch);
|
||||
vd chan = dmat_column(&bandpassed, ch);
|
||||
/// Inplace squaring of the signal
|
||||
for (us sample = 0; sample < bandpassed.n_rows; sample++) {
|
||||
tmp = getdmatval(&bandpassed, sample, ch);
|
||||
*tmp = *tmp * *tmp;
|
||||
}
|
||||
|
||||
/// Output resulting levels at a lower interval
|
||||
us i = 0;
|
||||
while (cur_offset < samples_bandpassed) {
|
||||
iVARTRACE(10, i);
|
||||
iVARTRACE(10, cur_offset);
|
||||
/// Compute level
|
||||
d level = 10 * d_log10(*getvdval(&power_filtered, cur_offset) /
|
||||
ref_level / ref_level);
|
||||
// Now that all data for the channel is squared, we can run it through
|
||||
// the low-pass filter
|
||||
cur_offset = slm->cur_offset;
|
||||
|
||||
*getdmatval(&levels, i++, ch) = level;
|
||||
cur_offset = cur_offset + downsampling_fac;
|
||||
}
|
||||
iVARTRACE(15, cur_offset);
|
||||
iVARTRACE(15, i);
|
||||
dbgassert(i == (int) nsamples_output, "BUG");
|
||||
/// Apply single-pole lowpass filter for current filterbank channel
|
||||
TRACE(15, "Start filtering");
|
||||
vd power_filtered;
|
||||
if(slm->splowpass) {
|
||||
power_filtered = Sosfilterbank_filter(slm->splowpass[ch], &chan);
|
||||
} else {
|
||||
power_filtered = dmat_foreign(&chan);
|
||||
}
|
||||
TRACE(15, "Filtering done");
|
||||
dbgassert(chan.n_rows == power_filtered.n_rows, "BUG");
|
||||
|
||||
vd_free(&chan);
|
||||
vd_free(&power_filtered);
|
||||
/// Output resulting levels at a lower interval
|
||||
us i = 0;
|
||||
while (cur_offset < samples_bandpassed) {
|
||||
iVARTRACE(10, i);
|
||||
iVARTRACE(10, cur_offset);
|
||||
/// Compute level
|
||||
d level = 10 * d_log10(*getvdval(&power_filtered, cur_offset) /
|
||||
ref_level / ref_level);
|
||||
|
||||
*getdmatval(&levels, i++, ch) = level;
|
||||
cur_offset = cur_offset + downsampling_fac;
|
||||
}
|
||||
iVARTRACE(15, cur_offset);
|
||||
iVARTRACE(15, i);
|
||||
dbgassert(i == (int) nsamples_output, "BUG");
|
||||
|
||||
vd_free(&power_filtered);
|
||||
vd_free(&chan);
|
||||
}
|
||||
}
|
||||
slm->cur_offset = cur_offset - samples_bandpassed;
|
||||
slm->cur_offset = cur_offset - samples_bandpassed;
|
||||
|
||||
if (!slm->splowpass) {
|
||||
/// Raw copy of to levels. Happens only when the low-pass filter does not
|
||||
/// have to come into action.
|
||||
dmat_copy(&levels, &bandpassed);
|
||||
}
|
||||
|
||||
vd_free(&prefiltered);
|
||||
dmat_free(&bandpassed);
|
||||
feTRACE(15);
|
||||
return levels;
|
||||
vd_free(&prefiltered);
|
||||
dmat_free(&bandpassed);
|
||||
feTRACE(15);
|
||||
return levels;
|
||||
}
|
||||
|
||||
void Slm_free(Slm *slm) {
|
||||
fsTRACE(15);
|
||||
assertvalidptr(slm);
|
||||
if (slm->prefilter) {
|
||||
Sosfilterbank_free(slm->prefilter);
|
||||
}
|
||||
|
||||
us filterbank_size;
|
||||
if (slm->bandpass) {
|
||||
filterbank_size = Sosfilterbank_getFilterbankSize(slm->bandpass);
|
||||
Sosfilterbank_free(slm->bandpass);
|
||||
} else {
|
||||
filterbank_size = 1;
|
||||
}
|
||||
if (slm->splowpass) {
|
||||
for (us ch = 0; ch < filterbank_size; ch++) {
|
||||
Sosfilterbank_free(slm->splowpass[ch]);
|
||||
fsTRACE(15);
|
||||
assertvalidptr(slm);
|
||||
if (slm->prefilter) {
|
||||
Sosfilterbank_free(slm->prefilter);
|
||||
}
|
||||
a_free(slm->splowpass);
|
||||
}
|
||||
a_free(slm);
|
||||
|
||||
feTRACE(15);
|
||||
us filterbank_size;
|
||||
if (slm->bandpass) {
|
||||
filterbank_size = Sosfilterbank_getFilterbankSize(slm->bandpass);
|
||||
Sosfilterbank_free(slm->bandpass);
|
||||
} else {
|
||||
filterbank_size = 1;
|
||||
}
|
||||
if (slm->splowpass) {
|
||||
for (us ch = 0; ch < filterbank_size; ch++) {
|
||||
Sosfilterbank_free(slm->splowpass[ch]);
|
||||
}
|
||||
a_free(slm->splowpass);
|
||||
}
|
||||
a_free(slm);
|
||||
|
||||
feTRACE(15);
|
||||
}
|
||||
|
||||
@ -67,6 +67,9 @@ class FilterBankDesigner:
|
||||
return bool(np.all(llim_full <= h_dB) and
|
||||
np.all(ulim_full >= h_dB))
|
||||
|
||||
def band_limits(self, x, filter_class):
|
||||
raise NotImplementedError()
|
||||
|
||||
def fm(self, x):
|
||||
"""Returns the exact midband frequency of the bandpass filter.
|
||||
|
||||
@ -569,3 +572,4 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
|
||||
else:
|
||||
raise ValueError('Unimplemented sampling frequency for SOS'
|
||||
'filter design')
|
||||
|
||||
|
||||
@ -30,32 +30,42 @@ class SLM:
|
||||
"""
|
||||
|
||||
def __init__(self,
|
||||
fs,
|
||||
fbdesigner,
|
||||
tw=TimeWeighting.fast,
|
||||
fw=FreqWeighting.A,
|
||||
xmin = None,
|
||||
xmax = None,
|
||||
include_overall=True):
|
||||
"""
|
||||
Initialize a sound level meter object.
|
||||
|
||||
Args:
|
||||
fs: Sampling frequency of input data [Hz]
|
||||
fbdesigner: FilterBankDesigner to use for creating the
|
||||
(fractional) octave bank filters
|
||||
(fractional) octave bank filters. Set this one to None to only do
|
||||
overalls
|
||||
fs: Sampling frequency [Hz]
|
||||
tw: Time Weighting to apply
|
||||
fw: Frequency weighting to apply
|
||||
xmin: Filter designator of lowest band
|
||||
xmax: Filter designator of highest band
|
||||
include_overall: If true, a non-functioning filter is added which
|
||||
is used to compute the overall level.
|
||||
|
||||
"""
|
||||
|
||||
self.fbdesigner = fbdesigner
|
||||
self.xs = fbdesigner.xs[:]
|
||||
if xmin is None:
|
||||
xmin = fbdesigner.xs[0]
|
||||
if xmax is None:
|
||||
xmax = fbdesigner.xs[-1]
|
||||
self.xs = list(range(xmin, xmax + 1))
|
||||
|
||||
nfilters = len(self.xs)
|
||||
if include_overall: nfilters +=1
|
||||
self.include_overall = include_overall
|
||||
|
||||
fs = fbdesigner.fs
|
||||
spld = SPLFilterDesigner(fs)
|
||||
if fw == FreqWeighting.A:
|
||||
prefilter = spld.A_Sos_design().flatten()
|
||||
@ -64,24 +74,33 @@ class SLM:
|
||||
elif fw == FreqWeighting.Z:
|
||||
prefilter = None
|
||||
else:
|
||||
raise ValueError('Not implemented prefilter')
|
||||
raise ValueError(f'Not implemented prefilter {fw}')
|
||||
|
||||
# 'Probe' size of filter coefficients
|
||||
self.nom_txt = []
|
||||
|
||||
sos0 = fbdesigner.createSOSFilter(self.xs[0]).flatten()
|
||||
sos = np.empty((nfilters, sos0.size), dtype=float, order='C')
|
||||
sos[0, :] = sos0
|
||||
if fbdesigner is not None:
|
||||
assert fbdesigner.fs == fs
|
||||
sos0 = fbdesigner.createSOSFilter(self.xs[0]).flatten()
|
||||
sos = np.empty((nfilters, sos0.size), dtype=float, order='C')
|
||||
sos[0, :] = sos0
|
||||
|
||||
for i, x in enumerate(self.xs[1:]):
|
||||
sos[i, :] = fbdesigner.createSOSFilter(x).flatten()
|
||||
self.nom_txt.append(fbdesigner.nominal_txt(x))
|
||||
for i, x in enumerate(self.xs[1:]):
|
||||
sos[i, :] = fbdesigner.createSOSFilter(x).flatten()
|
||||
self.nom_txt.append(fbdesigner.nominal_txt(x))
|
||||
|
||||
if include_overall:
|
||||
# Create a unit impulse response filter, every third index equals
|
||||
# 1, so b0 = 1 and a0 is 1 (by definition)
|
||||
sos[-1,:] = 0
|
||||
sos[-1,::3] = 1
|
||||
if include_overall:
|
||||
# Create a unit impulse response filter, every third index equals
|
||||
# 1, so b0 = 1 and a0 is 1 (by definition)
|
||||
sos[-1,:] = 0
|
||||
sos[-1,::3] = 1
|
||||
self.nom_txt.append('overall')
|
||||
else:
|
||||
# No filterbank, means we do only compute the overall values. This
|
||||
# means that in case of include_overall, it creates two overall
|
||||
# channels. That would be confusing, so we do not allow it.
|
||||
assert include_overall == False
|
||||
sos = None
|
||||
self.nom_txt.append('overall')
|
||||
|
||||
self.slm = pyxSlm(prefilter, sos,
|
||||
@ -142,7 +161,7 @@ class SLM:
|
||||
output[self.nom_txt[i]] = {'t': t,
|
||||
'data': levels[:, i],
|
||||
'x': x}
|
||||
if self.include_overall:
|
||||
if self.include_overall and self.fbdesigner is not None:
|
||||
output['overall'] = {'t': t, 'data': levels[:, i+1], 'x': 0}
|
||||
return output
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user