diff --git a/src/lasp/dsp/lasp_avpowerspectra.cpp b/src/lasp/dsp/lasp_avpowerspectra.cpp index 7832f66..83c7f05 100644 --- a/src/lasp/dsp/lasp_avpowerspectra.cpp +++ b/src/lasp/dsp/lasp_avpowerspectra.cpp @@ -81,6 +81,7 @@ AvPowerSpectra::AvPowerSpectra(const us nfft, const Window::WindowType w, void AvPowerSpectra::reset() { _timeBuf.reset(); _est.reset(); + n_averages=0; } diff --git a/src/lasp/dsp/lasp_biquadbank.cpp b/src/lasp/dsp/lasp_biquadbank.cpp index 6ee8060..115f6e3 100644 --- a/src/lasp/dsp/lasp_biquadbank.cpp +++ b/src/lasp/dsp/lasp_biquadbank.cpp @@ -2,23 +2,25 @@ #include "lasp_biquadbank.h" #include "debugtrace.hpp" #include "lasp_thread.h" +#include #include using std::cerr; using std::endl; using rte = std::runtime_error; +using lock = std::scoped_lock; SeriesBiquad::SeriesBiquad(const vd &filter_coefs) { DEBUGTRACE_ENTER; if (filter_coefs.n_cols != 1) { throw rte("Expected filter coefficients for a single SeriesBiquad as a " - "single column with length 6 x n_filters"); + "single column with length 6 x n_filters"); } if (filter_coefs.n_rows % 6 != 0) { cerr << "Number of rows given: " << filter_coefs.n_rows << endl; throw rte("filter_coefs should be multiple of 6, given: " + - std::to_string(filter_coefs.n_rows)); + std::to_string(filter_coefs.n_rows)); } us nfilters = filter_coefs.n_rows / 6; @@ -32,7 +34,7 @@ SeriesBiquad::SeriesBiquad(const vd &filter_coefs) { /// Check if third row in this matrix equals unity. if (!arma::approx_equal(sos.row(3), arma::rowvec(nfilters, arma::fill::ones), - "absdiff", 1e-9)) { + "absdiff", 1e-9)) { std::cerr << "Read row: " << sos.row(3) << endl; throw rte( @@ -82,6 +84,7 @@ BiquadBank::BiquadBank(const dmat &filters, const vd *gains) { * @brief Make sure the pool is created once, such that all threads are ready * for use. */ + lock lck(_mtx); getPool(); for (us i = 0; i < filters.n_cols; i++) { @@ -96,6 +99,8 @@ BiquadBank::BiquadBank(const dmat &filters, const vd *gains) { } void BiquadBank::setGains(const vd &gains) { DEBUGTRACE_ENTER; + lock lck(_mtx); + const us nfilters = _filters.size(); if (gains.size() != nfilters) { throw rte("Invalid number of gain values given."); @@ -105,34 +110,46 @@ void BiquadBank::setGains(const vd &gains) { void BiquadBank::filter(vd &inout) { - dmat res(inout.n_rows, _filters.size()); + lock lck(_mtx); std::vector> futs; +#if 1 auto &pool = getPool(); - for (us i = 0; i < res.n_cols; i++) { + vd inout_cpy = inout; + for (us i = 0; i < _filters.size(); i++) { futs.emplace_back(pool.submit( - [&](us i) { - // Copy column - vd col = inout; - _filters[i].filter(col); - return col; + [&](vd inout, us i) { + _filters[i].filter(inout); + return inout; }, // Launch a task to filter. - i // Column i as argument to the lambda function above. + inout_cpy, i // Column i as argument to the lambda function above. )); } // Zero-out in-out and sum-up the filtered values inout.zeros(); - for (us i = 0; i < res.n_cols; i++) { + for (us i = 0; i < _filters.size(); i++) { inout += futs[i].get() * _gains[i]; } +#else + /// Testing, unthreaded version + vd inout_cpy = inout; + inout.zeros(); + for (us i = 0; i < _filters.size(); i++) { + vd cpy_again = inout_cpy; + _filters[i].filter(cpy_again); + inout += cpy_again * _gains(i); + } +#endif } void BiquadBank::reset() { DEBUGTRACE_ENTER; + lock lck(_mtx); for (auto &f : _filters) { f.reset(); } } std::unique_ptr BiquadBank::clone() const { + lock lck(_mtx); return std::make_unique(_filters, _gains); } diff --git a/src/lasp/dsp/lasp_biquadbank.h b/src/lasp/dsp/lasp_biquadbank.h index 0fd6fab..c37df12 100644 --- a/src/lasp/dsp/lasp_biquadbank.h +++ b/src/lasp/dsp/lasp_biquadbank.h @@ -49,7 +49,7 @@ public: class BiquadBank : public Filter { std::vector _filters; vd _gains; - + mutable std::mutex _mtx; public: /** diff --git a/src/lasp/dsp/lasp_rtaps.cpp b/src/lasp/dsp/lasp_rtaps.cpp index d7d2de1..4f3af10 100644 --- a/src/lasp/dsp/lasp_rtaps.cpp +++ b/src/lasp/dsp/lasp_rtaps.cpp @@ -34,21 +34,21 @@ bool RtAps::inCallback_threaded(const DaqData &data) { if (_filterPrototype) { // Adjust number of filters, if necessary - if (nchannels > _freqWeightingFilter.size()) { - while (nchannels > _freqWeightingFilter.size()) { - _freqWeightingFilter.emplace_back(_filterPrototype->clone()); + if (nchannels > _freqWeightingFilters.size()) { + while (nchannels > _freqWeightingFilters.size()) { + _freqWeightingFilters.emplace_back(_filterPrototype->clone()); } - for (auto &filter : _freqWeightingFilter) { + for (auto &filter : _freqWeightingFilters) { filter->reset(); } } // Apply filtering - #pragma omp parallel for + /* #pragma omp parallel for */ for (us i = 0; i < nchannels; i++) { vd col = fltdata.col(i); - _freqWeightingFilter.at(i)->filter(col); + _freqWeightingFilters.at(i)->filter(col); fltdata.col(i) = col; } } // End of if(_filterPrototype) diff --git a/src/lasp/dsp/lasp_rtaps.h b/src/lasp/dsp/lasp_rtaps.h index 26a0510..9443c75 100644 --- a/src/lasp/dsp/lasp_rtaps.h +++ b/src/lasp/dsp/lasp_rtaps.h @@ -23,7 +23,7 @@ class RtAps : public ThreadedInDataHandler { std::mutex _mtx; std::unique_ptr _filterPrototype; - std::vector> _freqWeightingFilter; + std::vector> _freqWeightingFilters; AvPowerSpectra _ps; diff --git a/src/lasp/dsp/lasp_siggen.cpp b/src/lasp/dsp/lasp_siggen.cpp index aa747e7..361e7be 100644 --- a/src/lasp/dsp/lasp_siggen.cpp +++ b/src/lasp/dsp/lasp_siggen.cpp @@ -1,11 +1,12 @@ /* #define DEBUGTRACE_ENABLED */ -#include "debugtrace.hpp" #include "lasp_siggen.h" +#include "debugtrace.hpp" #include "lasp_mathtypes.h" #include #include using std::cerr; using std::endl; +using rte = std::runtime_error; inline d level_amp(d level_dB) { return pow(10, level_dB / 20); } @@ -21,18 +22,24 @@ vd Siggen::genSignal(const us nframes) { if (!_muted) { vd signal_dynamic = _level_linear * genSignalUnscaled(nframes); - if (_filter) { - _filter->filter(signal_dynamic); + for (auto f : _filters) { + assert(f.second); + f.second->filter(signal_dynamic); } signal += signal_dynamic; } return signal; } -void Siggen::setFilter(std::shared_ptr filter) { +void Siggen::setFilter(const std::string &name, + std::shared_ptr filter) { DEBUGTRACE_ENTER; mutexlock lck(_mtx); - _filter = filter; + if (filter) { + _filters[name] = filter; + } else if (_filters.find(name) != _filters.end()) { + _filters.extract(name); + } } void Siggen::setDCOffset(const d offset) { DEBUGTRACE_ENTER; @@ -48,5 +55,9 @@ void Siggen::reset(const d newFs) { DEBUGTRACE_ENTER; mutexlock lck(_mtx); _fs = newFs; + for (auto &f : _filters) { + assert(f.second); + f.second->reset(); + } resetImpl(); } diff --git a/src/lasp/dsp/lasp_siggen.h b/src/lasp/dsp/lasp_siggen.h index 969ade9..1cd2e16 100644 --- a/src/lasp/dsp/lasp_siggen.h +++ b/src/lasp/dsp/lasp_siggen.h @@ -22,7 +22,7 @@ class DaqData; */ class Siggen { private: - std::shared_ptr _filter; + std::map> _filters; d _dc_offset = 0, _level_linear = 1; bool _muted = false; @@ -37,12 +37,14 @@ public: virtual ~Siggen() = default; /** - * @brief Set a post-filter on the signal. For example to EQ the signal, or - * otherwise to shape the spectrum. + * @brief Set a filter on the signal. For example to EQ the signal, or + * otherwise to shape the spectrum. Filters are stored in a map, and + * setFilter can be used to overwrite them, if a null pointer is given. * + * @param name The name of the filter (i.e. "noiseshaper", or "eq" * @param f The filter to install. */ - void setFilter(std::shared_ptr f); + void setFilter(const std::string& name, std::shared_ptr f); /** * @brief Set a linear DC offset value to the signal diff --git a/src/lasp/lasp_octavefilter.py b/src/lasp/lasp_octavefilter.py index fdcb39e..632cd78 100644 --- a/src/lasp/lasp_octavefilter.py +++ b/src/lasp/lasp_octavefilter.py @@ -12,6 +12,7 @@ __all__ = ['FirOctaveFilterBank', 'FirThirdOctaveFilterBank', from .filter.filterbank_design import (OctaveBankDesigner, ThirdOctaveBankDesigner) +from .lasp_cpp import BiquadBank import numpy as np @@ -33,12 +34,11 @@ class OverallFilterBank: """ Filter input data """ - if data.ndim == 1: - data = data[:, None] - elif data.ndim == 2: - assert data.shape[1] == 1, "invalid number of channels, should be 1" + if data.ndim == 2 and data.shape[1] != 1: + raise RuntimeError("invalid number of channels, should be 1") if data.shape[0] == 0: + # No new time samples return {} # Output given as a dictionary with x as the key @@ -235,12 +235,20 @@ class SosFilterBank: self.fs = fs self.xs = list(range(xmin, xmax + 1)) + + # The number of parallel filters nfilt = len(self.xs) + self.nfilt = nfilt - self._fb = pyxSosFilterBank(nfilt, 5) + + sos = None for i, x in enumerate(self.xs): - sos = self.designer.createSOSFilter(x) - self._fb.setFilter(i, sos) + channel = self.designer.createSOSFilter(x) + if sos is None: + sos = np.empty((channel.size, len(self.xs))) + sos[:, i] = channel.flatten() + + self._fb = BiquadBank(sos) self.xmin = xmin self.xmax = xmax @@ -250,13 +258,14 @@ class SosFilterBank: """ Filter input data """ - assert data.ndim == 2 - assert data.shape[1] == 1, "invalid number of channels, should be 1" + + if data.ndim > 1 and data.shape[1] != 1: + raise RuntimeError("invalid number of channels, should be 1") if data.shape[0] == 0: return {} - filtered_data = self._fb.filter_(data) + filtered_data = self._fb.filter(data) # Output given as a dictionary with nom_txt as the key output = {} diff --git a/src/lasp/pybind11/lasp_dsp_pybind.cpp b/src/lasp/pybind11/lasp_dsp_pybind.cpp index e8cdfea..48d0273 100644 --- a/src/lasp/pybind11/lasp_dsp_pybind.cpp +++ b/src/lasp/pybind11/lasp_dsp_pybind.cpp @@ -79,11 +79,20 @@ void init_dsp(py::module &m) { }); /// BiquadBank - py::class_> bqb(m, "BiquadBank"); - bqb.def(py::init()); + py::class_> bqb(m, + "BiquadBank"); + bqb.def(py::init([](dpyarray coefs) { + return std::make_shared(NpyToMat(coefs)); + })); + bqb.def(py::init([](dpyarray coefs, dpyarray gains) { + vd gains_arma = NpyToMat(gains); + return std::make_shared(NpyToMat(coefs), &gains_arma); + })); + bqb.def("setGains", [](BiquadBank &b, dpyarray gains) { b.setGains(NpyToCol(gains)); }); bqb.def("filter", [](BiquadBank &b, dpyarray input) { + // Yes, a copy here vd inout = NpyToCol(input); b.filter(inout); return ColToNpy(inout); diff --git a/src/lasp/pybind11/lasp_siggen.cpp b/src/lasp/pybind11/lasp_siggen.cpp index ecad69f..9a6e424 100644 --- a/src/lasp/pybind11/lasp_siggen.cpp +++ b/src/lasp/pybind11/lasp_siggen.cpp @@ -32,7 +32,7 @@ void init_siggen(py::module &m) { siggen.def("setMute", &Siggen::setMute); siggen.def("setLevel", &Siggen::setLevel); siggen.def("setLevel", &Siggen::setLevel, py::arg("newLevel"), py::arg("dB") = true); - /* siggen.def("genSignal", &Siggen::genSignal); */ + siggen.def("setFilter", &Siggen::setFilter); py::class_> noise(m, "Noise", siggen);