diff --git a/img/subsampling_slm.pdf b/img/subsampling_slm.pdf new file mode 100644 index 0000000..e0cdeac Binary files /dev/null and b/img/subsampling_slm.pdf differ diff --git a/lasp/c/CMakeLists.txt b/lasp/c/CMakeLists.txt index fed41fc..fa3e715 100644 --- a/lasp/c/CMakeLists.txt +++ b/lasp/c/CMakeLists.txt @@ -19,7 +19,7 @@ add_library(lasp_lib lasp_firfilterbank.c lasp_sosfilterbank.c lasp_decimation.c - lasp_sp_lowpass.c + lasp_slm.c ) diff --git a/lasp/c/lasp_math_raw.h b/lasp/c/lasp_math_raw.h index 06cf60d..66a2901 100644 --- a/lasp/c/lasp_math_raw.h +++ b/lasp/c/lasp_math_raw.h @@ -33,6 +33,7 @@ #define d_sin sin #define d_cos cos #define d_pow pow +#define d_log10 log10 #else // LASP_DOUBLE_PRECISION not defined #define c_conj conjf @@ -47,7 +48,7 @@ #define d_sin sinf #define d_cos cosf #define d_pow powf - +#define d_log10 log10f #endif // LASP_DOUBLE_PRECISION #ifdef M_PI diff --git a/lasp/c/lasp_slm.c b/lasp/c/lasp_slm.c new file mode 100644 index 0000000..4fd3795 --- /dev/null +++ b/lasp/c/lasp_slm.c @@ -0,0 +1,185 @@ +#define TRACERPLUS (10) +#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 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); + + Slm *slm = NULL; + if (tau <= 0) { + WARN("Invalid time constant"); + return NULL; + } else 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'. + 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; + + /// Create the single pole lowpass + us filterbank_size; + if (bandpass) { + filterbank_size = Sosfilterbank_getFilterbankSize(bandpass); + } else { + filterbank_size = 1; + } + + 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); + + 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); + us cur_offset = slm->cur_offset; + + 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 = dmat_alloc(nsamples_output, filterbank_size); + + for (us ch = 0; ch < bandpassed.n_cols; ch++) { + iVARTRACE(15, ch); + /// Inplace squaring of the signal + for (us sample = 0; sample < bandpassed.n_rows; sample++) { + tmp = getdmatval(&bandpassed, sample, ch); + *tmp = *tmp * *tmp; + } + + // Now that all data for the channel is squared, we can run it through + // the low-pass filter + vd chan = dmat_column(&bandpassed, ch); + + /// 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"); + + /// Output resulting levels at a lower interval + us i = 0; + d level; + while (cur_offset < samples_bandpassed) { + iVARTRACE(10, i); + iVARTRACE(10, cur_offset); + 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 == nsamples_output, "BUG"); + slm->cur_offset = cur_offset - samples_bandpassed; + + vd_free(&chan); + vd_free(&power_filtered); + } + + 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); + } + assertvalidptr(slm->splowpass); + + us filterbank_size; + if (slm->bandpass) { + filterbank_size = Sosfilterbank_getFilterbankSize(slm->bandpass); + Sosfilterbank_free(slm->bandpass); + } else { + filterbank_size = 1; + } + for (us ch = 0; ch < filterbank_size; ch++) { + Sosfilterbank_free(slm->splowpass[ch]); + } + a_free(slm->splowpass); + + a_free(slm); + + feTRACE(15); +} diff --git a/lasp/c/lasp_slm.h b/lasp/c/lasp_slm.h index 7df37ba..67970ca 100644 --- a/lasp/c/lasp_slm.h +++ b/lasp/c/lasp_slm.h @@ -2,60 +2,71 @@ // // Author: J.A. de Jong - ASCEE // -// Description: Implementation of a real time Sound Level Meter, based on -// given slm. +// 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 +#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. + * 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 fb: Sosfiterbank handle, used to pre-filter data - * @param T: Single pole low pass filter time constant to use. + * @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* fb); +Slm* Slm_create(Sosfilterbank* weighting,Sosfilterbank* bandpass, + const d fs, + const d tau, + const d ref_level, + us* downsampling_fac); /** - * Initialize the filter coeficients in the slm + * Run the sound level meter on a piece of time data. * - * @param fb: slm handle - * @param filter_no: Filter number in the bank - * @param coefss: Array of filter coefficients. Should have a length of - * nsections x 6, for each of the sections, it contains (b0, b1, b2, a0, - * a1, a2), where a are the numerator coefficients and b are the denominator - * coefficients. + * @param[in] slm: Slm handle + * @param[in] input_data: Vector of input data samples. * -*/ -void slm_setFilter(slm* fb,const us filter_no, - const vd coefs); - -/** - * Filters x using h, returns y - * - * @param x Input time sequence block. Should have at least one sample. - - * @return Filtered output in an allocated array. The number of - * columns in this array equals the number of filters in the - * slm. The number of output samples is equal to the number of - * input samples in x. + * @return Output result of Sound Level values in [dB], for each bank in the filter + * bank. */ -dmat slm_filter(slm* fb, - const vd* x); - +dmat Slm_run(Slm* slm, + vd* input_data); /** - * Cleans up an existing filter bank. + * Cleans up an existing Slm * * @param f slm handle */ -void slm_free(slm* f); +void Slm_free(Slm* f); -#endif // LASP_slm_H +#endif // LASP_SLM_H ////////////////////////////////////////////////////////////////////// diff --git a/lasp/lasp_slm.py b/lasp/lasp_slm.py index 930e265..8ac595c 100644 --- a/lasp/lasp_slm.py +++ b/lasp/lasp_slm.py @@ -4,7 +4,7 @@ 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) @@ -24,11 +24,11 @@ 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, tw=TimeWeighting.default, ): """ Initialize a sound level meter object. Number of channels comes from the weighcal object. diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx index e05a2ef..52d2103 100644 --- a/lasp/wrappers.pyx +++ b/lasp/wrappers.pyx @@ -455,6 +455,76 @@ 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": + 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) + + +tau_fast = TAU_FAST +tau_slow = TAU_SLOW +tau_impulse = TAU_IMPULSE + +cdef class Slm: + cdef: + c_Slm* c_slm + 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 = sos_prefilter.shape[1] // 6 + us bandpass_nsections = sos_bandpass.shape[1] // 6 + us bandpass_nchannels = sos_bandpass.shape[0] + c_Sosfilterbank* prefilter = NULL + c_Sosfilterbank* bandpass = NULL + + if sos_prefilter is not None: + if sos_prefilter.shape[0] != 1: + raise ValueError('Prefilter should have only one channel') + prefilter = Sosfilterbank_create(1,prefilter_nsections) + if prefilter is NULL: + raise RuntimeError('Error creating pre-filter') + + if sos_bandpass is not None: + bandpass = Sosfilterbank_create(bandpass_nchannels, + bandpass_nsections) + if bandpass == NULL: + if prefilter: + Sosfilterbank_free(prefilter) + raise RuntimeError('Error creating bandpass filter') + + 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): + cdef vd data_vd = dmat_foreign_data(data.shape[0], 1, + &data[0], False) + cdef dmat res = Slm_run(self.c_slm, &data_vd) + result = dmat_to_ndarray(&res,True) + return result + + + def __dealloc__(self): + if self.c_slm: + Slm_free(self.c_slm) + + + + cdef class Decimator: cdef: c_Decimator* dec