fixed threshold at 1e-13

cpp_src/device/CMakeLists.txt

@@ -2,7 +2,7 @@
 include_directories(uldaq)
 include_directories(portaudio)
 
-add_library(lasp_device_lib OBJECT 
+add_library(lasp_device_lib OBJECT
     lasp_daq.cpp
     lasp_daqconfig.cpp
     lasp_daqdata.cpp
@@ -10,12 +10,8 @@ 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.
@@ -28,7 +24,9 @@ target_include_directories(lasp_device_lib INTERFACE
     ${CMAKE_CURRENT_SOURCE_DIR})
 
 if(LASP_HAS_ULDAQ)
-    target_link_libraries(lasp_device_lib uldaq)
+    add_subdirectory(uldaq)
+    target_include_directories(lasp_device_lib INTERFACE uldaq)
+    target_link_libraries(lasp_device_lib uldaq_backend uldaq)
 endif()
 if(LASP_HAS_RTAUDIO)
     target_link_libraries(lasp_device_lib rtaudio)

cpp_src/device/uldaq/CMakeLists.txt

@@ -0,0 +1,6 @@
+
+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_uldaq_common.h"
 #include "lasp_daq.h"
+#include "lasp_uldaq_common.h"
 
 string getErrMsg(UlError err) {
   string errstr;
@@ -21,11 +21,9 @@ 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,7 +10,8 @@ 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();
 
@@ -21,7 +22,8 @@ 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; */
@@ -34,11 +36,13 @@ 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)));
 
@@ -52,37 +56,46 @@ 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);
@@ -91,32 +104,42 @@ 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;
@@ -128,6 +151,136 @@ ccube AvPowerSpectra::compute(const dmat &timedata) {
   /// Othewise, we return an empty ccube.
   return run_once ? _est : ccube();
 }
-ccube AvPowerSpectra::get_est() const {
+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
+{
   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,3 +165,107 @@ 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 tre
+   * @param mute if true
    */
   void setMute(bool mute = true) { _muted = mute; _interruption_frame_count=0; }
 

cpp_src/dsp/lasp_siggen_impl.cpp

@@ -24,60 +24,84 @@ 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");
   }
-  for (us i = 0; i < nframes; i++) {
-    res(i) = _signal[_cur_pos];
+  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];
     _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);
 
@@ -86,7 +110,7 @@ void Sweep::resetImpl() {
   bool forward_sweep = flags & ForwardSweep;
   bool backward_sweep = flags & BackwardSweep;
 
-  const d Dt = 1 / _fs;  // Deltat
+  const d Dt = 1 / _fs; // Deltat
 
   // Estimate N, the number of samples in the sweep part (non-quiescent part):
   const us Ns = (us)(Ts * _fs);
@@ -94,14 +118,18 @@ 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_;
@@ -110,8 +138,10 @@ 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)));
@@ -120,12 +150,16 @@ 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"); */
 
@@ -150,18 +184,24 @@ 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;
@@ -169,9 +209,12 @@ 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);
@@ -180,7 +223,8 @@ 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;
@@ -191,12 +235,16 @@ 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;
@@ -212,7 +260,8 @@ 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);
@@ -236,29 +285,37 @@ 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 {
+  } // End of log sweep
+  else
+  {
     // Either log or linear sweep had to be given as flags.
     assert(false);
   }

cpp_src/dsp/lasp_siggen_impl.h

@@ -58,6 +58,7 @@ 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;
@@ -68,8 +69,14 @@ 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;
 
 };
@@ -81,7 +88,8 @@ class Sweep : public Periodic {
   d fl_, fu_, Ts, Tq;
   us index;
   us flags;
-
+  vd fn_ { 1, arma::fill::zeros};  
+  
   void resetImpl() override;
 
  public:
@@ -90,6 +98,8 @@ 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,6 +129,29 @@ 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,11 +43,22 @@ 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/filter/filterbank_design.py

@@ -15,6 +15,7 @@ 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
@@ -254,9 +255,9 @@ class FilterBankDesigner:
         fuu = self.fu(xu)
 
         if scale == 'lin':
-            freq = np.linspace(fll, fuu, npoints)
+            freq = np.linspace(fll, fuu, npoints, dtype=LASP_NUMPY_FLOAT_TYPE)
         elif scale == 'log':
-            freq = np.logspace(np.log10(fll), np.log10(fuu), npoints)
+            freq = np.logspace(np.log10(fll), np.log10(fuu), npoints, dtype=LASP_NUMPY_FLOAT_TYPE)
         else:
             raise ValueError(f'Invalid scale parameter: {scale}')
 

python_src/lasp/lasp_common.py

@@ -7,6 +7,7 @@ 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.
@@ -395,7 +396,9 @@ 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)
+    return np.linspace(
+        start, start + N / fs, N, endpoint=False, dtype=LASP_NUMPY_FLOAT_TYPE
+    )
 
 
 def getFreq(fs, nfft):
@@ -406,6 +409,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)
+    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)

python_src/lasp/lasp_config.py

@@ -5,9 +5,12 @@ 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
@@ -16,28 +19,28 @@ else:
     LASP_NUMPY_COMPLEX_TYPE = np.float64
 
 
-def zeros(shape, dtype=float, order='F'):
-    if dtype == float:
+def zeros(shape, dtype=float, order="F"):
+    if dtype is float:
         return np.zeros(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order=order)
-    elif dtype == complex:
+    elif dtype is 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 == float:
+def ones(shape, dtype=float, order="F"):
+    if dtype is float:
         return np.ones(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order=order)
-    elif dtype == complex:
+    elif dtype is 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 == float:
+
+def empty(shape, dtype=float, order="F"):
+    if dtype is float:
         return np.empty(shape, dtype=LASP_NUMPY_FLOAT_TYPE, order=order)
-    elif dtype == complex:
+    elif dtype is 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
+from .lasp_cpp import Window, DaqChannel, AvPowerSpectra, AvSweepPowerSpectra
 from typing import List
 from functools import lru_cache
 from .lasp_config import ones
@@ -814,9 +814,11 @@ class Measurement:
             channels = list(range(self.nchannels))
 
         nchannels = len(channels)
-        aps = AvPowerSpectra(nfft, window, overlap)
-        freq = getFreq(self.samplerate, nfft)
 
+        # aps = AvPowerSpectra(nfft, window, overlap)
+        aps = AvSweepPowerSpectra(nfft, window, overlap)
+        freq = getFreq(self.samplerate, nfft)
+        
         for data in self.iterData(channels, **kwargs):
             CS = aps.compute(data)
 
@@ -1005,7 +1007,8 @@ class Measurement:
 
         # Convert range to [-1, 1]
         # TODO: this is wrong for float data where full scale > 1
-        sensone = np.ones_like(self.sensitivity)
+        sensone = np.ones(data.shape[1])
+        
         data = scaleBlockSens(data, sensone)
 
         if dtype == "int16" or dtype == "int32":

python_src/lasp/lasp_octavefilter.py

@@ -11,6 +11,7 @@ __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
 
 
@@ -46,7 +47,7 @@ class OverallFilterBank:
         Ncur = data.shape[0]
         tend = tstart + Ncur / self.fs
 
-        t = np.linspace(tstart, tend, Ncur, endpoint=False)
+        t = np.linspace(tstart, tend, Ncur, endpoint=False, dtype=LASP_NUMPY_FLOAT_TYPE)
         self.N += Ncur
 
         output["Overall"] = {"t": t, "data": data, "x": 0}
@@ -114,7 +115,7 @@ class SosFilterBank:
         Ncur = data.shape[0]
         tend = tstart + Ncur / self.fs
 
-        t = np.linspace(tstart, tend, Ncur, endpoint=False)
+        t = np.linspace(tstart, tend, Ncur, endpoint=False, dtype=LASP_NUMPY_FLOAT_TYPE)
         self.N += Ncur
 
         for i, x in enumerate(self.xs):

python_src/lasp/lasp_slm.py

@@ -4,14 +4,15 @@
 Sound level meter implementation
 @author: J.A. de Jong - ASCEE
 """
+
 from .lasp_cpp import cppSLM
-from .lasp_config import empty
+from .lasp_config import empty, LASP_NUMPY_FLOAT_TYPE
 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:
@@ -89,24 +90,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,6 +7,7 @@ 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
 
@@ -39,7 +40,7 @@ class WeighCal:
         self.calfile = calfile
 
         # Frequencies used for the filter design
-        freq_design = np.linspace(0, 17e3, 3000)
+        freq_design = np.linspace(0, 17e3, 3000, dtype=LASP_NUMPY_FLOAT_TYPE)
         freq_design[-1] = fs/2
 
         # Objective function for the frequency response
@@ -140,6 +141,7 @@ 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]),