fully implemented AvSweepPowerSpectra

2024-12-03 10:47:10 +01:00 · 2024-12-03 10:47:10 +01:00 · 4cd390871a
commit 4cd390871a
parent 6cdc182b57
4 changed files with 131 additions and 75 deletions
--- a/cpp_src/dsp/lasp_avpowerspectra.cpp
+++ b/cpp_src/dsp/lasp_avpowerspectra.cpp
@ -158,7 +158,8 @@ ccube AvPowerSpectra::get_est() const

 AvSweepPowerSpectra::AvSweepPowerSpectra(const us nfft, const Window::WindowType w,
                                         const d overlap_percentage,
-                                         const d time_constant_times_fs)
+                                         const d time_constant_times_fs,
+                                         const vd *A)
    : _nfft(nfft), _ps(nfft, w)
 {

@ -167,6 +168,14 @@ AvSweepPowerSpectra::AvSweepPowerSpectra(const us nfft, const Window::WindowType
  {
    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);
@ -235,25 +244,24 @@ ccube AvSweepPowerSpectra::compute(const dmat &timedata)
        _n_averages = vd(nfreq, arma::fill::zeros);
      }
            
-      // d threshold = arma::mean(arma::real(temp.slice(0).col(0)));
-      d threshold = 1e-13;
+      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)) > threshold)
+        if (real(temp(f, 0, 0)) / _A(f) > 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];
          }
        }
      }
@ -284,3 +292,8 @@ ccube AvSweepPowerSpectra::get_est() const
 {
  return _est;
 }
+
+vd AvSweepPowerSpectra::get_A() const
+{
+  return _A;
+}
--- a/cpp_src/dsp/lasp_avpowerspectra.h
+++ b/cpp_src/dsp/lasp_avpowerspectra.h
@ -186,6 +186,7 @@ class AvSweepPowerSpectra {
  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;
  /**
@ -225,10 +226,12 @@ public:
  AvSweepPowerSpectra(const us nfft = 2048,
                 const Window::WindowType w = Window::WindowType::Hann,
                 const d overlap_percentage = 50.,
-                 const d fs_tau = -1);
+                 const d fs_tau = -1,
+                 const vd *A = nullptr); // amplitude envelope
                 
-  AvSweepPowerSpectra(const AvSweepPowerSpectra &) = delete;
-  AvSweepPowerSpectra &operator=(const AvSweepPowerSpectra &) = delete;
+  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
@ -260,6 +263,8 @@ public:
   * */
  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.
--- a/cpp_src/pybind11/lasp_dsp_pybind.cpp
+++ b/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,22 +122,36 @@ 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> 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);
+  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));

-  asps.def("compute", [](AvSweepPowerSpectra &asps, dpyarray timedata) {
+    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);
    {
@ -145,44 +159,60 @@ void init_dsp(py::module &m) {
      res = asps.compute(timedata_mat);
    }

-    return CubeToNpy<c>(res.value_or(ccube(0, 0, 0)));
-  });
-  asps.def("get_est", [](const AvSweepPowerSpectra &sps) {
+    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);
-  });
+    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)); });
 }
 /** @} */
--- a/python_src/lasp/lasp_measurement.py
+++ b/python_src/lasp/lasp_measurement.py
@ -19,6 +19,9 @@ 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

@ -842,10 +845,15 @@ class Measurement:

        nchannels = len(channels)

-        # aps = AvPowerSpectra(nfft, window, overlap)
-        aps = AvSweepPowerSpectra(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)