Compare commits

...

18 Commits

Author SHA1 Message Date
d67fce8d06 Bugfixes in slm 2020-01-24 20:44:12 +01:00
3dd7e23587 Added statistics to sound level meter 2020-01-24 14:08:24 +01:00
6a756bf1b1 Merged master 2020-01-23 17:23:33 +01:00
0bbfe22d83 Bugfix for lowest frequency point in converting octave band to narrow band 2020-01-23 17:23:13 +01:00
22f1cb59ff Update SLM API. Now works properly with different start and stop bands, and only overall band 2020-01-22 21:11:20 +01:00
685329d307 Almost-working SLM implementation in LASP, only statistics not yet done 2020-01-21 21:10:38 +01:00
98421893f0 Bugfix in tracer. Found using clang_check 2020-01-21 21:07:07 +01:00
bd88882f25 Negative time constant for single pole lowpass filter means no time-weighting at all. 2020-01-21 21:06:17 +01:00
29dc70fad3 Some tracing 2020-01-20 20:19:09 +01:00
b4360951a8 Updated wrappers. SLM almost fully functioning. Only does not produce statistics (max, min, eq) yet. 2020-01-20 20:18:37 +01:00
388fa0ed0a Bugfix in SLM for resetting offset with multi-channel filter banks 2020-01-20 20:15:52 +01:00
b4ed3138d4 Reversed poles of weighting filters to make filter stable 2020-01-20 20:13:24 +01:00
b3197a3dee Midband frequencies can now also be called with list of xs 2020-01-20 19:47:09 +01:00
d91007dbaa Updated sound pressure weighting filters API 2020-01-20 19:46:41 +01:00
1258ff4919 Merged master 2020-01-20 12:11:50 +01:00
438c657f65 First work on SLM. Seems to be working properly without pre-filtering and bandpass bank 2020-01-20 12:10:24 +01:00
b5088bef14 Merge branch 'splweighting_sos' into slm_c_impl 2020-01-18 13:45:28 +01:00
31c6fcdc22 Initial commit. Deleted single pole lowpass header file 2020-01-17 20:57:23 +01:00
15 changed files with 714 additions and 245 deletions

BIN
img/subsampling_slm.pdf Normal file

Binary file not shown.

View File

@ -1,5 +1,6 @@
if(!LASP_DEBUG)
SET_SOURCE_FILES_PROPERTIES(lasp_sosfilterbank.c PROPERTIES COMPILE_FLAGS -O3)
SET_SOURCE_FILES_PROPERTIES(lasp_slm.c PROPERTIES COMPILE_FLAGS -O3)
endif(!LASP_DEBUG)
add_library(lasp_lib
@ -19,7 +20,7 @@ add_library(lasp_lib
lasp_firfilterbank.c
lasp_sosfilterbank.c
lasp_decimation.c
lasp_sp_lowpass.c
lasp_slm.c
)

View File

@ -15,6 +15,11 @@
* else when required. For example for debugging purposes.
*/
static inline void* a_malloc(size_t nbytes) {
#if LASP_DEBUG
if(nbytes == 0) {
FATAL("Tried to allocate 0 bytes");
}
#endif
void* ptr = malloc(nbytes);
if(!ptr) {
FATAL("Memory allocation failed. Exiting");

View File

@ -445,10 +445,22 @@ static inline void dmat_copy(dmat* to,const dmat* from) {
dbgassert(to->n_cols==from->n_cols,SIZEINEQUAL);
for(us col=0;col<to->n_cols;col++) {
d_copy(getdmatval(to,0,col),
getdmatval(from,0,col),
to->n_rows,1,1);
getdmatval(from,0,col),
to->n_rows,1,1);
}
}
/**
* Allocate a new array, with size based on other.
*
* @param[in] from: Array to copy
*/
static inline dmat dmat_alloc_from_dmat(const dmat* from) {
assertvalidptr(from);
dmat thecopy = dmat_alloc(from->n_rows, from->n_cols);
return thecopy;
}
/**
* Copy contents of one matrix to another. Sizes should be equal
*
@ -547,17 +559,36 @@ static inline void cmat_conj_inplace(cmat* x) {
}
}
/**
* Computes the maximum value for each row, returns a vector with maximum
* values for each column in the matrix.
*
* @param x
*/
static inline vd dmat_max(const dmat x) {
vd max_vals = vd_alloc(x.n_cols);
d max_val = -d_inf;
for(us j=0; j< x.n_cols; j++) {
for(us i=0; i< x.n_rows; i++) {
max_val = *getdmatval(&x, i, j) < max_val? max_val : *getdmatval(&x, i,j);
}
*getvdval(&max_vals, j) = max_val;
}
return max_vals;
}
#ifdef LASP_DEBUG
void print_cmat(const cmat* m);
void print_dmat(const dmat* m);
#define print_vc print_cmat
#define print_vd print_dmat
#define print_vc(x) assert_vx(x) print_cmat(x)
#define print_vd(x) assert_vx(x) print_dmat(x)
#else
#define print_cmat(m)
#define print_vc(m)
#define print_dmat(m)
#define print_vc(m)
#define print_vd(m)
#endif
#endif // LASP_MAT_H

View File

@ -11,6 +11,7 @@
#include "lasp_assert.h"
#include "lasp_tracer.h"
#include "lasp_types.h"
#include <float.h>
#include <math.h>
#if LASP_USE_BLAS == 1
@ -33,6 +34,8 @@
#define d_sin sin
#define d_cos cos
#define d_pow pow
#define d_log10 log10
#define d_epsilon (DBL_EPSILON)
#else // LASP_DOUBLE_PRECISION not defined
#define c_conj conjf
@ -47,9 +50,13 @@
#define d_sin sinf
#define d_cos cosf
#define d_pow powf
#define d_log10 log10f
#define d_epsilon (FLT_EPSILON)
#endif // LASP_DOUBLE_PRECISION
/// Positive infinite
#define d_inf (INFINITY)
#ifdef M_PI
static const d number_pi = M_PI;
#else

293
lasp/c/lasp_slm.c Normal file
View File

@ -0,0 +1,293 @@
#define TRACERPLUS (-5)
#include "lasp_slm.h"
#include "lasp_assert.h"
#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 Pm; /// Storage for the computing the mean of the square of the signal.
vd Pmax; /// Storage for maximum computed signal power so far.
vd Ppeak; /// Storage for computing peak powers so far.
us N; /// Counter for the number of time samples counted that came in
} 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);
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'.
us ds_fac;
if (tau > 0) {
// A reasonable 'framerate' for the sound level meter, based on the
// filtering time constant.
d fs_slm = 10 / tau;
dVARTRACE(15, fs_slm);
if(fs_slm < 30) {
fs_slm = 30;
}
ds_fac = (us)(fs / fs_slm);
if (ds_fac == 0) {
// If we get 0, it should be 1
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;
}
/// Initialize statistics gatherers
slm->Ppeak = vd_alloc(filterbank_size);
slm->Pmax = vd_alloc(filterbank_size);
slm->Pm = vd_alloc(filterbank_size);
slm->N = 0;
vd_set(&(slm->Ppeak), 0);
vd_set(&(slm->Pmax), 0);
vd_set(&(slm->Pm), 0);
feTRACE(15);
return slm;
}
dmat Slm_run(Slm *slm, vd *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;
/// 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);
iVARTRACE(15, downsampling_fac);
us cur_offset = slm->cur_offset;
/// Compute the number of samples output
us nsamples_output = samples_bandpassed;
if (downsampling_fac > 1) {
nsamples_output = (samples_bandpassed - cur_offset) / downsampling_fac;
if(nsamples_output > samples_bandpassed) {
// This means overflow of unsigned number calculations
nsamples_output = 0;
}
while(nsamples_output * downsampling_fac + cur_offset < samples_bandpassed) {
nsamples_output++;
}
}
iVARTRACE(15, nsamples_output);
iVARTRACE(15, cur_offset);
dmat levels = dmat_alloc(nsamples_output, filterbank_size);
us N, ch;
for (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;
*getvdval(&(slm->Ppeak), ch) = d_max(*getvdval(&(slm->Ppeak), ch), *tmp);
}
// Now that all data for the channel is squared, we can run it through
// the low-pass filter
cur_offset = slm->cur_offset;
/// 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");
/// Output resulting levels at a lower interval
us i = 0;
N = slm->N;
d *Pm = getvdval(&(slm->Pm), ch);
while (cur_offset < samples_bandpassed) {
iVARTRACE(10, i);
iVARTRACE(10, cur_offset);
/// Filtered power.
const d P = *getvdval(&power_filtered, cur_offset);
dVARTRACE(15, P);
/// Compute maximum, compare to current maximum
*getvdval(&(slm->Pmax), ch) = d_max(*getvdval(&(slm->Pmax), ch), P);
/// Update mean power
d Nd = (d) N;
*Pm = (*Pm*Nd + P ) / (Nd+1);
N++;
dVARTRACE(15, *Pm);
/// Compute level
d level = 10 * d_log10((P + d_epsilon ) / 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);
}
/// Update sample counter
dbgassert(ch >0, "BUG");
slm->N = N;
slm->cur_offset = cur_offset - samples_bandpassed;
vd_free(&prefiltered);
dmat_free(&bandpassed);
feTRACE(15);
return levels;
}
static inline vd levels_from_power(const vd* power,const d ref_level){
fsTRACE(15);
vd levels = dmat_alloc_from_dmat(power);
for(us i=0; i< levels.n_rows; i++) {
*getvdval(&levels, i) = 10 * d_log10(
(*getvdval(power, i) + d_epsilon) / ref_level / ref_level);
}
feTRACE(15);
return levels;
}
vd Slm_Lpeak(Slm* slm) {
fsTRACE(15);
assertvalidptr(slm);
vd Lpeak = levels_from_power(&(slm->Ppeak), slm->ref_level);
feTRACE(15);
return Lpeak;
}
vd Slm_Lmax(Slm* slm) {
fsTRACE(15);
assertvalidptr(slm);
vd Lmax = levels_from_power(&(slm->Pmax), slm->ref_level);
feTRACE(15);
return Lmax;
}
vd Slm_Leq(Slm* slm) {
fsTRACE(15);
assertvalidptr(slm);
print_vd(&(slm->Pm));
vd Leq = levels_from_power(&(slm->Pm), slm->ref_level);
feTRACE(15);
return Leq;
}
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]);
}
a_free(slm->splowpass);
}
vd_free(&(slm->Ppeak));
vd_free(&(slm->Pmax));
vd_free(&(slm->Pm));
a_free(slm);
feTRACE(15);
}

96
lasp/c/lasp_slm.h Normal file
View File

@ -0,0 +1,96 @@
// lasp_slm.h
//
// Author: J.A. de Jong - ASCEE
//
// Description: Multi-purpose implementation of a real time Sound Level Meter,
// can be used for full-signal filtering, (fractional) octave band filtering,
// etc.
// //////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_SLM_H
#define LASP_SLM_H
#include "lasp_types.h"
#include "lasp_mat.h"
#include "lasp_sosfilterbank.h"
typedef struct Slm Slm;
#define TAU_FAST (1.0/8.0)
#define TAU_SLOW (1.0)
#define TAU_IMPULSE (35e-3)
/**
* Initializes a Sound level meter. NOTE: Sound level meter takes over
* ownership of pointers to the filterbanks for prefiltering and band-pass
* filtering! After passing them to the constructor of the Slm, they should not
* be touched anymore.
*
* @param[in] weighting: Sosfiterbank handle, used to pre-filter data. This is
* in most cases an A-weighting filter, or C-weighting. If NULL, no
* pre-filtering is done. That can be the case in situations of Z-weighting.
* @param[in] fb: Sosfiterbank handle, bandpass filters.
* @param[in] fs: Sampling frequency in [Hz], used for computing the
* downsampling factor and size of output arrays.
* @param[in] tau: Time constant of single pole low pass filter for squared data.
* Three pre-defined values can be used: TAU_FAST, for fast filtering,
* TAU_SLOW, for slow filtering, TAU_IMPULSE for impulse filtering
* @param[in] ref_level: Reference level when computing dB's. I.e. P_REF_AIR for
* sound pressure levels in air
* @param[out] downsampling_fac: Here, the used downsampling factor is stored
* which is used on returning data. If the value is for example 10, the
* 'sampling' frequency of output data from `Slm_run` is 4800 is fs is set to
* 48000. This downsampling factor is a function of the used time weighting.
*
* @return Slm: Handle of the sound level meter, NULL on error.
*/
Slm* Slm_create(Sosfilterbank* weighting,Sosfilterbank* bandpass,
const d fs,
const d tau,
const d ref_level,
us* downsampling_fac);
/**
* Run the sound level meter on a piece of time data.
*
* @param[in] slm: Slm handle
* @param[in] input_data: Vector of input data samples.
*
* @return Output result of Sound Level values in [dB], for each bank in the filter
* bank.
*/
dmat Slm_run(Slm* slm,
vd* input_data);
/**
* Cleans up an existing Slm
*
* @param f slm handle
*/
void Slm_free(Slm* f);
/**
* Returns the (raw) peak level computed so far.
*
* @param[in] slm: Slm handle
* @return Vector of peak level for each channel in the filterbank.
*/
vd Slm_Lpeak(Slm* slm);
/**
* Returns the equivalent level computed so far.
*
* @param[in] slm: Slm handle
* @return Vector of equivalent levels for each channel in the filterbank.
*/
vd Slm_Leq(Slm* slm);
/**
* Returns the maximum level computed so far.
*
* @param[in] slm: Slm handle
* @return Vector of maximum levels for each channel in the filterbank.
*/
vd Slm_Lmax(Slm* slm);
#endif // LASP_SLM_H
//////////////////////////////////////////////////////////////////////

View File

@ -73,6 +73,8 @@ void Sosfilterbank_setFilter(Sosfilterbank* fb,const us filter_no,
fsTRACE(15);
assertvalidptr(fb);
assert_vx(&filter_coefs);
iVARTRACE(15, filter_coefs.n_rows);
iVARTRACE(15, filter_no);
dbgassert(filter_no < fb->filterbank_size, "Illegal filter number");
dbgassert(filter_coefs.n_rows == fb->nsections * 6,
"Illegal filter coefficient length");

View File

@ -1,91 +0,0 @@
// lasp_sp_lowpass.c
//
// Author: J.A. de Jong -ASCEE
//
// Description:
// Single pole lowpass IIR filter implementation.
//////////////////////////////////////////////////////////////////////
#include "lasp_sp_lowpass.h"
#include "lasp_alloc.h"
typedef struct SPLowpass_s {
d a;
d b;
d xlast,ylast;
} SPLowpass;
SPLowpass* SPLowpass_create(d fs,d tau) {
fsTRACE(15);
d tau_i = 1/tau;
d T = 1/fs;
if(fs <= 0) {
WARN("Invalid sampling frequency");
return NULL;
}
else if(T >= tau ) {
WARN("Invalid tau, should be (much) larger than sampling"
" time (1/fs)");
return NULL;
}
SPLowpass* lp = a_malloc(sizeof(SPLowpass));
lp->xlast = 0;
lp->ylast = 0;
lp->a = tau_i/(2*fs+tau_i);
lp->b = (2*fs-tau_i)/(2*fs+tau_i);
feTRACE(15);
return lp;
}
vd SPLowpass_filter(SPLowpass* lp,
const vd* input) {
fsTRACE(15);
dbgassert(lp && input,NULLPTRDEREF);
assert_vx(input);
us input_size = input->n_rows;
if(input_size == 0) {
return vd_alloc(0);
}
vd output = vd_alloc(input_size);
const d xlast = lp->xlast;
const d ylast = lp->ylast;
const d a = lp->a;
const d b = lp->b;
*getvdval(&output,0) = a*(xlast + *getvdval(input,0)) +
b*ylast;
for(us i=1;i<input_size;i++) {
*getvdval(&output,i) = a*(*getvdval(input,i-1) +
*getvdval(input,i)) +
b*(*getvdval(&output,i-1));
}
lp->xlast = *getvdval(input ,input_size-1);
lp->ylast = *getvdval(&output,input_size-1);
feTRACE(15);
return output;
}
void SPLowpass_free(SPLowpass* lp) {
fsTRACE(15);
dbgassert(lp,NULLPTRDEREF);
a_free(lp);
feTRACE(15);
}
//////////////////////////////////////////////////////////////////////

View File

@ -1,51 +0,0 @@
// lasp_sp_lowpass.h
//
// Author: J.A. de Jong - ASCEE
//
// Description: Single-pole low pass IIR filter. Used for example for
// fast and fast and slow averaging of the square of the squared
// A-weighted acoustic pressure. This class uses the bilinear
// transform to derive the filter coefficients.
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_SP_LOWPASS_H
#define LASP_SP_LOWPASS_H
#include "lasp_types.h"
#include "lasp_mat.h"
typedef struct SPLowpass_s SPLowpass;
/**
* Create a single pole lowpass IIR filter.
*
* @param fs Sampling frequency [Hz]
* @param tau Time constant of the lowpass filter. Should be >0,
* otherwise an error occurs
*
* @return Pointer to dynamically allocated lowpass filter. NULL on
* error.
*/
SPLowpass* SPLowpass_create(d fs,d tau);
/**
* Use the filter to filter input data
*
* @param lp Lowpass filter handlen
* @param input Vector of input samples.
*
* @return Output data. Length is equal to input at all cases.
*/
vd SPLowpass_filter(SPLowpass* lp,
const vd* input);
/**
* Free a single pole lowpass filter
*
* @param lp
*/
void SPLowpass_free(SPLowpass* lp);
#endif // LASP_SP_LOWPASS_H
//////////////////////////////////////////////////////////////////////

View File

@ -8,6 +8,7 @@
#pragma once
#ifndef ASCEE_TRACER_H
#define ASCEE_TRACER_H
#include "lasp_types.h"
#include <string.h>
#include <stdio.h>
#include <stdlib.h>

View File

@ -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.
@ -256,7 +259,7 @@ class FilterBankDesigner:
for x in self.xs:
xl = x
fl = self.fl(x)
if self.fl(x+1) > freq[0]:
if fl >= freq[0]:
break
# Find upper frequency xu
@ -569,3 +572,4 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
else:
raise ValueError('Unimplemented sampling frequency for SOS'
'filter design')

View File

@ -76,15 +76,15 @@ class Window:
return Window.types[cb.currentIndex()]
class TimeWeighting:
none = (None, 'Raw (no time weighting)')
none = (-1, 'Raw (no time weighting)')
uufast = (1e-4, '0.1 ms')
ufast = (30e-3, '30 ms')
ufast = (35e-3, 'Impulse (35 ms)')
fast = (0.125, 'Fast (0.125 s)')
slow = (1.0, 'Slow (1.0 s)')
tens = (10, '10 s')
infinite = (np.Inf, 'Infinite (Leq)')
types = (none, uufast, ufast, fast, slow, tens, infinite)
default = 3
tens = (10., '10 s')
types = (none, uufast, ufast, fast, slow, tens)
default = fast
default_index = 3
@staticmethod
def fillComboBox(cb):
@ -97,8 +97,9 @@ class TimeWeighting:
cb.clear()
for tw in TimeWeighting.types:
cb.addItem(tw[1], tw)
cb.setCurrentIndex(TimeWeighting.default)
cb.setCurrentIndex(TimeWeighting.default_index)
@staticmethod
def getCurrent(cb):
return TimeWeighting.types[cb.currentIndex()]
@ -110,7 +111,8 @@ class FreqWeighting:
C = ('C', 'C-weighting')
Z = ('Z', 'Z-weighting')
types = (A, C, Z)
default = 0
default = A
default_index = 0
@staticmethod
def fillComboBox(cb):
@ -123,8 +125,9 @@ class FreqWeighting:
cb.clear()
for fw in FreqWeighting.types:
cb.addItem(fw[1], fw)
cb.setCurrentIndex(FreqWeighting.default)
cb.setCurrentIndex(FreqWeighting.default_index)
@staticmethod
def getCurrent(cb):
return FreqWeighting.types[cb.currentIndex()]

View File

@ -4,9 +4,10 @@
Sound level meter implementation
@author: J.A. de Jong - ASCEE
"""
from .wrappers import SPLowpass
from .wrappers import Slm as pyxSlm
import numpy as np
from .lasp_common import (TimeWeighting, P_REF)
from .lasp_common import (TimeWeighting, P_REF, FreqWeighting)
from .filter import SPLFilterDesigner
__all__ = ['SLM', 'Dummy']
@ -24,78 +25,171 @@ class SLM:
"""
Multi-channel Sound Level Meter. Input data: time data with a certain
sampling frequency. Output: time-weighted (fast/slow) sound pressure
levels in dB(A/C/Z).
levels in dB(A/C/Z). Possibly in octave bands.
"""
def __init__(self, fs, tw=TimeWeighting.default):
def __init__(self,
fs,
fbdesigner,
tw=TimeWeighting.fast,
fw=FreqWeighting.A,
xmin = None,
xmax = None,
include_overall=True):
"""
Initialize a sound level meter object.
Number of channels comes from the weighcal object.
Args:
fs: Sampling frequency of input data [Hz]
fbdesigner: FilterBankDesigner to use for creating the
(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.
if tw[0] is not TimeWeighting.none[0]:
# Initialize the single-pole low-pass filter for given time-
# weighting value.
self._lp = SPLowpass(fs, tw[0])
"""
self.fbdesigner = fbdesigner
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
spld = SPLFilterDesigner(fs)
if fw == FreqWeighting.A:
prefilter = spld.A_Sos_design().flatten()
elif fw == FreqWeighting.C:
prefilter = spld.C_Sos_design().flatten()
elif fw == FreqWeighting.Z:
prefilter = None
else:
self._lp = Dummy()
self._Lmax = 0.
raise ValueError(f'Not implemented prefilter {fw}')
# Storage for computing the equivalent level
self._sq = 0. # Square of the level data, storage
self._N = 0
self._Leq = 0.
# 'Probe' size of filter coefficients
self.nom_txt = []
@property
def Leq(self):
"""
Returns the equivalent level of the recorded signal so far
"""
return self._Leq
if fbdesigner is not None:
assert fbdesigner.fs == fs
sos0 = fbdesigner.createSOSFilter(self.xs[0]).flatten()
self.nom_txt.append(fbdesigner.nominal_txt(self.xs[0]))
sos = np.empty((nfilters, sos0.size), dtype=float, order='C')
sos[0, :] = sos0
@property
def Lmax(self):
"""
Returns the currently maximum recorded level
"""
return self._Lmax
for i, x in enumerate(self.xs[1:]):
sos[i+1, :] = fbdesigner.createSOSFilter(x).flatten()
self.nom_txt.append(fbdesigner.nominal_txt(x))
def addData(self, data):
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,
fs, tw[0], P_REF)
dsfac = self.slm.downsampling_fac
if dsfac > 0:
# Not unfiltered data
self.fs_slm = fs / self.slm.downsampling_fac
else:
self.fs_slm = fs
# Initialize counter to 0
self.N = 0
def run(self, data):
"""
Add new fresh timedata to the Sound Level Meter
Args:
data:
returns:
level values as a function of time
data: one-dimensional input data
"""
assert data.ndim == 2
assert data.shape[1] == 1
assert data.shape[1] == 1, "invalid number of channels, should be 1"
P_REFsq = P_REF**2
if data.shape[0] == 0:
return {}
# Squared
sq = data**2
levels = self.slm.run(data)
# Update equivalent level
N1 = sq.shape[0]
self._sq = (np.sum(sq) + self._sq*self._N)/(N1+self._N)
self._N += N1
self._Leq = 10*np.log10(self._sq/P_REFsq)
tstart = self.N / self.fs_slm
Ncur = levels.shape[0]
tend = tstart + Ncur / self.fs_slm
# Time-weight the signal
tw = self._lp.filter_(sq)
t = np.linspace(tstart, tend, Ncur, endpoint=False)
self.N += Ncur
Level = 10*np.log10(tw/P_REFsq)
output = {}
# Update maximum level
curmax = np.max(Level)
if curmax > self._Lmax:
self._Lmax = curmax
for i, x in enumerate(self.xs):
# '31.5' to '16k'
output[self.nom_txt[i]] = {'t': t,
'data': levels[:, [i]],
'x': x}
if self.include_overall and self.fbdesigner is not None:
output['overall'] = {'t': t, 'data': levels[:, i+1], 'x': 0}
return Level
return output
def return_as_dict(self, dat):
"""
Helper function used to
"""
output = {}
for i, x in enumerate(self.xs):
# '31.5' to '16k'
output[self.nom_txt[i]] = { 'data': dat[i],
'x': x}
if self.include_overall and self.fbdesigner is not None:
output['overall'] = {'data': dat[i+1], 'x': 0}
return output
def Leq(self):
"""
Returns the computed equivalent levels for each filter channel
"""
return self.return_as_dict(self.slm.Leq())
def Lmax(self):
"""
Returns the computed max levels for each filter channel
"""
return self.return_as_dict(self.slm.Lmax())
def Lpeak(self):
"""
Returns the computed peak levels for each filter channel
"""
return self.return_as_dict(self.slm.Lpeak())
def Leq_array(self):
return self.slm.Leq()
def Lmax_array(self):
return self.slm.Lmax()
def Lpeak_array(self):
return self.slm.Lpeak()

View File

@ -455,6 +455,117 @@ cdef extern from "lasp_decimation.h":
dmat Decimator_decimate(c_Decimator* dec,const dmat* samples) nogil
void Decimator_free(c_Decimator* dec) nogil
cdef extern from "lasp_slm.h" nogil:
ctypedef struct c_Slm "Slm"
d TAU_FAST, TAU_SLOW, TAU_IMPULSE
c_Slm* Slm_create(c_Sosfilterbank* prefilter,
c_Sosfilterbank* bandpass,
d fs, d tau, d ref_level,
us* downsampling_fac)
dmat Slm_run(c_Slm* slm,vd* input_data)
void Slm_free(c_Slm* slm)
vd Slm_Leq(c_Slm*)
vd Slm_Lmax(c_Slm*)
vd Slm_Lpeak(c_Slm*)
tau_fast = TAU_FAST
tau_slow = TAU_SLOW
tau_impulse = TAU_IMPULSE
cdef class Slm:
cdef:
c_Slm* c_slm
public us downsampling_fac
def __cinit__(self, d[::1] sos_prefilter,
d[:, ::1] sos_bandpass,
d fs, d tau, d ref_level):
cdef:
us prefilter_nsections
us bandpass_nsections
us bandpass_nchannels
c_Sosfilterbank* prefilter = NULL
c_Sosfilterbank* bandpass = NULL
vd coefs_vd
d[:] coefs
if sos_prefilter is not None:
assert sos_prefilter.size % 6 == 0
prefilter_nsections = sos_prefilter.size // 6
prefilter = Sosfilterbank_create(1,prefilter_nsections)
coefs = sos_prefilter
coefs_vd = dmat_foreign_data(prefilter_nsections*6,1,
&coefs[0],False)
Sosfilterbank_setFilter(prefilter, 0, coefs_vd)
if prefilter is NULL:
raise RuntimeError('Error creating pre-filter')
if sos_bandpass is not None:
assert sos_bandpass.shape[1] % 6 == 0
bandpass_nsections = sos_bandpass.shape[1] // 6
bandpass_nchannels = sos_bandpass.shape[0]
bandpass = Sosfilterbank_create(bandpass_nchannels,
bandpass_nsections)
if bandpass == NULL:
if prefilter:
Sosfilterbank_free(prefilter)
raise RuntimeError('Error creating bandpass filter')
for i in range(bandpass_nchannels):
coefs = sos_bandpass[i, :]
coefs_vd = dmat_foreign_data(bandpass_nsections*6,1,
&coefs[0],False)
Sosfilterbank_setFilter(bandpass, i, coefs_vd)
self.c_slm = Slm_create(prefilter, bandpass,
fs, tau, ref_level,
&self.downsampling_fac)
if self.c_slm is NULL:
Sosfilterbank_free(prefilter)
Sosfilterbank_free(bandpass)
raise RuntimeError('Error creating sound level meter')
def run(self, d[:, ::1] data):
assert data.shape[1] == 1
cdef vd data_vd = dmat_foreign_data(data.shape[0], 1,
&data[0,0], False)
cdef dmat res
with nogil:
res = Slm_run(self.c_slm, &data_vd)
result = dmat_to_ndarray(&res,True)
return result
def Leq(self):
cdef vd res
res = Slm_Leq(self.c_slm)
# True below, means transfer ownership of allocated data to Numpy
return dmat_to_ndarray(&res,True)
def Lmax(self):
cdef vd res
res = Slm_Lmax(self.c_slm)
# True below, means transfer ownership of allocated data to Numpy
return dmat_to_ndarray(&res,True)
def Lpeak(self):
cdef vd res
res = Slm_Lpeak(self.c_slm)
# True below, means transfer ownership of allocated data to Numpy
return dmat_to_ndarray(&res,True)
def __dealloc__(self):
if self.c_slm:
Slm_free(self.c_slm)
cdef class Decimator:
cdef:
c_Decimator* dec
@ -485,43 +596,6 @@ cdef class Decimator:
if self.dec != NULL:
Decimator_free(self.dec)
cdef extern from "lasp_sp_lowpass.h":
ctypedef struct c_SPLowpass "SPLowpass"
c_SPLowpass* SPLowpass_create(d fs,d tau)
vd SPLowpass_filter(c_SPLowpass* lp,
const vd* _input)
void SPLowpass_free(c_SPLowpass* lp)
cdef class SPLowpass:
cdef:
c_SPLowpass* lp
def __cinit__(self,d fs,d tau):
self.lp = SPLowpass_create(fs,tau)
if not self.lp:
raise RuntimeError('Error creating lowpass filter')
def __dealloc__(self):
if self.lp:
SPLowpass_free(self.lp)
def filter_(self,d[::1,:] input_):
assert input_.shape[1] == 1
if input_.shape[0] == 0:
return np.array([],dtype=NUMPY_FLOAT_TYPE)
cdef vd input_vd = dmat_foreign_data(input_.shape[0],1,
&input_[0,0],False)
cdef dmat output = SPLowpass_filter(self.lp,&input_vd)
# # Steal the pointer from output
result = dmat_to_ndarray(&output,True)
dmat_free(&output)
vd_free(&input_vd)
return result
cdef extern from "lasp_siggen.h":
ctypedef struct c_Siggen "Siggen"