Measurement.fromnpy(): accept sensitivity as scalar or 0-dim numpy.ndarray

cpp_src/device/CMakeLists.txt

@@ -10,8 +10,12 @@ add_library(lasp_device_lib OBJECT
     lasp_rtaudiodaq.cpp
     lasp_streammgr.cpp
     lasp_indatahandler.cpp
+    lasp_uldaq.cpp
+    uldaq/lasp_uldaq_impl.cpp
+    uldaq/lasp_uldaq_bufhandler.cpp
+    uldaq/lasp_uldaq_common.cpp
     portaudio/lasp_portaudiodaq.cpp
-)
+    )
 
 # Callback requires certain arguments that are not used by code. This disables
 # a compiler warning about it.
@@ -24,9 +28,7 @@ target_include_directories(lasp_device_lib INTERFACE
     ${CMAKE_CURRENT_SOURCE_DIR})
 
 if(LASP_HAS_ULDAQ)
-    add_subdirectory(uldaq)
-    target_include_directories(lasp_device_lib INTERFACE uldaq)
-    target_link_libraries(lasp_device_lib uldaq_backend uldaq)
+    target_link_libraries(lasp_device_lib uldaq)
 endif()
 if(LASP_HAS_RTAUDIO)
     target_link_libraries(lasp_device_lib rtaudio)

cpp_src/device/uldaq/CMakeLists.txt

@@ -1,6 +0,0 @@
-
-add_library(uldaq_backend lasp_uldaq.cpp lasp_uldaq_bufhandler.cpp lasp_uldaq_common.cpp lasp_uldaq_impl.cpp)
-target_include_directories(uldaq_backend PUBLIC ../)
-target_include_directories(uldaq_backend PUBLIC ../../)
-target_include_directories(uldaq_backend PUBLIC ../../dsp)
-target_include_directories(uldaq_backend INTERFACE ${CMAKE_CURRENT_SOURCE_DIR})

cpp_src/device/uldaq/lasp_uldaq_common.cpp

@@ -3,8 +3,8 @@
 #include "lasp_config.h"
 
 #if LASP_HAS_ULDAQ == 1
-#include "lasp_daq.h"
 #include "lasp_uldaq_common.h"
+#include "lasp_daq.h"
 
 string getErrMsg(UlError err) {
   string errstr;
@@ -21,9 +21,11 @@ void showErr(string errstr) {
   std::cerr << "***********************************************\n\n";
 }
 void showErr(UlError err) {
-  if (err != ERR_NO_ERROR) showErr(getErrMsg(err));
+  if (err != ERR_NO_ERROR)
+    showErr(getErrMsg(err));
 }
 
+
 void throwOnPossibleUlException(UlError err) {
   if (err == ERR_NO_ERROR) {
     return;

cpp_src/dsp/lasp_avpowerspectra.cpp

@@ -10,8 +10,7 @@ PowerSpectra::PowerSpectra(const us nfft, const Window::WindowType w)
     : PowerSpectra(Window::create(w, nfft)) {}
 
 PowerSpectra::PowerSpectra(const vd &window)
-    : nfft(window.size()), _fft(nfft), _window(window)
-{
+    : nfft(window.size()), _fft(nfft), _window(window) {
 
   d win_pow = arma::sum(window % window) / window.size();
 
@@ -22,8 +21,7 @@ PowerSpectra::PowerSpectra(const vd &window)
   DEBUGTRACE_PRINT(win_pow);
 }
 
-ccube PowerSpectra::compute(const dmat &input)
-{
+ccube PowerSpectra::compute(const dmat &input) {
 
   /// Run very often. Silence this one.
   /* DEBUGTRACE_ENTER; */
@@ -36,13 +34,11 @@ ccube PowerSpectra::compute(const dmat &input)
   cmat rfft = _fft.fft(input_tmp);
 
   ccube output(rfft.n_rows, input.n_cols, input.n_cols);
-  for (us i = 0; i < input.n_cols; i++)
-  {
+  for (us i = 0; i < input.n_cols; i++) {
     /// This one can be run in parallel without any problem. Note that it is
     /// the inner loop that is run in parallel.
 #pragma omp parallel for
-    for (us j = 0; j < input.n_cols; j++)
-    {
+    for (us j = 0; j < input.n_cols; j++) {
       output.slice(j).col(i) =
           _scale_fac * (rfft.col(i) % arma::conj(rfft.col(j)));
 
@@ -56,46 +52,37 @@ ccube PowerSpectra::compute(const dmat &input)
 AvPowerSpectra::AvPowerSpectra(const us nfft, const Window::WindowType w,
                                const d overlap_percentage,
                                const d time_constant_times_fs)
-    : _ps(nfft, w)
-{
+    : _ps(nfft, w) {
 
   DEBUGTRACE_ENTER;
-  if (overlap_percentage >= 100 || overlap_percentage < 0)
-  {
+  if (overlap_percentage >= 100 || overlap_percentage < 0) {
     throw rte("Overlap percentage should be >= 0 and < 100");
   }
 
   _overlap_keep = (nfft * overlap_percentage) / 100;
   DEBUGTRACE_PRINT(_overlap_keep);
-  if (_overlap_keep >= nfft)
-  {
+  if (_overlap_keep >= nfft) {
     throw rte("Overlap is too high. Results in no jump. Please "
               "choose a smaller overlap percentage or a higher nfft");
   }
-  if (time_constant_times_fs < 0)
-  {
+  if (time_constant_times_fs < 0) {
     _mode = Mode::Averaging;
-  }
-  else if (time_constant_times_fs == 0)
-  {
+  } else if (time_constant_times_fs == 0) {
     _mode = Mode::Spectrogram;
-  }
-  else
-  {
+  } else {
     _mode = Mode::Leaking;
-    _alpha = d_exp(-static_cast<d>((nfft - _overlap_keep) / time_constant_times_fs));
+    _alpha = d_exp(-static_cast<d>((nfft - _overlap_keep)/time_constant_times_fs));
     DEBUGTRACE_PRINT(_alpha);
   }
 }
-void AvPowerSpectra::reset()
-{
+void AvPowerSpectra::reset() {
   _timeBuf.reset();
   _est.reset();
-  _n_averages = 0;
+  n_averages=0;
+
 }
 
-ccube AvPowerSpectra::compute(const dmat &timedata)
-{
+ccube AvPowerSpectra::compute(const dmat &timedata) {
 
   DEBUGTRACE_ENTER;
   DEBUGTRACE_PRINT(timedata.n_rows);
@@ -104,42 +91,32 @@ ccube AvPowerSpectra::compute(const dmat &timedata)
   _timeBuf.push(timedata);
 
   bool run_once = false;
-  while (_timeBuf.size() >= _ps.nfft)
-  {
+  while (_timeBuf.size() >= _ps.nfft) {
 
     DEBUGTRACE_PRINT((int)_mode);
 
     dmat samples = _timeBuf.pop(_ps.nfft, _overlap_keep);
 
-    switch (_mode)
-    {
+    switch (_mode) {
     case (Mode::Spectrogram):
       DEBUGTRACE_PRINT("Spectrogram mode");
       _est = _ps.compute(samples);
-
       break;
     case (Mode::Averaging):
       DEBUGTRACE_PRINT("Averaging mode");
-      _n_averages++;
-      if (_n_averages == 1)
-      {
+      n_averages++;
+      if (n_averages == 1) {
         _est = _ps.compute(samples);
-
-      }
-      else
-      {
-        _est = (static_cast<d>(_n_averages - 1) / _n_averages) * _est +
-               _ps.compute(samples) / _n_averages;
+      } else {
+        _est = (static_cast<d>(n_averages - 1) / n_averages) * _est +
+               _ps.compute(samples) / n_averages;
       }
       break;
     case (Mode::Leaking):
       DEBUGTRACE_PRINT("Leaking mode");
-      if (arma::size(_est) == arma::size(0, 0, 0))
-      {
+      if (arma::size(_est) == arma::size(0, 0, 0)) {
         _est = _ps.compute(samples);
-      }
-      else
-      {
+      } else {
         _est = _alpha * _est + (1 - _alpha) * _ps.compute(samples);
       }
       break;
@@ -151,136 +128,6 @@ ccube AvPowerSpectra::compute(const dmat &timedata)
   /// Othewise, we return an empty ccube.
   return run_once ? _est : ccube();
 }
-ccube AvPowerSpectra::get_est() const
-{
-  return _est;
-}
-
-AvSweepPowerSpectra::AvSweepPowerSpectra(const us nfft, const Window::WindowType w,
-                                         const d overlap_percentage,
-                                         const d time_constant_times_fs)
-    : _nfft(nfft), _ps(nfft, w)
-{
-
-  DEBUGTRACE_ENTER;
-  if (overlap_percentage >= 100 || overlap_percentage < 0)
-  {
-    throw rte("Overlap percentage should be >= 0 and < 100");
-  }
-
-  _overlap_keep = (nfft * overlap_percentage) / 100;
-  DEBUGTRACE_PRINT(_overlap_keep);
-  if (_overlap_keep >= nfft)
-  {
-    throw rte("Overlap is too high. Results in no jump. Please "
-              "choose a smaller overlap percentage or a higher nfft");
-  }
-  if (time_constant_times_fs < 0)
-  {
-    _mode = Mode::Averaging;
-  }
-  else if (time_constant_times_fs == 0)
-  {
-    _mode = Mode::Spectrogram;
-  }
-  else
-  {
-    _mode = Mode::Leaking;
-    _alpha = d_exp(-static_cast<d>((nfft - _overlap_keep) / time_constant_times_fs));
-    DEBUGTRACE_PRINT(_alpha);
-  }
-}
-void AvSweepPowerSpectra::reset()
-{
-  _timeBuf.reset();
-  _est.reset();
-  _n_averages.reset();
-}
-
-ccube AvSweepPowerSpectra::compute(const dmat &timedata)
-{
-
-  DEBUGTRACE_ENTER;
-  DEBUGTRACE_PRINT(timedata.n_rows);
-  DEBUGTRACE_PRINT(_ps.nfft);
-
-  _timeBuf.push(timedata);
-
-  bool run_once = false;
-  while (_timeBuf.size() >= _ps.nfft)
-  {
-
-    DEBUGTRACE_PRINT((int)_mode);
-
-    dmat samples = _timeBuf.pop(_ps.nfft, _overlap_keep);
-
-    us ncols = timedata.n_cols;
-    
-    switch (_mode)
-    {
-    case (Mode::Spectrogram):
-      DEBUGTRACE_PRINT("Spectrogram mode");
-      _est = _ps.compute(samples);
-      break;
-    case (Mode::Averaging):
-    {
-      DEBUGTRACE_PRINT("Averaging mode");
-
-      ccube temp = _ps.compute(samples);
-      us nfreq = arma::size(temp)[0];
-
-      if (_est.size() == 0) {
-        // Initialize empty
-        _est = ccube(nfreq, ncols, ncols, arma::fill::zeros);
-        _n_averages = vd(nfreq, arma::fill::zeros);
-      }
-            
-      // d threshold = arma::mean(arma::real(temp.slice(0).col(0)));
-      d threshold = 1e-13;
-      
-      for (us f=0; f < nfreq; f++)
-      {
-        if (real(temp(f, 0, 0)) > threshold)
-        {
-          _n_averages[f]++;
-          if (_n_averages[f] == 1)
-          {
-            _est.row(f) = temp.row(f);
-            // _est.subcube(f, 0, 0, f, ncols - 1, ncols - 1) = temp.subcube(f, 0, 0, f, ncols - 1, ncols - 1);
-          }
-          else
-          {
-            _est.row(f) = (static_cast<d>(_n_averages[f] - 1) / _n_averages[f]) * _est.row(f) +
-                                                             temp.row(f) / _n_averages[f];
-            // _est.subcube(f, 0, 0, f, ncols - 1, ncols - 1) = (static_cast<d>(_n_averages[f] - 1) / _n_averages[f]) * _est.subcube(f, 0, 0, f, ncols - 1, ncols - 1) +
-                                                            //  temp.subcube(f, 0, 0, f, ncols - 1, ncols - 1) / _n_averages[f];
-          }
-        }
-      }
-    }
-    break;
-
-    case (Mode::Leaking):
-      DEBUGTRACE_PRINT("Leaking mode");
-      if (arma::size(_est) == arma::size(0, 0, 0))
-      {
-        _est = _ps.compute(samples);
-      }
-      else
-      {
-        _est = _alpha * _est + (1 - _alpha) * _ps.compute(samples);
-      }
-      break;
-    } // end switch mode
-    run_once = true;
-  }
- 
-  /// Return a copy of current estimator in case we have done one update.
-  /// Othewise, we return an empty ccube.
-  return run_once ? _est : ccube();
-}
-
-ccube AvSweepPowerSpectra::get_est() const
-{
+ccube AvPowerSpectra::get_est() const {
   return _est;
 }

cpp_src/dsp/lasp_avpowerspectra.h

@@ -81,7 +81,7 @@ class AvPowerSpectra {
 
   Mode _mode;
   d _alpha;          // Only valid in case of 'Leaking'
-  us _n_averages = 0; // Only valid in case of 'Averaging'
+  us n_averages = 0; // Only valid in case of 'Averaging'
 
   PowerSpectra _ps;
   /**
@@ -165,107 +165,3 @@ public:
   us exactOverlapSamples() const { return _ps.nfft - _overlap_keep; }
 };
 /** @} */
-
-
-/**
- * @brief Estimate cross-power spectra using Welch' method of spectral
- * estimation. The exact amount of overlap in Welch' method is rounded up to a
- * certain amount of samples. This class is specifically for sweeps and will ignore
- * fft blocks where the auto-power of the "reference" signal is below the average 
- * when in averaging mode.
- */
-class AvSweepPowerSpectra {
-
-  enum class Mode {
-    Averaging = 0,  // Averaging all time date
-    Leaking = 1,    // Exponential weighting of an "instantaneous cps"
-    Spectrogram = 2 // Instantenous spectrum, no averaging
-  };
-
-  Mode _mode;
-  d _alpha;          // Only valid in case of 'Leaking'
-  us _nfft;
-  vd _n_averages { 1, arma::fill::zeros}; // Only valid in case of 'Averaging'
-
-  PowerSpectra _ps;
-  /**
-   * @brief Current estimate of cross-spectral density
-   */
-  ccube _est;
-
-  /**
-   * @brief Buffer of storage of time data.
-   */
-  TimeBuffer _timeBuf;
-  /**
-   * @brief The amount of samples to keep in the overlap
-   */
-  us _overlap_keep;
-
-public:
-  /**
-   * @brief Initalize averaged power spectra computer. If a time constant is
-   * given > 0, it is used in a kind of exponential weighting.
-   *
-   * @param nfft The fft length
-   * @param w The window type.
-   * @param overlap_percentage A number 0 < overlap_percentage <= 100. It
-   * determines the amount of overlap used in Welch' method. A typical value is
-   * 50 %, i.e. 50.
-   * @param fs_tau Value should either be < 0, indicating that the
-   * estimate is averages over all time data.
-   * For a value = 0 the instantaneous power spectrum is returned, which can be
-   * interpreted as the spectrogram. For a value > 0 a exponential forgetting is
-   * used, where the value is used as the time constant such that the decay
-   * follows approximately the trend exp(-n/fs_tau), where n is the
-   * sample number in the power spectra. To choose 'fast' time weighting, set
-   * time_constant to the value of fs*0.125, where fs denotes the sampling
-   * frequency.
-   **/
-  AvSweepPowerSpectra(const us nfft = 2048,
-                 const Window::WindowType w = Window::WindowType::Hann,
-                 const d overlap_percentage = 50.,
-                 const d fs_tau = -1);
-                 
-  AvSweepPowerSpectra(const AvSweepPowerSpectra &) = delete;
-  AvSweepPowerSpectra &operator=(const AvSweepPowerSpectra &) = delete;
-
-  /**
-   * @brief Reset to empty state. Clears the time buffer and sets estimator to
-   * empty.
-   */
-  void reset();
-
-  /**
-   * @brief Compute an update of the power spectra based on given time data.
-   * Note that the number of channels is determined from the first time this
-   * function is called. If a later call has an incompatible number of
-   * channels, a runtime error is thrown.
-   *
-   * @param timedata
-   *
-   * @return a copy of the latest estimate of the power spectra. an
-   * update is only given if the amount of new time data is enough to compute a
-   * new estimate. if no new estimate is available, it returns an empty ccube.
-   * Note that the latest available estimate can be obtained using get_est().
-   * */
-  ccube compute(const dmat &timedata);
-
-  /**
-   * @brief Returns the latest estimate of cps (cross-power spectra.
-   *
-   * @return a copy of the latest estimate of the power spectra. an
-   * update is only given if the amount of new time data is enough to compute a
-   * new estimate. If no estimate is available, it returns an empty ccube.
-   * */
-  ccube get_est() const;
-
-  /**
-   * @brief The overlap is rounded to a certain amount of time samples. This
-   * function returns that value.
-   *
-   * @return The amount of samples in overlapping.
-   */
-  us exactOverlapSamples() const { return _ps.nfft - _overlap_keep; }
-};
-/** @} */

cpp_src/dsp/lasp_siggen.h

@@ -73,7 +73,7 @@ public:
    * @brief Mute the signal. Passes through the DC offset. No lock is hold. If
    * it just works one block later, than that is just the case.
    *
-   * @param mute if true
+   * @param mute if tre
    */
   void setMute(bool mute = true) { _muted = mute; _interruption_frame_count=0; }
 

cpp_src/dsp/lasp_siggen_impl.cpp

@@ -24,84 +24,60 @@ DEBUGTRACE_VARIABLES;
 
 Noise::Noise(){DEBUGTRACE_ENTER}
 
-vd Noise::genSignalUnscaled(us nframes)
-{
+vd Noise::genSignalUnscaled(us nframes) {
   return arma::randn<vd>(nframes);
 }
 void Noise::resetImpl() {}
 
 Sine::Sine(const d freq) : omg(2 * arma::datum::pi * freq) { DEBUGTRACE_ENTER; }
 
-vd Sine::genSignalUnscaled(const us nframes)
-{
+vd Sine::genSignalUnscaled(const us nframes) {
   /* DEBUGTRACE_ENTER; */
   slock lck(_mtx);
   const d pi = arma::datum::pi;
   vd phase_vec =
       arma::linspace(phase, phase + omg * (nframes - 1) / _fs, nframes);
   phase += omg * nframes / _fs;
-  while (phase > 2 * arma::datum::pi)
-  {
+  while (phase > 2 * arma::datum::pi) {
     phase -= 2 * pi;
   }
   return arma::sin(phase_vec);
 }
 
-vd Periodic::genSignalUnscaled(const us nframes)
-{
+vd Periodic::genSignalUnscaled(const us nframes) {
   vd res(nframes);
   slock lck(_mtx);
-  if (_signal.size() == 0)
-  {
+  if (_signal.size() == 0) {
     throw rte("No signal defined while calling");
   }
-  if (_signal.size() != A_.size())
-  {
-    std::cout << "Seq size: " << _signal.size() << ", A size: " << A_.size() << "\n";
-    throw rte("Sequence and amplitude envelopes have different lengths");
-  }
-  for (us i = 0; i < nframes; i++)
-  {
-    res(i) = A_[_cur_pos] * _signal[_cur_pos];
+  for (us i = 0; i < nframes; i++) {
+    res(i) = _signal[_cur_pos];
     _cur_pos++;
     _cur_pos %= _signal.size();
   }
   return res;
 }
 
-void Periodic::setA(const vd &A)
-{
-  A_ = A;
-}
-
 Sweep::Sweep(const d fl, const d fu, const d Ts, const d Tq, const us flags)
-    : fl_(fl), fu_(fu), Ts(Ts), Tq(Tq), flags(flags)
-{
-  if (fl <= 0 || fu < fl || Ts <= 0)
-  {
+    : fl_(fl), fu_(fu), Ts(Ts), Tq(Tq), flags(flags) {
+  if (fl <= 0 || fu < fl || Ts <= 0) {
     throw rte("Invalid sweep parameters");
   }
-  if ((flags & ForwardSweep) && (flags & BackwardSweep))
-  {
+  if ((flags & ForwardSweep) && (flags & BackwardSweep)) {
     throw rte(
         "Both forward and backward sweep flag set. Please only set either one "
         "or none for a continuous sweep");
   }
-  if ((flags & LinearSweep) && (flags & LogSweep))
-  {
+  if ((flags & LinearSweep) && (flags & LogSweep)) {
     throw rte(
         "Both logsweep and linear sweep flag set. Please only set either one.");
   }
-  if (!((flags & LinearSweep) || (flags & LogSweep)))
-  {
+  if (!((flags & LinearSweep) || (flags & LogSweep))) {
     throw rte("Either LinearSweep or LogSweep should be given as flag");
   }
-
-  resetImpl();
 }
 
-void Sweep::resetImpl()
-{
+void Sweep::resetImpl() {
   DEBUGTRACE_ENTER;
   slock lck(_mtx);
 
@@ -118,18 +94,14 @@ void Sweep::resetImpl()
   const us N = Ns + Nq;
 
   _signal = vd(N, arma::fill::zeros);
-  fn_ = vd(N, arma::fill::zeros);
   index = 0;
 
   d fl, fu;
   /* Swap fl and fu for a backward sweep */
-  if (backward_sweep)
-  {
+  if (backward_sweep) {
     fu = fl_;
     fl = fu_;
-  }
-  else
-  {
+  } else {
     /* Case of continuous sweep, or forward sweep */
     fl = fl_;
     fu = fu_;
@@ -138,10 +110,8 @@ void Sweep::resetImpl()
   d phase = 0;
 
   /* Linear sweep */
-  if (flags & LinearSweep)
-  {
-    if (forward_sweep || backward_sweep)
-    {
+  if (flags & LinearSweep) {
+    if (forward_sweep || backward_sweep) {
       /* Forward or backward sweep */
       /* TRACE(15, "Forward or backward sweep"); */
       us K = (us)(Dt * (fl * Ns + 0.5 * (Ns - 1) * (fu - fl)));
@@ -150,16 +120,12 @@ void Sweep::resetImpl()
       /* iVARTRACE(15, K); */
       /* dVARTRACE(15, eps); */
 
-      for (us n = 0; n < Ns; n++)
-      {
+      for (us n = 0; n < Ns; n++) {
         _signal(n) = d_sin(phase);
         d fn = fl + ((d)n) / Ns * (fu + eps - fl);
-        fn_(n) = fn;
         phase += 2 * arma::datum::pi * Dt * fn;
       }
-    }
-    else
-    {
+    } else {
       /* Continous sweep */
       /* TRACE(15, "continuous sweep"); */
 
@@ -184,24 +150,18 @@ void Sweep::resetImpl()
       /* dVARTRACE(15, eps); */
       d phase = 0;
 
-      for (us n = 0; n <= Ns; n++)
-      {
+      for (us n = 0; n <= Ns; n++) {
         /* iVARTRACE(17, n); */
-        if (n < N)
-        {
+        if (n < N) {
           _signal[n] = d_sin(phase);
         }
 
         d fn;
-        if (n <= Nf)
-        {
+        if (n <= Nf) {
           fn = fl + ((d)n) / Nf * (fu - fl);
-        }
-        else
-        {
+        } else {
           fn = fu - ((d)n - Nf) / Nb * (fu + eps - fl);
         }
-        fn_(n) = fn;
         /* dbgassert(fn >= 0, "BUG"); */
 
         phase += 2 * number_pi * Dt * fn;
@@ -209,12 +169,9 @@ void Sweep::resetImpl()
       /* This should be a very small number!! */
       /* dVARTRACE(15, phase); */
     }
-  }
-  else if (flags & LogSweep)
-  {
+  } else if (flags & LogSweep) {
     DEBUGTRACE_PRINT("Log sweep");
-    if (forward_sweep || backward_sweep)
-    {
+    if (forward_sweep || backward_sweep) {
       /* Forward or backward sweep */
       DEBUGTRACE_PRINT("Forward or backward sweep");
       d k1 = (fu / fl);
@@ -223,8 +180,7 @@ void Sweep::resetImpl()
 
       /* Iterate k to the right solution */
       d E;
-      for (us iter = 0; iter < 10; iter++)
-      {
+      for (us iter = 0; iter < 10; iter++) {
         E = 1 + K / (Dt * fl) * (d_pow(k, 1.0 / Ns) - 1) - k;
         d dEdk = K / (Dt * fl) * d_pow(k, 1.0 / Ns) / (Ns * k) - 1;
         k -= E / dEdk;
@@ -235,16 +191,12 @@ void Sweep::resetImpl()
       DEBUGTRACE_PRINT(k);
       DEBUGTRACE_PRINT(E);
 
-      for (us n = 0; n < Ns; n++)
-      {
+      for (us n = 0; n < Ns; n++) {
         _signal[n] = d_sin(phase);
         d fn = fl * d_pow(k, ((d)n) / Ns);
-        fn_(n) = fn;
         phase += 2 * number_pi * Dt * fn;
       }
-    }
-    else
-    {
+    } else {
       DEBUGTRACE_PRINT("Continuous sweep");
 
       const us Nf = Ns / 2;
@@ -260,8 +212,7 @@ void Sweep::resetImpl()
 
       /* Newton iterations to converge k to the value such that the sweep is
        * continuous */
-      for (us iter = 0; iter < NITER_NEWTON; iter++)
-      {
+      for (us iter = 0; iter < NITER_NEWTON; iter++) {
         E = (k - 1) / (d_pow(k, 1.0 / Nf) - 1) +
             (k - 1) / (1 - d_pow(k, -1.0 / Nb)) - K / Dt / fl;
         DEBUGTRACE_PRINT(E);
@@ -285,37 +236,29 @@ void Sweep::resetImpl()
       DEBUGTRACE_PRINT(k);
       DEBUGTRACE_PRINT(E);
 
-      for (us n = 0; n <= Ns; n++)
-      {
+      for (us n = 0; n <= Ns; n++) {
         /* iVARTRACE(17, n); */
-        if (n < Ns)
-        {
+        if (n < Ns) {
           _signal[n] = d_sin(phase);
         }
 
         d fn;
-        if (n <= Nf)
-        {
+        if (n <= Nf) {
           fn = fl * d_pow(k, ((d)n) / Nf);
-        }
-        else
-        {
+        } else {
           fn = fl * k * d_pow(1 / k, ((d)n - Nf) / Nb);
         }
-        fn_(n) = fn;
         /* dbgassert(fn >= 0, "BUG"); */
 
         phase += 2 * number_pi * Dt * fn;
-        while (phase > 2 * number_pi)
-          phase -= 2 * number_pi;
+        while (phase > 2 * number_pi) phase -= 2 * number_pi;
         /* dVARTRACE(17, phase); */
       }
       /* This should be a very small number!! */
       DEBUGTRACE_PRINT(phase);
     }
   }  // End of log sweep
-  else
-  {
+  else {
     // Either log or linear sweep had to be given as flags.
     assert(false);
   }

cpp_src/dsp/lasp_siggen_impl.h

@@ -58,7 +58,6 @@ class Sine : public Siggen {
  * periodic as the frequency can be any floating point value.
  */
 class Periodic: public Siggen {
-  
   protected:
     vd _signal { 1, arma::fill::zeros};
     us _cur_pos = 0;
@@ -69,14 +68,8 @@ class Periodic: public Siggen {
      * @return As stated above
      */
     vd getSequence() const { return _signal; }
-    vd A_ { 1, arma::fill::ones};
-  
-    void setA(const vd& A);
 
     virtual vd genSignalUnscaled(const us nframes) override final;
-
-    vd getA() const { return A_; }
-
     ~Periodic() = default;
 
 };
@@ -88,7 +81,6 @@ class Sweep : public Periodic {
   d fl_, fu_, Ts, Tq;
   us index;
   us flags;
-  vd fn_ { 1, arma::fill::zeros};  
 
   void resetImpl() override;
 
@@ -98,8 +90,6 @@ class Sweep : public Periodic {
   static constexpr int LinearSweep = 1 << 2;
   static constexpr int LogSweep = 1 << 3;
 
-  vd getfn() const { return fn_; }
-
   /**
    * Create a sweep signal
    *

cpp_src/pybind11/lasp_dsp_pybind.cpp

@@ -129,29 +129,6 @@ void init_dsp(py::module &m) {
     return CubeToNpy<c>(est);
   });
 
-  /// AvSweepPowerSpectra
-  py::class_<AvSweepPowerSpectra> asps(m, "AvSweepPowerSpectra");
-  asps.def(py::init<const us, const Window::WindowType, const d, const d>(),
-          py::arg("nfft") = 2048,
-          py::arg("windowType") = Window::WindowType::Hann,
-          py::arg("overlap_percentage") = 50.0,
-          py::arg("time_constant") = -1);
-
-  asps.def("compute", [](AvSweepPowerSpectra &asps, dpyarray timedata) {
-    std::optional<ccube> res;
-    dmat timedata_mat = NpyToMat<d, false>(timedata);
-    {
-      py::gil_scoped_release release;
-      res = asps.compute(timedata_mat);
-    }
-
-    return CubeToNpy<c>(res.value_or(ccube(0, 0, 0)));
-  });
-  asps.def("get_est", [](const AvSweepPowerSpectra &sps) {
-    ccube est = sps.get_est();
-    return CubeToNpy<c>(est);
-  });
-
   py::class_<SLM> slm(m, "cppSLM");
 
   slm.def_static("fromBiquads", [](const d fs, const d Lref, const us ds,

cpp_src/pybind11/lasp_siggen.cpp

@@ -43,22 +43,11 @@ void init_siggen(py::module &m) {
 
   py::class_<Periodic, std::shared_ptr<Periodic>> periodic(m, "Periodic",
                                                            siggen);
-  periodic.def("setA", 
-               [](Periodic &p, const dpyarray A) {
-                p.setA(NpyToCol<d, false>(A));
-                });
-
   periodic.def("getSequence",
                [](const Sweep &s) { return ColToNpy<d>(s.getSequence()); });
 
-  periodic.def("getA", 
-               [](const Sweep &s) { return ColToNpy<d>(s.getA()); });  
-  
-
   py::class_<Sweep, std::shared_ptr<Sweep>> sweep(m, "Sweep", periodic);
   sweep.def(py::init<const d, const d, const d, const d, const us>());
-  sweep.def("getfn", 
-               [](const Sweep &s) { return ColToNpy<d>(s.getfn()); });
   sweep.def_readonly_static("ForwardSweep", &Sweep::ForwardSweep);
   sweep.def_readonly_static("BackwardSweep", &Sweep::BackwardSweep);
   sweep.def_readonly_static("LinearSweep", &Sweep::LinearSweep);

python_src/lasp/__init__.py

@@ -10,11 +10,11 @@ from .lasp_cpp import *
 
 # from .lasp_imptube import *  # TwoMicImpedanceTube
 from .lasp_measurement import *  # Measurement, scaleBlockSens
-from .lasp_octavefilter import *
+from .lasp_octavefilter import *  # OverallFilterBank, SosOctaveFilterBank, SosThirdOctaveFilterBank
 from .lasp_slm import *  # SLM, Dummy
 from .lasp_record import *  # RecordStatus, Recording
-from .lasp_daqconfigs import *
-from .lasp_measurementset import *
+from .lasp_daqconfigs import *  # DaqConfigurations
+from .lasp_measurementset import *  # MeasurementSet
 
 # from .lasp_siggen import *  # SignalType, NoiseType, SiggenMessage, SiggenData, Siggen
 # from .lasp_weighcal import *  # WeighCal

python_src/lasp/filter/filterbank_design.py

@@ -15,7 +15,6 @@ import warnings
 import numpy as np
 # For designing second-order sections
 from scipy.signal import butter
-from ..lasp_config import LASP_NUMPY_FLOAT_TYPE
 
 from .fir_design import bandpass_fir_design
 from .fir_design import freqResponse as firFreqResponse
@@ -255,9 +254,9 @@ class FilterBankDesigner:
         fuu = self.fu(xu)
 
         if scale == 'lin':
-            freq = np.linspace(fll, fuu, npoints, dtype=LASP_NUMPY_FLOAT_TYPE)
+            freq = np.linspace(fll, fuu, npoints)
         elif scale == 'log':
-            freq = np.logspace(np.log10(fll), np.log10(fuu), npoints, dtype=LASP_NUMPY_FLOAT_TYPE)
+            freq = np.logspace(np.log10(fll), np.log10(fuu), npoints)
         else:
             raise ValueError(f'Invalid scale parameter: {scale}')
 

python_src/lasp/lasp_common.py

@@ -7,7 +7,6 @@ from dataclasses import dataclass
 from dataclasses_json import dataclass_json
 from enum import Enum, unique, auto
 from .lasp_cpp import DaqChannel
-from .lasp_config import LASP_NUMPY_FLOAT_TYPE
 
 """
 Common definitions used throughout the code.
@@ -396,9 +395,7 @@ def getTime(fs, N, start=0):
         start: Optional start ofset in number of samples
     """
     assert N > 0 and fs > 0
-    return np.linspace(
-        start, start + N / fs, N, endpoint=False, dtype=LASP_NUMPY_FLOAT_TYPE
-    )
+    return np.linspace(start, start + N/fs, N, endpoint=False)
 
 
 def getFreq(fs, nfft):
@@ -409,6 +406,6 @@ def getFreq(fs, nfft):
         fs: Sampling frequency [Hz]
         nfft: Fft length (int)
     """
-    df = fs / nfft  # frequency resolution
-    K = nfft // 2 + 1  # number of frequency bins
-    return np.linspace(0, (K - 1) * df, K, dtype=LASP_NUMPY_FLOAT_TYPE)
+    df = fs/nfft   # frequency resolution
+    K = nfft//2+1  # number of frequency bins
+    return np.linspace(0, (K-1)*df, K)

python_src/lasp/lasp_config.py

@@ -5,12 +5,9 @@ Author: J.A. de Jong - ASCEE
 
 Description: LASP configuration
 """
-
 import numpy as np
 from .lasp_cpp import LASP_DOUBLE_PRECISION
 
-__all__ = ["zeros", "ones", "empty", "LASP_NUMPY_FLOAT_TYPE", "LASP_NUMPY_COMPLEX_TYPE"]
-
 if LASP_DOUBLE_PRECISION:
     LASP_NUMPY_FLOAT_TYPE = np.float64
     LASP_NUMPY_COMPLEX_TYPE = np.complex128
@@ -19,28 +16,28 @@ else:
     LASP_NUMPY_COMPLEX_TYPE = np.float64
 
 
-def zeros(shape, dtype=float, order="F"):
-    if dtype is float:
+def zeros(shape, dtype=float, order='F'):
+    if dtype == float:
         return np.zeros(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order=order)
-    elif dtype is complex:
+    elif dtype == complex:
         return np.zeros(shape, dtype=LASP_NUMPY_COMPLEX_TYPE, order=order)
     else:
         raise RuntimeError(f"Unknown dtype: {dtype}")
 
 
-def ones(shape, dtype=float, order="F"):
-    if dtype is float:
+def ones(shape, dtype=float, order='F'):
+    if dtype == float:
         return np.ones(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order=order)
-    elif dtype is complex:
+    elif dtype == complex:
         return np.ones(shape, dtype=LASP_NUMPY_COMPLEX_TYPE, order=order)
     else:
         raise RuntimeError(f"Unknown dtype: {dtype}")
 
-
-def empty(shape, dtype=float, order="F"):
-    if dtype is float:
+def empty(shape, dtype=float, order='F'):
+    if dtype == float:
         return np.empty(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order=order)
-    elif dtype is complex:
+    elif dtype == complex:
         return np.empty(shape, dtype=LASP_NUMPY_COMPLEX_TYPE, order=order)
     else:
         raise RuntimeError(f"Unknown dtype: {dtype}")
+

python_src/lasp/lasp_measurement.py

@@ -15,7 +15,7 @@ from scipy.io import wavfile
 import os, time, wave, logging
 from .lasp_common import SIQtys, Qty, getFreq
 from .lasp_version import LASP_VERSION_MAJOR, LASP_VERSION_MINOR
-from .lasp_cpp import Window, DaqChannel, AvPowerSpectra, AvSweepPowerSpectra
+from .lasp_cpp import Window, DaqChannel, AvPowerSpectra
 from typing import List
 from functools import lru_cache
 from .lasp_config import ones
@@ -814,9 +814,7 @@ class Measurement:
             channels = list(range(self.nchannels))
 
         nchannels = len(channels)
-
-        # aps = AvPowerSpectra(nfft, window, overlap)
-        aps = AvSweepPowerSpectra(nfft, window, overlap)
+        aps = AvPowerSpectra(nfft, window, overlap)
         freq = getFreq(self.samplerate, nfft)
 
         for data in self.iterData(channels, **kwargs):
@@ -1007,8 +1005,7 @@ class Measurement:
 
         # Convert range to [-1, 1]
         # TODO: this is wrong for float data where full scale > 1
-        sensone = np.ones(data.shape[1])
-        
+        sensone = np.ones_like(self.sensitivity)
         data = scaleBlockSens(data, sensone)
 
         if dtype == "int16" or dtype == "int32":
@@ -1159,11 +1156,8 @@ class Measurement:
         if data.ndim != 2:
             data = data[:, np.newaxis]
 
-        try:
-            len(sensitivity)
-        except:
-            raise ValueError("Sensitivity should be given as array-like data type")
-        sensitivity = np.asarray(sensitivity)
+        if not (isinstance(sensitivity, np.ndarray) and sensitivity.ndim >= 1):
+            sensitivity = np.asarray(sensitivity)[np.newaxis]
 
         nchannels = data.shape[1]
         if nchannels != sensitivity.shape[0]:

python_src/lasp/lasp_octavefilter.py

@@ -11,7 +11,6 @@ __all__ = ["OverallFilterBank", "SosOctaveFilterBank", "SosThirdOctaveFilterBank
 
 from .filter.filterbank_design import OctaveBankDesigner, ThirdOctaveBankDesigner
 from .lasp_cpp import BiquadBank
-from .lasp_config import empty, LASP_NUMPY_FLOAT_TYPE
 import numpy as np
 
 
@@ -47,7 +46,7 @@ class OverallFilterBank:
         Ncur = data.shape[0]
         tend = tstart + Ncur / self.fs
 
-        t = np.linspace(tstart, tend, Ncur, endpoint=False, dtype=LASP_NUMPY_FLOAT_TYPE)
+        t = np.linspace(tstart, tend, Ncur, endpoint=False)
         self.N += Ncur
 
         output["Overall"] = {"t": t, "data": data, "x": 0}
@@ -115,7 +114,7 @@ class SosFilterBank:
         Ncur = data.shape[0]
         tend = tstart + Ncur / self.fs
 
-        t = np.linspace(tstart, tend, Ncur, endpoint=False, dtype=LASP_NUMPY_FLOAT_TYPE)
+        t = np.linspace(tstart, tend, Ncur, endpoint=False)
         self.N += Ncur
 
         for i, x in enumerate(self.xs):

python_src/lasp/lasp_slm.py

@@ -4,15 +4,14 @@
 Sound level meter implementation
 @author: J.A. de Jong - ASCEE
 """
-
 from .lasp_cpp import cppSLM
-from .lasp_config import empty, LASP_NUMPY_FLOAT_TYPE
+from .lasp_config import empty
 import numpy as np
-from .lasp_common import TimeWeighting, FreqWeighting, P_REF
+from .lasp_common import (TimeWeighting, FreqWeighting, P_REF)
 from .filter import SPLFilterDesigner
 import logging
 
-__all__ = ["SLM", "Dummy"]
+__all__ = ['SLM', 'Dummy']
 
 
 class Dummy:
@@ -90,24 +89,24 @@ class SLM:
         elif fw == FreqWeighting.Z:
             prefilter = None
         else:
-            raise ValueError(f"Not implemented prefilter {fw}")
+            raise ValueError(f'Not implemented prefilter {fw}')
 
         # 'Probe' size of filter coefficients
         self.nom_txt = []
 
         # This is a bit of a hack, as the 5 is hard-encoded here, but should in
         # fact be coming from somewhere else..
-        sos_overall = np.array([1, 0, 0, 1, 0, 0] * 5, dtype=float)
+        sos_overall = np.array([1, 0, 0, 1, 0, 0]*5, dtype=float)
 
         if fbdesigner is not None:
             assert fbdesigner.fs == fs
             sos_firstx = fbdesigner.createSOSFilter(self.xs[0]).flatten()
             self.nom_txt.append(fbdesigner.nominal_txt(self.xs[0]))
-            sos = empty((sos_firstx.size, nfilters), dtype=float, order="C")
+            sos = 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:

python_src/lasp/lasp_weighcal.py

@@ -7,7 +7,6 @@ Weighting and calibration filter in one
 from .lasp_common import FreqWeighting
 from .filter import SPLFilterDesigner
 from lasp.lasp_config import ones, empty
-from ..lasp_config import LASP_NUMPY_FLOAT_TYPE
 from .wrappers import FilterBank
 import numpy as np
 
@@ -40,7 +39,7 @@ class WeighCal:
         self.calfile = calfile
 
         # Frequencies used for the filter design
-        freq_design = np.linspace(0, 17e3, 3000, dtype=LASP_NUMPY_FLOAT_TYPE)
+        freq_design = np.linspace(0, 17e3, 3000)
         freq_design[-1] = fs/2
 
         # Objective function for the frequency response
@@ -141,7 +140,6 @@ class WeighCal:
         """
         Returns the frequency response of the designed FIR filter
         """
-        raise RuntimeError('This code bugs. TODO: Re-implement?')
         if freq is None:
             freq = np.logspace(1, np.log10(self.fs/2), 500)
         return (freq, frp(self.fs, freq, self._firs[chan]),