Made it such that the signal generator always starts from the start (pos 0)

INSTALLATION.md

@@ -1,145 +0,0 @@
-# Installation of LASP - Linux (Ubuntu-based)
-
-## Prerequisites
-
-Run the following on the command line to install all prerequisites on
-Debian-based Linux, x86-64:
-
-- `sudo apt install python3-pip libfftw3-3 libopenblas-base libusb-1.0-0 libpulse0`
-
-
-## Installation from wheel (recommended for non-developers)
-
-Go to: [LASP releases](https://code.ascee.nl/ASCEE/lasp/releases/latest/) and
-download the latest `.whl`. Then run:
-
-- `pip install lasp-*-linux_x86_64.whl`
-
-## From source (Ubuntu-based)
-
-### Prerequisites
-
-Run the following one-liner:
-
-- `sudo apt install -y git python3 python3-virtualenv python3-venv libopenblas-dev python3-pip libfftw3-dev libusb-1.0-0-dev libpulse-dev python3-build`
-
-If building RtAudio with the ALSA backend, you will also require the following packages:
-
-- `sudo apt install libclalsadrv-dev`
-
-If building RtAudio with the Jack Audio Connection Kit (JACK) backend, you will also require the following packages:
-
-- `sudo apt install libjack-jackd2-dev`
-
-### Download & build
-
-- `$ git clone --recursive https://code.ascee.nl/ASCEE/lasp.git`
-- `$ cd lasp`
-- `pip install -e .`
-
-# Building and installation for Raspberry Pi (Raspberry Pi OS)
-
-Run the following on the command line to install all prerequisites on
-Raspberry Pi OS:
-
-- `sudo apt install libfftw3-dev libopenblas64-dev libhdf5-dev libclalsadrv-dev`
-
-In a virtualenv: install `build`
-
-- `$ pip install build`
-
-Then run:
-
-- `$ git clone --recursive https://code.ascee.nl/ASCEE/lasp.git`
-- `$ cd lasp`
-- `$ pyproject-build`
-
-Which will generate a `whl` in the `dist` folder, that is redistributable for Raspberry Pis that run Raspberry Pi OS.
-
-When installing the `whl`, it appears that H5PY takes quite some time to install. To follow this process, run it it verbose mode.
-
-
-# Installation - (x86_64) Windows (with WinPython), build with MSYS2
-
-## Prerequisites
-
-- Download and install [WinPython](https://winpython.github.io)
-
-## From wheel
-
-- Download latest wheel from [LASP releases](https://code.ascee.nl/ASCEE/lasp/releases/latest/) and
-download the latest `.whl`. Then install with `pip`.
-
-## From source
-
-- Download and install [MSYS2](https://msys2.org). Make sure to install the
-  x86_64 version.
-- When unzipping WinPython, make sure to choose a proper and simple path, i.e.
-  C:\winpython
-- Download and install [Git for Windows](https://git-scm.com)
-- Open an MSYS2 **MINGW64** terminal, and install some tools we require:
-  - `$ pacman -S git`
-- Create a new virtualenv:
-  - `$ /c/winpython/<py-distr-dir>/python.exe -m venv venv`
-- Add the venv-python to the path (eases a lot of commands)
-  - `$ export PATH=$PATH:~/venv/Scripts`
-- Install `build`:
-  - `$ pip install build`
-- Clone LASP:
-  - `$ git clone --recurse-submodules https://code.ascee.nl/ascee/lasp && cd lasp`
-- If run for the first time, we have to install the libraries we depend on in
-  MSYS2 (this only has to be done on a fresh MSYS2 installation):
-  - `$ scripts/install_msys2_builddeps.sh`
-- Copy over required DLL's to be included in distribution:
-  - `scripts/copy_windows_dlls.sh`
-- And... build!
-  - `pyproject-build`
-- Lastly: the generated wheel can be installed in the current virtualenv:
-  - `pip install dist/lasp*.whl`
-
-
-# Documentation
-
-## Online
-
-[Online LASP documentation](https://lasp.ascee.nl/).
-
-## In directory (Linux/Debian)
-
-`$ sudo apt install doxygen graphviz`
-`$ pip install doxypypy`
-
-While still in lasp dir:
-
-`$ doxygen`
-
-This will build the documentation. It can be read by:
-
-`$ <YOUR-BROWSER> doc/html/index.html`
-
-# Usage
-
-- See examples directories for IPython notebooks.
-- Please refer to the [documentation](https://lasp.ascee.nl/) for features.
-
-
-# Development docs
-
-## Bumping version number
-
-When bumping the version number, please update the number in
-
-- `pyproject.toml`
-- `CMakeLists.txt`
-
-Then, create a commit with tag `vX.X.X`, and push it.
-
-## Updating to latest version (editable mode)
-
-When updating to the latest version of LASP in editable mode:
-
-```bash
-- $ git pull
-- $ git submodule update
-- $ pip install -e . -v
-```

README.md

@@ -1,11 +1,13 @@
 Library for Acoustic Signal Processing
 ======================================
 
+
 Welcome to LASP: Library for Acoustic Signal Processing. LASP is a C++ library
 with a Python interface which is supposed to acquire and process (multi) sensor data in real time on a PC and output results.
 
 Current features that are implemented:
 
+
 - Communication with data acquisition (DAQ) devices, of which:
     - Internal sound cards via the [RtAudio](http://www.music.mcgill.ca/~gary/rtaudio) backend. Many thanks to Gary P. Scavone et al.
     - [Measurement Computing](https://www.mccdaq.com) [DT9838A](https://www.mccdaq.com/Products/Sound-Vibration-DAQ/DT9837) signal analyzer.
@@ -29,8 +31,10 @@ Current features that are implemented:
 - Spectra data smoothing algorithms
 - Sensor calibration for microphones
 
+
 Future features (wish-list)
 
+
 - Conventional and delay-and-sum beam-forming algorithms
 - Impedance tube measurement processing
 
@@ -38,7 +42,151 @@ For now, the source code is well-documented on [lasp.ascee.nl](https://lasp.asce
 additional documentation (the math behind it). This is maintained
 in a sister repository [lasp-doc](https://code.ascee.nl/ascee/lasp-doc).
 
-For installation instructions, please visit
-[INSTALLATON.md](https://code.ascee.nl/ASCEE/lasp/src/branch/master/INSTALLATION.md).
-If you have any question(s), please feel free to contact us:
-[email](mailto:info@ascee.nl).
+If you have any question(s), please feel free to contact us: [email](info@ascee.nl).
+
+
+# Installation - Linux (Ubuntu-based)
+
+## Prerequisites
+
+Run the following on the command line to install all prerequisites on
+Debian-based Linux, x86-64:
+
+- `sudo apt install python3-pip libfftw3-3 libopenblas-base libusb-1.0-0 libpulse0`
+
+
+## Installation from wheel (recommended for non-developers)
+
+Go to: [LASP releases](https://code.ascee.nl/ASCEE/lasp/releases/latest/) and
+download the latest `.whl`. Then run:
+
+- `pip install lasp-*-linux_x86_64.whl`
+
+## From source (Ubuntu-based)
+
+### Prerequisites
+
+Run the following one-liner:
+
+- `sudo apt install -y git python3 python3-virtualenv python3-venv libopenblas-dev python3-pip libfftw3-dev libusb-1.0-0-dev libpulse-dev python3-build`
+
+If building RtAudio with the ALSA backend, you will also require the following packages:
+
+- `sudo apt install libclalsadrv-dev`
+
+If building RtAudio with the Jack Audio Connection Kit (JACK) backend, you will also require the following packages:
+
+- `sudo apt install libjack-jackd2-dev`
+
+### Download & build
+
+- `$ git clone --recursive https://code.ascee.nl/ASCEE/lasp.git`
+- `$ cd lasp`
+- `pip install -e .`
+
+# Building and installation for Raspberry Pi (Raspberry Pi OS)
+
+Run the following on the command line to install all prerequisites on
+Raspberry Pi OS:
+
+- `sudo apt install libfftw3-dev libopenblas64-dev libhdf5-dev libclalsadrv-dev`
+
+In a virtualenv: install `build`
+
+- `$ pip install build`
+
+Then run:
+
+- `$ git clone --recursive https://code.ascee.nl/ASCEE/lasp.git`
+- `$ cd lasp`
+- `$ pyproject-build`
+
+Which will generate a `whl` in the `dist` folder, that is redistributable for Raspberry Pis that run Raspberry Pi OS.
+
+When installing the `whl`, it appears that H5PY takes quite some time to install. To follow this process, run it it verbose mode.
+
+
+# Installation - (x86_64) Windows (with WinPython), build with MSYS2
+
+## Prerequisites
+
+- Download and install [WinPython](https://winpython.github.io)
+
+## From wheel
+
+- Download latest wheel from [LASP releases](https://code.ascee.nl/ASCEE/lasp/releases/latest/) and
+download the latest `.whl`. Then install with `pip`.
+
+## From source
+
+- Download and install [MSYS2](https://msys2.org). Make sure to install the
+  x86_64 version.
+- When unzipping WinPython, make sure to choose a proper and simple path, i.e.
+  C:\winpython
+- Download and install [Git for Windows](https://git-scm.com)
+- Open an MSYS2 **MINGW64** terminal, and install some tools we require:
+  - `$ pacman -S git`
+- Create a new virtualenv:
+  - `$ /c/winpython/<py-distr-dir>/python.exe -m venv venv`
+- Add the venv-python to the path (eases a lot of commands)
+  - `$ export PATH=$PATH:~/venv/Scripts`
+- Install `build`:
+  - `$ pip install build`
+- Clone LASP:
+  - `$ git clone --recurse-submodules https://code.ascee.nl/ascee/lasp && cd lasp`
+- If run for the first time, we have to install the libraries we depend on in
+  MSYS2 (this only has to be done on a fresh MSYS2 installation):
+  - `$ scripts/install_msys2_builddeps.sh`
+- Copy over required DLL's to be included in distribution:
+  - `scripts/copy_windows_dlls.sh`
+- And... build!
+  - `pyproject-build`
+- Lastly: the generated wheel can be installed in the current virtualenv:
+  - `pip install dist/lasp*.whl`
+
+
+# Documentation
+
+## Online
+
+[Online LASP documentation](https://lasp.ascee.nl/).
+
+## In directory (Linux/Debian)
+
+`$ sudo apt install doxygen graphviz`
+`$ pip install doxypypy`
+
+While still in lasp dir:
+
+`$ doxygen`
+
+This will build the documentation. It can be read by:
+
+`$ <YOUR-BROWSER> doc/html/index.html`
+
+# Usage
+
+- See examples directories for IPython notebooks.
+- Please refer to the [documentation](https://lasp.ascee.nl/) for features.
+
+
+# Development docs
+
+## Bumping version number
+
+When bumping the version number, please update the number in
+
+- `pyproject.toml`
+- `CMakeLists.txt`
+
+Then, create a commit with tag `vX.X.X`, and push it.
+
+## Updating to latest version (editable mode)
+
+When updating to the latest version of LASP in editable mode:
+
+```bash
+- $ git pull
+- $ git submodule update
+- $ pip install -e . -v
+```

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));
     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,149 @@ 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,
+                                         const vd *A)
+    : _nfft(nfft), _ps(nfft, w)
+{
+
+  DEBUGTRACE_ENTER;
+  if (overlap_percentage >= 100 || overlap_percentage < 0)
+  {
+    throw rte("Overlap percentage should be >= 0 and < 100");
+  }
+  if(A == nullptr) {
+    _A = vd(nfft, arma::fill::ones);
+  } else {
+    if(A->n_elem != (nfft / 2 + 1)) {
+      throw rte("Invalid length given for amplitude coefficients");
+    }
+    _A = arma::square(*A);
+  }
+
+  _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);
+      }
+            
+      vd corr_pow = arma::real(temp.slice(0).col(0)) / _A;
+
+      d threshold = pow(10.0, log10(arma::max(corr_pow)*arma::median(corr_pow))/2.0);
+      // d threshold = 1e-13;
+      
+      for (us f=0; f < nfreq; f++)
+      {
+        if (real(temp(f, 0, 0)) / _A(f) > threshold)
+        {
+          _n_averages[f]++;
+          if (_n_averages[f] == 1)
+          {
+            _est.row(f) = temp.row(f);
+          }
+          else
+          {
+            _est.row(f) = (static_cast<d>(_n_averages[f] - 1) / _n_averages[f]) * _est.row(f) +
+                                                             temp.row(f) / _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;
+}
+
+vd AvSweepPowerSpectra::get_A() const
+{
+  return _A;
+}

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,112 @@ 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'
+  vd _A; // squared amplitude envelope
+
+  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,
+                 const vd *A = nullptr); // amplitude envelope
+                 
+  AvSweepPowerSpectra(const AvSweepPowerSpectra &) = default;
+  AvSweepPowerSpectra &operator=(const AvSweepPowerSpectra &) = default;
+  AvSweepPowerSpectra (AvSweepPowerSpectra&& ) = default;
+
+  /**
+   * @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;
+
+  vd get_A() 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

@@ -38,6 +38,7 @@ protected:
   int _interruption_frame_count = 0;
 
   virtual void resetImpl() = 0;
+  virtual void resetPos() = 0;
   virtual vd genSignalUnscaled(const us nframes) = 0;
 
 public:
@@ -73,9 +74,9 @@ 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; }
+  void setMute(bool mute = true) { _muted = mute; _interruption_frame_count=0; resetPos();}
 
   /**
    * @brief Set the level of the signal generator

cpp_src/dsp/lasp_siggen_impl.cpp

@@ -24,60 +24,85 @@ DEBUGTRACE_VARIABLES;
 
 Noise::Noise(){DEBUGTRACE_ENTER}
 
-vd Noise::genSignalUnscaled(us nframes) {
+vd Noise::genSignalUnscaled(us nframes)
+{
   return arma::randn<vd>(nframes);
 }
 void Noise::resetImpl() {}
+void Noise::resetPos() {}
 
 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);
 
@@ -94,14 +119,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 +139,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 +151,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 +185,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 +210,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 +224,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 +236,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 +261,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 +286,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 {
+  else
+  {
     // Either log or linear sweep had to be given as flags.
     assert(false);
   }

cpp_src/dsp/lasp_siggen_impl.h

@@ -18,6 +18,7 @@ class Noise : public Siggen {
   d level_linear;
   virtual vd genSignalUnscaled(const us nframes) override;
   void resetImpl() override;
+  void resetPos() override;
 
  public:
   /**
@@ -39,6 +40,7 @@ class Sine : public Siggen {
 
  protected:
   void resetImpl() override final { phase = 0; }
+  void resetPos() override final { phase = 0; }
   virtual vd genSignalUnscaled(const us nframes) override final;
 
  public:
@@ -58,6 +60,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 +71,16 @@ class Periodic: public Siggen {
      * @return As stated above
      */
     vd getSequence() const { return _signal; }
+    vd A_ { 1, arma::fill::ones};
+  
+    void setA(const vd& A);
+
+    void resetPos() override final { _cur_pos = 0; }
 
     virtual vd genSignalUnscaled(const us nframes) override final;
+
+    vd getA() const { return A_; }
+
     ~Periodic() = default;
 
 };
@@ -81,6 +92,7 @@ class Sweep : public Periodic {
   d fl_, fu_, Ts, Tq;
   us index;
   us flags;
+  vd fn_ { 1, arma::fill::zeros};  
   
   void resetImpl() override;
   
@@ -90,6 +102,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

@@ -29,27 +29,28 @@ using rte = std::runtime_error;
  * @param m The Python module to add classes and methods to
  */
 
-void init_dsp(py::module &m) {
+void init_dsp(py::module &m)
+{
   py::class_<Fft> fft(m, "Fft");
   fft.def(py::init<us>());
-  fft.def("fft", [](Fft &f, dpyarray dat) {
+  fft.def("fft", [](Fft &f, dpyarray dat)
+          {
     if (dat.ndim() == 1) {
       return ColToNpy<c>(f.fft(NpyToCol<d, false>(dat)));
     } else if (dat.ndim() == 2) {
       return MatToNpy<c>(f.fft(NpyToMat<d, false>(dat)));
     } else {
       throw rte("Invalid dimensions of array");
-    }
-  });
-  fft.def("ifft", [](Fft &f, cpyarray dat) {
+    } });
+  fft.def("ifft", [](Fft &f, cpyarray dat)
+          {
     if (dat.ndim() == 1) {
       return ColToNpy<d>(f.ifft(NpyToCol<c, false>(dat)));
     } else if (dat.ndim() == 2) {
       return MatToNpy<d>(f.ifft(NpyToMat<c, false>(dat)));
     } else {
       throw rte("Invalid dimensions of array");
-    }
-  });
+    } });
 
   fft.def_static("load_fft_wisdom", &Fft::load_fft_wisdom);
   fft.def_static("store_fft_wisdom", &Fft::store_fft_wisdom);
@@ -71,41 +72,39 @@ void init_dsp(py::module &m) {
   /// SeriesBiquad
   py::class_<SeriesBiquad, std::shared_ptr<SeriesBiquad>> sbq(m, "SeriesBiquad",
                                                               filter);
-  sbq.def(py::init([](dpyarray filter) {
-    return std::make_shared<SeriesBiquad>(NpyToCol<d, false>(filter));
-  }));
-  sbq.def("filter", [](SeriesBiquad &s, dpyarray input) {
+  sbq.def(py::init([](dpyarray filter)
+                   { return std::make_shared<SeriesBiquad>(NpyToCol<d, false>(filter)); }));
+  sbq.def("filter", [](SeriesBiquad &s, dpyarray input)
+          {
     vd res = NpyToCol<d, true>(input);
     s.filter(res);
-    return ColToNpy<d>(res);
-  });
+    return ColToNpy<d>(res); });
 
   /// BiquadBank
   py::class_<BiquadBank, Filter, std::shared_ptr<BiquadBank>> bqb(m,
                                                                   "BiquadBank");
-  bqb.def(py::init([](dpyarray coefs) {
-    return std::make_shared<BiquadBank>(NpyToMat<d, false>(coefs));
-  }));
-  bqb.def(py::init([](dpyarray coefs, dpyarray gains) {
+  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);
-  }));
+    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) {
+          [](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);
-  });
+    return ColToNpy(inout); });
 
   /// PowerSpectra
   py::class_<PowerSpectra> ps(m, "PowerSpectra");
   ps.def(py::init<const us, const Window::WindowType>());
-  ps.def("compute", [](PowerSpectra &ps, dpyarray input) {
-    return CubeToNpy<c>(ps.compute(NpyToMat<d, false>(input)));
-  });
+  ps.def("compute", [](PowerSpectra &ps, dpyarray input)
+         { return CubeToNpy<c>(ps.compute(NpyToMat<d, false>(input))); });
 
   /// AvPowerSpectra
   py::class_<AvPowerSpectra> aps(m, "AvPowerSpectra");
@@ -114,7 +113,8 @@ void init_dsp(py::module &m) {
           py::arg("windowType") = Window::WindowType::Hann,
           py::arg("overlap_percentage") = 50.0, py::arg("time_constant") = -1);
 
-  aps.def("compute", [](AvPowerSpectra &aps, dpyarray timedata) {
+  aps.def("compute", [](AvPowerSpectra &aps, dpyarray timedata)
+          {
     std::optional<ccube> res;
     dmat timedata_mat = NpyToMat<d, false>(timedata);
     {
@@ -122,44 +122,97 @@ void init_dsp(py::module &m) {
       res = aps.compute(timedata_mat);
     }
 
-    return CubeToNpy<c>(res.value_or(ccube(0, 0, 0)));
-  });
-  aps.def("get_est", [](const AvPowerSpectra &ps) {
+    return CubeToNpy<c>(res.value_or(ccube(0, 0, 0))); });
+  aps.def("get_est", [](const AvPowerSpectra &ps)
+          {
     ccube est = ps.get_est();
-    return CubeToNpy<c>(est);
-  });
+    return CubeToNpy<c>(est); });
+
+  /// AvSweepPowerSpectra
+  py::class_<AvSweepPowerSpectra, std::unique_ptr<AvSweepPowerSpectra>> asps(m, "AvSweepPowerSpectra");
+  asps.def(py::init([](const us nfft, 
+                       const Window::WindowType windowtype, 
+                       const d overlap_percentage, 
+                       const d time_constant, 
+                       const dpyarray A)
+                       {
+    std::unique_ptr<vd> theA = std::make_unique<vd>(NpyToCol<d, false>(A));
+
+    return std::make_unique<AvSweepPowerSpectra>(
+        nfft,
+        windowtype,
+        overlap_percentage,
+        time_constant,
+        theA.get()); 
+
+                         }
+
+  ));
+
+
+  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); });
+
+  asps.def("get_A", [](const AvSweepPowerSpectra &sps)
+          {
+  vd A = sps.get_A();
+  return ColToNpy<d>(A); });
 
   py::class_<SLM> slm(m, "cppSLM");
 
   slm.def_static("fromBiquads", [](const d fs, const d Lref, const us ds,
-                                   const d tau, dpyarray bandpass) {
-    return SLM::fromBiquads(fs, Lref, ds, tau, NpyToMat<d, false>(bandpass));
-  });
+                                   const d tau, dpyarray bandpass)
+                 {
+    return SLM::fromBiquads(fs, Lref, ds, tau, NpyToMat<d, false>(bandpass)); });
 
   slm.def_static("fromBiquads", [](const d fs, const d Lref, const us ds,
                                    const d tau, dpyarray prefilter,
-                                   py::array_t<d> bandpass) {
+                                   py::array_t<d> bandpass)
+                 {
     return SLM::fromBiquads(fs, Lref, ds, tau, NpyToCol<d, false>(prefilter),
-                            NpyToMat<d, false>(bandpass));
-  });
+                            NpyToMat<d, false>(bandpass)); });
 
-  slm.def("run", [](SLM &slm, dpyarray in) {
-    return MatToNpy<d>(slm.run(NpyToCol<d, false>(in)));
-  });
-  slm.def("Pm", [](const SLM &slm) { return ColToNpy<d>(slm.Pm); });
-  slm.def("Pmax", [](const SLM &slm) { return ColToNpy<d>(slm.Pmax); });
-  slm.def("Ppeak", [](const SLM &slm) { return ColToNpy<d>(slm.Ppeak); });
+  slm.def("run", [](SLM &slm, dpyarray in)
+          {
+    return MatToNpy<d>(slm.run(NpyToCol<d, false>(in))); });
+  slm.def("Pm", [](const SLM &slm)
+          {
+    return ColToNpy<d>(slm.Pm); });
+  slm.def("Pmax", [](const SLM &slm)
+          {
+    return ColToNpy<d>(slm.Pmax); });
+  slm.def("Ppeak", [](const SLM &slm)
+          {
+    return ColToNpy<d>(slm.Ppeak); });
 
-  slm.def("Leq", [](const SLM &slm) { return ColToNpy<d>(slm.Leq()); });
-  slm.def("Lmax", [](const SLM &slm) { return ColToNpy<d>(slm.Lmax()); });
-  slm.def("Lpeak", [](const SLM &slm) { return ColToNpy<d>(slm.Lpeak()); });
+  slm.def("Leq", [](const SLM &slm)
+          {
+    return ColToNpy<d>(slm.Leq()); });
+  slm.def("Lmax", [](const SLM &slm)
+          {
+    return ColToNpy<d>(slm.Lmax()); });
+  slm.def("Lpeak", [](const SLM &slm)
+          {
+    return ColToNpy<d>(slm.Lpeak()); });
   slm.def_static("suggestedDownSamplingFac", &SLM::suggestedDownSamplingFac);
 
   // Frequency smoother
-  m.def("freqSmooth", [](dpyarray freq, dpyarray X, unsigned w) {
+  m.def("freqSmooth", [](dpyarray freq, dpyarray X, unsigned w)
+        {
     vd freqa = NpyToCol<d, false>(freq);
     vd Xa = NpyToCol<d, false>(X);
-    return ColToNpy(freqSmooth(freqa, Xa, w));
-  });
+    return ColToNpy(freqSmooth(freqa, Xa, w)); });
 }
 /** @} */

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/__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 *  # OverallFilterBank, SosOctaveFilterBank, SosThirdOctaveFilterBank
+from .lasp_octavefilter import *
 from .lasp_slm import *  # SLM, Dummy
 from .lasp_record import *  # RecordStatus, Recording
-from .lasp_daqconfigs import *  # DaqConfigurations
-from .lasp_measurementset import *  # MeasurementSet
+from .lasp_daqconfigs import *
+from .lasp_measurementset import *
 
 # from .lasp_siggen import *  # SignalType, NoiseType, SiggenMessage, SiggenData, Siggen
 # from .lasp_weighcal import *  # WeighCal

python_src/lasp/lasp_measurement.py

@@ -15,10 +15,13 @@ 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
+
+from scipy.interpolate import CubicSpline
+
 """!
 Author: J.A. de Jong - ASCEE
 
@@ -306,6 +309,16 @@ class Measurement:
             except KeyError:
                 self._type_int = 0
 
+            try:
+                dat = {}
+                
+                for key in f['siggen']['Data'].attrs:
+                    dat[key] = f['siggen']['Data'].attrs[key]
+
+                self._siggen = {'Type': f['siggen'].attrs["Type"], 'Data': dat}
+            except KeyError:
+                self._siggen = {'Type': 'Muted', 'Data': {}}
+
             # Due to a previous bug, the channel names were not stored
             # consistently, i.e. as 'channel_names' and later camelcase.
             try:
@@ -571,12 +584,29 @@ class Measurement:
         """
         self.setAttribute("type_int", type_.value)
 
+    def setSiggen(self, siggen_: dict):
+        """
+        Set the signal generator data
+        """
+        with self.file("r+") as f:
+            siggen_group = f.create_group("siggen", track_order=True)
+            siggen_group.attrs['Type'] = siggen_['Type']
+            data_group = siggen_group.create_group("Data", track_order=True)
+            for key in siggen_['Data'].keys():
+                data_group.attrs[key] = siggen_['Data'][key]
+
     def measurementType(self):
         """
         Returns type of measurement
         """
         return MeasurementType(self._type_int)
     
+    def signalGenerator(self):
+        """
+        Returns signal generator data
+        """
+        return self._siggen
+
     @property
     def name(self):
         """Returns filename base without extension."""
@@ -814,9 +844,16 @@ class Measurement:
             channels = list(range(self.nchannels))
 
         nchannels = len(channels)
-        aps = AvPowerSpectra(nfft, window, overlap)
+
         freq = getFreq(self.samplerate, nfft)
 
+        if self._siggen['Type'] == 'Sweep':
+            iA = CubicSpline(self._siggen['Data']['fA'], self._siggen['Data']['A'])
+
+            aps = AvSweepPowerSpectra(nfft, window, overlap, -1, iA(freq))
+        else:
+            aps = AvPowerSpectra(nfft, window, overlap)
+
         for data in self.iterData(channels, **kwargs):
             CS = aps.compute(data)
 
@@ -1005,7 +1042,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":
@@ -1156,8 +1194,11 @@ class Measurement:
         if data.ndim != 2:
             data = data[:, np.newaxis]
 
-        if not (isinstance(sensitivity, np.ndarray) and sensitivity.ndim >= 1):
-            sensitivity = np.asarray(sensitivity)[np.newaxis]
+        try:
+            len(sensitivity)
+        except:
+            raise ValueError("Sensitivity should be given as array-like data type")
+        sensitivity = np.asarray(sensitivity)
 
         nchannels = data.shape[1]
         if nchannels != sensitivity.shape[0]: