Compare commits
18 Commits
183404338e
...
d67fce8d06
Author | SHA1 | Date | |
---|---|---|---|
d67fce8d06 | |||
3dd7e23587 | |||
6a756bf1b1 | |||
0bbfe22d83 | |||
22f1cb59ff | |||
685329d307 | |||
98421893f0 | |||
bd88882f25 | |||
29dc70fad3 | |||
b4360951a8 | |||
388fa0ed0a | |||
b4ed3138d4 | |||
b3197a3dee | |||
d91007dbaa | |||
1258ff4919 | |||
438c657f65 | |||
b5088bef14 | |||
31c6fcdc22 |
BIN
img/subsampling_slm.pdf
Normal file
BIN
img/subsampling_slm.pdf
Normal file
Binary file not shown.
@ -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
|
||||
)
|
||||
|
||||
|
||||
|
@ -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");
|
||||
|
@ -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
|
||||
|
@ -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
293
lasp/c/lasp_slm.c
Normal 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
96
lasp/c/lasp_slm.h
Normal 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
|
||||
//////////////////////////////////////////////////////////////////////
|
@ -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");
|
||||
|
@ -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);
|
||||
}
|
||||
|
||||
//////////////////////////////////////////////////////////////////////
|
@ -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
|
||||
//////////////////////////////////////////////////////////////////////
|
@ -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>
|
||||
|
@ -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')
|
||||
|
||||
|
@ -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()]
|
||||
|
||||
|
192
lasp/lasp_slm.py
192
lasp/lasp_slm.py
@ -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()
|
||||
|
@ -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"
|
||||
|
Loading…
x
Reference in New Issue
Block a user