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

Added suggested SLM down sampling factor, improved quite a lot of documentation. Measurement object can now work with old, as well as new measurement files.

Browse Source
This commit is contained in:
Anne de Jong 2022-10-01 19:59:35 +02:00
parent 59085b091f
commit bb26fc6bcc
15 changed files with 214 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/lasp/__init__.py
View File

@ -11,7 +11,7 @@ __version__ = lasp_cpp.__version__
# from .lasp_imptube import * # TwoMicImpedanceTube
from .lasp_measurement import * # Measurement, scaleBlockSens
from .lasp_octavefilter import *
# from .lasp_slm import * # SLM, Dummy
from .lasp_slm import * # SLM, Dummy
from .lasp_record import * # RecordStatus, Recording
from .lasp_daqconfigs import *
# from .lasp_siggen import * # SignalType, NoiseType, SiggenMessage, SiggenData, Siggen

13
src/lasp/device/lasp_daqconfig.h
View File

@ -200,6 +200,19 @@ public:
* output.
*/
double digitalHighPassCutOn = -1;
/**
* @brief Compare two channels to eachother. They are equal when the name,
* sensitivity and quantity are the same.
*
* @param other The DaqChannel to compare with
*
* @return true if equal
*/
bool operator==(const DaqChannel& other) const {
return other.name == name && other.sensitivity == sensitivity &&
other.qty == qty;
}
};
/**

4
src/lasp/dsp/lasp_avpowerspectra.h
View File

@ -7,6 +7,10 @@
#include <optional>
/** \defgroup dsp Digital Signal Processing utilities
* These are classes and functions used for processing raw signal data, to
* obtain statistics.
*
*
* @{
*/

5
src/lasp/dsp/lasp_biquadbank.h
View File

@ -2,7 +2,7 @@
#include "lasp_filter.h"
/**
* \ingroup dsp
* \addtogroup dsp
* @{
*/
@ -53,7 +53,8 @@ public:
/**
* @brief Initialize biquadbank.
*
* @param filters Filters for each filter in the bank
* @param filters Filters for each filter in the bank. First axis isis the
* coefficient index, second axis is the filter index.
* @param gains Gain values. Given as pointer, if not given (nulltpr), gains
* are initialized with unity gain.
*/

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

@ -9,29 +9,29 @@
using std::cerr;
using std::endl;
using std::runtime_error;
using rte = std::runtime_error;
using std::unique_ptr;
SLM::SLM(const d fs, const d Lref, const us downsampling_fac, const d tau,
std::unique_ptr<Filter> pre_filter,
std::vector<std::unique_ptr<Filter>> bandpass)
: _pre_filter(std::move(pre_filter)), _bandpass(std::move(bandpass)),
_alpha(exp(-1 / (fs * tau))),
_sp_storage(_bandpass.size(), arma::fill::zeros), // Storage for
// components of
// single pole low pass
// filter
Lrefsq(Lref*Lref), // Reference level
downsampling_fac(downsampling_fac),
std::unique_ptr<Filter> pre_filter,
std::vector<std::unique_ptr<Filter>> bandpass)
: _pre_filter(std::move(pre_filter)), _bandpass(std::move(bandpass)),
_alpha(exp(-1 / (fs * tau))),
_sp_storage(_bandpass.size(), arma::fill::zeros), // Storage for
// components of
// single pole low pass
// filter
Lrefsq(Lref*Lref), // Reference level
downsampling_fac(downsampling_fac),
// Initalize mean square
Pm(_bandpass.size(), arma::fill::zeros),
// Initalize mean square
Pm(_bandpass.size(), arma::fill::zeros),
// Initalize max
Pmax(_bandpass.size(), arma::fill::zeros),
// Initalize max
Pmax(_bandpass.size(), arma::fill::zeros),
// Initalize peak
Ppeak(_bandpass.size(), arma::fill::zeros)
// Initalize peak
Ppeak(_bandpass.size(), arma::fill::zeros)
{
DEBUGTRACE_ENTER;
@ -40,13 +40,13 @@ SLM::SLM(const d fs, const d Lref, const us downsampling_fac, const d tau,
getPool();
if (Lref <= 0) {
throw runtime_error("Invalid reference level");
throw rte("Invalid reference level");
}
if (tau <= 0) {
throw runtime_error("Invalid time constant for Single pole lowpass filter");
throw rte("Invalid time constant for Single pole lowpass filter");
}
if (fs <= 0) {
throw runtime_error("Invalid sampling frequency");
throw rte("Invalid sampling frequency");
}
}
SLM::~SLM() {}
@ -66,25 +66,40 @@ std::vector<unique_ptr<Filter>> createBandPass(const dmat &coefs) {
}
return bf;
}
us SLM::suggestedDownSamplingFac(const d fs,const d tau) {
if(tau<0) throw rte("Invalid time weightin time constant");
if(fs<=0) throw rte("Invalid sampling frequency");
// A reasonable 'framerate' for the sound level meter, based on the
// filtering time constant.
if (tau > 0) {
d fs_slm = 10 / tau;
if(fs_slm < 30) {
fs_slm = 30;
}
return std::max((us) 1, static_cast<us>(fs / fs_slm));
} else {
return 1;
}
}
SLM SLM::fromBiquads(const d fs, const d Lref, const us downsampling_fac,
const d tau, const vd &pre_filter_coefs,
const dmat &bandpass_coefs) {
const d tau, const vd &pre_filter_coefs,
const dmat &bandpass_coefs) {
DEBUGTRACE_ENTER;
return SLM(fs, Lref, downsampling_fac, tau,
std::make_unique<SeriesBiquad>(pre_filter_coefs),
createBandPass(bandpass_coefs));
std::make_unique<SeriesBiquad>(pre_filter_coefs),
createBandPass(bandpass_coefs));
}
SLM SLM::fromBiquads(const d fs, const d Lref, const us downsampling_fac,
const d tau, const dmat &bandpass_coefs) {
const d tau, const dmat &bandpass_coefs) {
DEBUGTRACE_ENTER;
return SLM(fs, Lref, downsampling_fac, tau,
nullptr, // Pre-filter
createBandPass(bandpass_coefs) // Bandpass coefficients
);
nullptr, // Pre-filter
createBandPass(bandpass_coefs) // Bandpass coefficients
);
}
vd SLM::run_single(vd work,const us i) {
@ -111,7 +126,7 @@ vd SLM::run_single(vd work,const us i) {
for (us j = 0; j < work.n_rows; j++) {
// Update mean square of signal, work is here still signal power
Pm(i) = (Pm(i) * static_cast<d>(N_local) + work(j)) /
(static_cast<d>(N_local) + 1);
(static_cast<d>(N_local) + 1);
N_local++;
@ -155,8 +170,8 @@ dmat SLM::run(const vd &input_orig) {
/* DEBUGTRACE_PRINT(res.n_rows); */
/* DEBUGTRACE_PRINT(futs.size()); */
// Update the total number of samples harvested so far. NOTE: *This should be
// done AFTER the threads are done!!!*
// Update the total number of samples harvested so far. NOTE: *This should be
// done AFTER the threads are done!!!*
}
N += input.n_rows;

17
src/lasp/dsp/lasp_slm.h
View File

@ -98,7 +98,8 @@ public:
* @param downsampling_fac Every 1/downsampling_fac value is returned from
* compute()
* @param tau Time consant of level meter
* @param bandpass_coefs Biquad filter coeffiecients for bandpass filter
* @param bandpass_coefs Biquad filter coefficients for bandpass filter. First axis isis the coefficient index, second axis is the filter index.
*
* @return Sound Level Meter object
*/
@ -149,6 +150,20 @@ public:
*/
vd Lmax() const { return 10 * arma::log10(Pmax / Lrefsq); };
/**
* @brief Comput a 'suggested' downsampling factor, i.e. a lower frame rate
* at which sound level meter values are returned from the computation. This
* is possible since the signal power is low-pas filtered with a single pole
* low pass filter. It can remove computational burden, especially for
* plotting, to have a value > 10.
*
* @param fs Sampling frequency of signal [Hz]
* @param tw Time weighting of SLM low pass filter
*
* @return Suggested downsampling factor, no unit. [-]
*/
static us suggestedDownSamplingFac(const d fs,const d tw);
private:
vd run_single(vd input, const us filter_no);
};

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

@ -460,7 +460,15 @@ class OctaveBankDesigner(FilterBankDesigner):
class ThirdOctaveBankDesigner(FilterBankDesigner):
def __init__(self, fs):
def __init__(self, fs: float):
"""
Initialize ThirdOctaveBankDesigner, a filter bank designer for
one-third octave bands
Args:
fs: Sampling frequency in [Hz]
"""
super().__init__(fs)
self.xs = list(range(-16, 14))
# Text corresponding to the nominal frequency

25
src/lasp/lasp_common.py
View File

@ -84,7 +84,6 @@ class Qty:
@unique
class SIQtys(Enum):
N = Qty(name='Number',
@ -113,26 +112,20 @@ class SIQtys(Enum):
cpp_enum = DaqChannel.Qty.Voltage
)
@staticmethod
def fillComboBox(cb):
"""
Fill to a combobox
Args:
cb: QComboBox to fill
"""
cb.clear()
for ty in SIQtys:
cb.addItem(f'{ty.value.unit_name}')
cb.setCurrentIndex(0)
@staticmethod
def default():
return SIQtys.N.value
@staticmethod
def getCurrent(cb):
return list(SIQtys)[cb.currentIndex()]
def fromCppEnum(enum):
"""
Convert enumeration index from - say - a measurement file back into
physical quantity information.
"""
for qty in SIQtys:
if qty.value.cpp_enum.value == enum:
return qty.value
raise RuntimeError(f'Qty corresponding to enum {enum} not found')
@dataclass

50
src/lasp/lasp_measurement.py
View File

@ -16,9 +16,13 @@ is assumed to be acoustic data.
'sensitivity': (Optionally) the stored sensitivity of the record channels.
This can be a single value, or a list of sensitivities for
each channel. Both representations are allowed.
For measurement files of LASP < v1.0
'qtys' : (Optionally): list of quantities that is recorded for each channel',
if this array is not found. Quantities are defaulted to 'Number / Full scale'
For measurement files of LASP >= 1.0
- Datasets:
'audio': 3-dimensional array of blocks of audio data. The first axis is the
@ -46,7 +50,7 @@ from .lasp_config import LASP_NUMPY_FLOAT_TYPE
from scipy.io import wavfile
import os, time, wave, logging
from .lasp_common import SIQtys, Qty, getFreq
from .lasp_cpp import Window, DaqChannel
from .lasp_cpp import Window, DaqChannel, LASP_VERSION_MAJOR
def getSampWidth(dtype):
@ -208,6 +212,13 @@ class Measurement:
self.N = (self.nblocks * self.blocksize)
self.T = self.N / self.samplerate
try:
self.version_major = f.attrs['LASP_VERSION_MAJOR']
self.version_minor = f.attrs['LASP_VERSION_MINOR']
except KeyError:
self.version_major = 0
self.version_minor = 1
# Due to a previous bug, the channel names were not stored
# consistently, i.e. as 'channel_names' and later camelcase.
try:
@ -239,6 +250,7 @@ class Measurement:
self._time = f.attrs['time']
self._qtys = None
try:
qtys_json = f.attrs['qtys']
# Load quantity data
@ -246,9 +258,20 @@ class Measurement:
except KeyError:
# If quantity data is not available, this is an 'old'
# measurement file.
logging.debug(f'Physical quantity data not available in measurement file. Assuming {SIQtys.default}')
self._qtys = [SIQtys.default() for i in range(self.nchannels)]
pass
try:
qtys_enum_idx = f.attrs['qtys_enum_idx']
self._qtys = [SIQtys.fromCppEnum(idx) for idx in qtys_enum_idx]
except KeyError:
pass
if self._qtys is None:
self._qtys = [SIQtys.default() for i in range(self.nchannels)]
logging.debug(f'Physical quantity data not available in measurement file. Assuming {SIQtys.default}')
logging.debug(f'Physical quantity data not available in measurement file. Assuming {SIQtys.default}')
def setAttribute(self, atrname, value):
"""
Set an attribute in the measurement file, and keep a local copy in
@ -276,14 +299,15 @@ class Measurement:
@property
def channelConfig(self):
chcfg = [DaqChannel(channel_enabled=True,
channel_name=chname,
sensitivity=sens,)
for chname, sens in zip(
self.channelNames,
self.sensitivity)]
for i, qty in enumerate(self.qtys):
chcfg[i].qty = qty
chcfg = []
for chname, sens, qty in zip(self.channelNames, self.sensitivity,
self.qtys):
ch = DaqChannel()
ch.enabled = True
ch.name = chname
ch.sensitivity = sens
ch.qty = qty.cpp_enum
chcfg.append(ch)
return chcfg
@channelConfig.setter
@ -292,9 +316,9 @@ class Measurement:
sens = []
qtys = []
for ch in chcfg:
chname.append(ch.channel_name)
chname.append(ch.name)
sens.append(ch.sensitivity)
qtys.append(ch.qty)
qtys.append(SIQtys.fromCppEnum(ch.qty))
self.channelNames = chname
self.sensitivity = sens

7
src/lasp/lasp_record.py
View File

@ -9,7 +9,8 @@ import os
import time
import h5py
import numpy as np
from .lasp_cpp import InDataHandler, StreamMgr
from .lasp_cpp import (InDataHandler, StreamMgr, LASP_VERSION_MAJOR,
LASP_VERSION_MINOR)
from .lasp_atomic import Atomic
@ -119,12 +120,16 @@ class Recording:
f = self.f
f.attrs['LASP_VERSION_MAJOR'] = LASP_VERSION_MAJOR
f.attrs['LASP_VERSION_MINOR'] = LASP_VERSION_MINOR
# Set the bunch of attributes
f.attrs['samplerate'] = daq.samplerate()
f.attrs['nchannels'] = daq.neninchannels()
f.attrs['blocksize'] = blocksize
f.attrs['sensitivity'] = [ch.sensitivity for ch in in_ch]
f.attrs['channelNames'] = [ch.name for ch in in_ch]
# f.attrs['
# Add the start delay here, as firstFrames() is called right after the
# constructor is called.

39
src/lasp/lasp_slm.py
View File

@ -4,7 +4,7 @@
Sound level meter implementation
@author: J.A. de Jong - ASCEE
"""
from .lasp_cpp import SLM as cppSLM
from .lasp_cpp import cppSLM
import numpy as np
from .lasp_common import (TimeWeighting, FreqWeighting, P_REF)
from .filter import SPLFilterDesigner
@ -33,8 +33,8 @@ class SLM:
def __init__(self,
fs,
fbdesigner=None,
tw=TimeWeighting.fast,
fw=FreqWeighting.A,
tw: TimeWeighting =TimeWeighting.fast,
fw: FreqWeighting =FreqWeighting.A,
xmin = None,
xmax = None,
include_overall=True,
@ -50,7 +50,7 @@ class SLM:
fs: Sampling frequency [Hz]
tw: Time Weighting to apply
fw: Frequency weighting to apply
xmin: Filter designator of lowest band
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.
@ -96,36 +96,42 @@ class SLM:
assert fbdesigner.fs == fs
sos_firstx = fbdesigner.createSOSFilter(self.xs[0]).flatten()
self.nom_txt.append(fbdesigner.nominal_txt(self.xs[0]))
sos = np.empty((nfilters, sos_firstx.size), dtype=float, order='C')
sos[0, :] = sos_firstx
sos = np.empty((sos_firstx.size, nfilters), dtype=float, order='C')
sos[:, 0] = sos_firstx
for i, x in enumerate(self.xs[1:]):
sos[i+1, :] = fbdesigner.createSOSFilter(x).flatten()
sos[:, i+1] = 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)
# a0 = 1, b0 = 1, rest is zero
sos[-1,:] = sos_overall
sos[:,-1] = sos_overall
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.
sos = sos_overall[np.newaxis,:]
sos = sos_overall[:,np.newaxis]
self.nom_txt.append('overall')
self.slm = cppSLM(fs, prefilter, sos,
fs, tw[0], level_ref_value)
# Downsampling factor, determine from single pole low pass filter time
# constant, such that aliasing is ~ allowed at 20 dB lower value
# and
dsfac = cppSLM.suggestedDownSamplingFac(fs, tw[0])
dsfac = self.slm.downsampling_fac
if dsfac > 0:
# Not unfiltered data
self.fs_slm = fs / self.slm.downsampling_fac
if prefilter is not None:
self.slm = cppSLM.fromBiquads(fs, level_ref_value, dsfac,
tw[0],
prefilter, sos)
else:
self.fs_slm = fs
self.slm = cppSLM.fromBiquads(fs, level_ref_value, dsfac,
tw[0],
sos)
self.fs_slm = fs / dsfac
# Initialize counter to 0
self.N = 0
@ -189,7 +195,6 @@ class SLM:
output['mid'] = self.fbdesigner.fm(list(self.xs))
logging.debug(list(self.xs))
logging.debug(output['mid'])
# Indiced at 0, as pyxSLM always returns 2D arrays.
if self.include_overall and self.fbdesigner is not None:
output['overall'] = dat[-1,0]

2
src/lasp/pybind11/lasp_daqconfig.cpp
View File

@ -69,6 +69,8 @@ void init_daqconfiguration(py::module &m) {
.value("UserDefined", DaqChannel::Qty::UserDefined);
daqchannel.def_readwrite("qty", &DaqChannel::qty);
daqchannel.def("__eq__", [](const DaqChannel& a, const DaqChannel& b) { return a==b;});
/// DaqConfiguration
daqconfig.def(py::init<>());
daqconfig.def(py::init<const DeviceInfo &>());

8
src/lasp/pybind11/lasp_dsp_pybind.cpp
View File

@ -69,15 +69,16 @@ void init_dsp(py::module &m) {
/// AvPowerSpectra
py::class_<AvPowerSpectra> aps(m, "AvPowerSpectra");
aps.def(py::init<const us, const Window::WindowType, const d, const int>(),
py::arg("nfft"), py::arg("WindowType"), py::arg("overlap_percentage"),
py::arg("time_constant"));
py::arg("nfft") = 2048, py::arg("windowType") =Window::WindowType::Hann, py::arg("overlap_percentage") = 50.0,
py::arg("time_constant") = -1);
aps.def("compute", [](AvPowerSpectra &aps, const dmat &timedata) {
std::optional<arma::cx_cube> res = aps.compute(timedata);
return res.value_or(arma::cx_cube(0,0,0));
}
);
py::class_<SLM> slm(m, "SLM");
py::class_<SLM> slm(m, "cppSLM");
slm.def_static(
"fromBiquads",
@ -94,5 +95,6 @@ void init_dsp(py::module &m) {
slm.def("Lpeak", &SLM::Lpeak);
slm.def("Leq", &SLM::Leq);
slm.def("Lmax", &SLM::Lmax);
slm.def_static("suggestedDownSamplingFac", &SLM::suggestedDownSamplingFac);
}
/** @} */

16
src/lasp/pybind11/lasp_streammgr.cpp
View File

@ -2,6 +2,7 @@
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <functional>
#include <stdint.h>
using std::cerr;
@ -15,14 +16,14 @@ void init_streammgr(py::module &m) {
m, "StreamMgr");
py::enum_<StreamMgr::StreamType>(smgr, "StreamType")
.value("input", StreamMgr::StreamType::input)
.value("output", StreamMgr::StreamType::output);
.value("input", StreamMgr::StreamType::input)
.value("output", StreamMgr::StreamType::output);
smgr.def("startStream", &StreamMgr::startStream);
smgr.def("stopStream", &StreamMgr::stopStream);
smgr.def_static("getInstance", []() {
return std::unique_ptr<StreamMgr, py::nodelete>(&StreamMgr::getInstance());
});
return std::unique_ptr<StreamMgr, py::nodelete>(&StreamMgr::getInstance());
});
smgr.def("stopAllStreams", &StreamMgr::stopAllStreams);
smgr.def("setSiggen", &StreamMgr::setSiggen);
@ -31,4 +32,11 @@ void init_streammgr(py::module &m) {
smgr.def("isStreamRunningOK", &StreamMgr::isStreamRunningOK);
smgr.def("isStreamRunning", &StreamMgr::isStreamRunning);
smgr.def("getDaq", &StreamMgr::getDaq, py::return_value_policy::reference);
smgr.def("rescanDAQDevices", [](StreamMgr& smgr, bool background) {
// A pure C++ callback is the second argument to rescanDAQDevices, which
// cannot be wrapped to Pybind11. Only the one without callback is
// forwarded here to Python code.
smgr.rescanDAQDevices(background);
},
py::arg("background")=false);
}

39
test/test_cppslm.py
View File

@ -1,6 +1,6 @@
#!/usr/bin/python3
import numpy as np
from lasp import SLM
from lasp import cppSLM
from lasp.filter import SPLFilterDesigner
import matplotlib.pyplot as plt
@ -11,7 +11,7 @@ def test_cppslm1():
fs = 48000
omg = 2*np.pi*1000
slm = SLM.fromBiquads(fs, 2e-5, 1, 0.125, [1.,0,0,1,0,0])
slm = cppSLM.fromBiquads(fs, 2e-5, 1, 0.125, [1.,0,0,1,0,0])
t = np.linspace(0, 10, 10*fs, endpoint=False)
@ -40,7 +40,7 @@ def test_cppslm2():
omg = 2*np.pi*1000
filt = SPLFilterDesigner(fs).A_Sos_design()
slm = SLM.fromBiquads(fs, 2e-5, 0, 0.125, filt.flatten(), [1.,0,0,1,0,0])
slm = cppSLM.fromBiquads(fs, 2e-5, 0, 0.125, filt.flatten(), [1.,0,0,1,0,0])
t = np.linspace(0, 10, 10*fs, endpoint=False)
@ -60,16 +60,33 @@ def test_cppslm2():
# (Fast, Slow etc)
assert np.isclose(out[-1,0], level, atol=1e-2)
def test_cppslm3():
fs = 48000
omg = 2*np.pi*1000
filt = SPLFilterDesigner(fs).A_Sos_design()
slm = cppSLM.fromBiquads(fs, 2e-5, 0, 0.125, filt.flatten(), [1.,0,0,1,0,0])
t = np.linspace(0, 10, 10*fs, endpoint=False)
in_ = 10*np.sin(omg*t) * np.sqrt(2)+np.random.randn()
# Compute overall RMS
rms = np.sqrt(np.sum(in_**2)/in_.size)
# Compute overall level
level = 20*np.log10(rms/2e-5)
# Output of SLM
out = slm.run(in_)
Lpeak = 20*np.log10(np.max(np.abs(in_)/2e-5))
Lpeak
slm.Lpeak()
assert np.isclose(out[-1,0], slm.Leq()[0][0], atol=1e-2)
assert np.isclose(Lpeak, slm.Lpeak()[0][0], atol=1e0)
if __name__ == '__main__':
test_cppslm1()
test_cppslm2()
# plt.plot(t,out[:,0])
# plt.close('all')
# from scipy.signal import sosfreqz
# omg, H = sosfreqz(filt)
# freq = omg / (2*np.pi)*fs
# plt.figure()
# plt.semilogx(freq[1:],20*np.log10(np.abs(H[1:])))
test_cppslm3()