Merge branch 'master' into develop

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;
 
 }
 

src/lasp/dsp/lasp_biquadbank.cpp

@@ -2,23 +2,25 @@
 #include "lasp_biquadbank.h"
 #include "debugtrace.hpp"
 #include "lasp_thread.h"
+#include <cassert>
 #include <vector>
 
 using std::cerr;
 using std::endl;
 using rte = std::runtime_error;
+using lock = std::scoped_lock<std::mutex>;
 
 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<std::future<vd>> 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) {
+          [&](vd inout, us i) {
-          // Copy column
+          _filters[i].filter(inout);
-          vd col = inout;
+          return inout;
-          _filters[i].filter(col);
-          return col;
           }, // 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<Filter> BiquadBank::clone() const {
+  lock lck(_mtx);
   return std::make_unique<BiquadBank>(_filters, _gains);
 }

src/lasp/dsp/lasp_biquadbank.h

@@ -49,7 +49,7 @@ public:
 class BiquadBank : public Filter {
   std::vector<SeriesBiquad> _filters;
   vd _gains;
-
+  mutable std::mutex _mtx;
 
 public:
   /**

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()) {
+    if (nchannels > _freqWeightingFilters.size()) {
-      while (nchannels > _freqWeightingFilter.size()) {
+      while (nchannels > _freqWeightingFilters.size()) {
-        _freqWeightingFilter.emplace_back(_filterPrototype->clone());
+        _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)

src/lasp/dsp/lasp_rtaps.h

@@ -23,7 +23,7 @@ class RtAps : public ThreadedInDataHandler {
 
   std::mutex _mtx;
   std::unique_ptr<Filter> _filterPrototype;
-  std::vector<std::unique_ptr<Filter>> _freqWeightingFilter;
+  std::vector<std::unique_ptr<Filter>> _freqWeightingFilters;
 
   AvPowerSpectra _ps;
 

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 <cassert>
 #include <type_traits>
 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) {
+    for (auto f : _filters) {
-      _filter->filter(signal_dynamic);
+      assert(f.second);
+      f.second->filter(signal_dynamic);
     }
     signal += signal_dynamic;
   }
 
   return signal;
 }
-void Siggen::setFilter(std::shared_ptr<Filter> filter) {
+void Siggen::setFilter(const std::string &name,
+                       std::shared_ptr<Filter> 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();
 }

src/lasp/dsp/lasp_siggen.h

@@ -22,7 +22,7 @@ class DaqData;
  */
 class Siggen {
 private:
-  std::shared_ptr<Filter> _filter;
+  std::map<std::string, std::shared_ptr<Filter>> _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
+   * @brief Set a filter on the signal. For example to EQ the signal, or
-   * otherwise to shape the spectrum.
+   * 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<Filter> f);
+  void setFilter(const std::string& name, std::shared_ptr<Filter> f);
 
   /**
    * @brief Set a linear DC offset value to the signal

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:
+        if data.ndim == 2 and data.shape[1] != 1:
-            data = data[:, None]
+            raise RuntimeError("invalid number of channels, should be 1")
-        elif data.ndim == 2:
-            assert data.shape[1] == 1, "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)
+            channel = self.designer.createSOSFilter(x)
-            self._fb.setFilter(i, sos)
+            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 = {}

src/lasp/pybind11/lasp_dsp_pybind.cpp

@@ -79,11 +79,20 @@ void init_dsp(py::module &m) {
   });
 
   /// BiquadBank
-  py::class_<BiquadBank, std::shared_ptr<BiquadBank>> bqb(m, "BiquadBank");
+  py::class_<BiquadBank, Filter, std::shared_ptr<BiquadBank>> bqb(m,
-  bqb.def(py::init<const dmat &, const vd *>());
+                                                                  "BiquadBank");
+  bqb.def(py::init([](dpyarray coefs) {
+    return std::make_shared<BiquadBank>(NpyToMat<d, false>(coefs));
+  }));
+  bqb.def(py::init([](dpyarray coefs, dpyarray gains) {
+    vd gains_arma = NpyToMat<d, false>(gains);
+    return std::make_shared<BiquadBank>(NpyToMat<d, false>(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<d, true>(input);
     b.filter(inout);
     return ColToNpy(inout);

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,std::shared_ptr<Noise>> noise(m, "Noise", siggen);