From 4cd390871aeeebc6458f26e4ee273892c8c0a41b Mon Sep 17 00:00:00 2001 From: Thijs Hekman Date: Tue, 3 Dec 2024 10:47:10 +0100 Subject: [PATCH] fully implemented AvSweepPowerSpectra --- cpp_src/dsp/lasp_avpowerspectra.cpp | 27 +++-- cpp_src/dsp/lasp_avpowerspectra.h | 11 +- cpp_src/pybind11/lasp_dsp_pybind.cpp | 154 ++++++++++++++++----------- python_src/lasp/lasp_measurement.py | 14 ++- 4 files changed, 131 insertions(+), 75 deletions(-) diff --git a/cpp_src/dsp/lasp_avpowerspectra.cpp b/cpp_src/dsp/lasp_avpowerspectra.cpp index 558a626..f78fb4d 100644 --- 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(_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(_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; +} diff --git a/cpp_src/dsp/lasp_avpowerspectra.h b/cpp_src/dsp/lasp_avpowerspectra.h index d972493..90b63a5 100644 --- 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. diff --git a/cpp_src/pybind11/lasp_dsp_pybind.cpp b/cpp_src/pybind11/lasp_dsp_pybind.cpp index d9bc9ca..b96e5eb 100644 --- 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(m, "Fft"); fft.def(py::init()); - fft.def("fft", [](Fft &f, dpyarray dat) { + fft.def("fft", [](Fft &f, dpyarray dat) + { if (dat.ndim() == 1) { return ColToNpy(f.fft(NpyToCol(dat))); } else if (dat.ndim() == 2) { return MatToNpy(f.fft(NpyToMat(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(f.ifft(NpyToCol(dat))); } else if (dat.ndim() == 2) { return MatToNpy(f.ifft(NpyToMat(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_> sbq(m, "SeriesBiquad", filter); - sbq.def(py::init([](dpyarray filter) { - return std::make_shared(NpyToCol(filter)); - })); - sbq.def("filter", [](SeriesBiquad &s, dpyarray input) { + sbq.def(py::init([](dpyarray filter) + { return std::make_shared(NpyToCol(filter)); })); + sbq.def("filter", [](SeriesBiquad &s, dpyarray input) + { vd res = NpyToCol(input); s.filter(res); - return ColToNpy(res); - }); + return ColToNpy(res); }); /// BiquadBank py::class_> bqb(m, "BiquadBank"); - bqb.def(py::init([](dpyarray coefs) { - return std::make_shared(NpyToMat(coefs)); - })); - bqb.def(py::init([](dpyarray coefs, dpyarray gains) { + bqb.def(py::init([](dpyarray coefs) + { return std::make_shared(NpyToMat(coefs)); })); + bqb.def(py::init([](dpyarray coefs, dpyarray gains) + { vd gains_arma = NpyToMat(gains); - return std::make_shared(NpyToMat(coefs), &gains_arma); - })); + return std::make_shared(NpyToMat(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(input); b.filter(inout); - return ColToNpy(inout); - }); + return ColToNpy(inout); }); /// PowerSpectra py::class_ ps(m, "PowerSpectra"); ps.def(py::init()); - ps.def("compute", [](PowerSpectra &ps, dpyarray input) { - return CubeToNpy(ps.compute(NpyToMat(input))); - }); + ps.def("compute", [](PowerSpectra &ps, dpyarray input) + { return CubeToNpy(ps.compute(NpyToMat(input))); }); /// AvPowerSpectra py::class_ 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 res; dmat timedata_mat = NpyToMat(timedata); { @@ -122,22 +122,36 @@ void init_dsp(py::module &m) { res = aps.compute(timedata_mat); } - return CubeToNpy(res.value_or(ccube(0, 0, 0))); - }); - aps.def("get_est", [](const AvPowerSpectra &ps) { + return CubeToNpy(res.value_or(ccube(0, 0, 0))); }); + aps.def("get_est", [](const AvPowerSpectra &ps) + { ccube est = ps.get_est(); - return CubeToNpy(est); - }); + return CubeToNpy(est); }); /// AvSweepPowerSpectra - py::class_ asps(m, "AvSweepPowerSpectra"); - asps.def(py::init(), - py::arg("nfft") = 2048, - py::arg("windowType") = Window::WindowType::Hann, - py::arg("overlap_percentage") = 50.0, - py::arg("time_constant") = -1); + py::class_> 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 theA = std::make_unique(NpyToCol(A)); - asps.def("compute", [](AvSweepPowerSpectra &asps, dpyarray timedata) { + return std::make_unique( + nfft, + windowtype, + overlap_percentage, + time_constant, + theA.get()); + + } + + )); + + + asps.def("compute", [](AvSweepPowerSpectra &asps, dpyarray timedata) + { std::optional res; dmat timedata_mat = NpyToMat(timedata); { @@ -145,44 +159,60 @@ void init_dsp(py::module &m) { res = asps.compute(timedata_mat); } - return CubeToNpy(res.value_or(ccube(0, 0, 0))); - }); - asps.def("get_est", [](const AvSweepPowerSpectra &sps) { + return CubeToNpy(res.value_or(ccube(0, 0, 0))); }); + asps.def("get_est", [](const AvSweepPowerSpectra &sps) + { ccube est = sps.get_est(); - return CubeToNpy(est); - }); + return CubeToNpy(est); }); + + asps.def("get_A", [](const AvSweepPowerSpectra &sps) + { + vd A = sps.get_A(); + return ColToNpy(A); }); py::class_ 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(bandpass)); - }); + const d tau, dpyarray bandpass) + { + return SLM::fromBiquads(fs, Lref, ds, tau, NpyToMat(bandpass)); }); slm.def_static("fromBiquads", [](const d fs, const d Lref, const us ds, const d tau, dpyarray prefilter, - py::array_t bandpass) { + py::array_t bandpass) + { return SLM::fromBiquads(fs, Lref, ds, tau, NpyToCol(prefilter), - NpyToMat(bandpass)); - }); + NpyToMat(bandpass)); }); - slm.def("run", [](SLM &slm, dpyarray in) { - return MatToNpy(slm.run(NpyToCol(in))); - }); - slm.def("Pm", [](const SLM &slm) { return ColToNpy(slm.Pm); }); - slm.def("Pmax", [](const SLM &slm) { return ColToNpy(slm.Pmax); }); - slm.def("Ppeak", [](const SLM &slm) { return ColToNpy(slm.Ppeak); }); + slm.def("run", [](SLM &slm, dpyarray in) + { + return MatToNpy(slm.run(NpyToCol(in))); }); + slm.def("Pm", [](const SLM &slm) + { + return ColToNpy(slm.Pm); }); + slm.def("Pmax", [](const SLM &slm) + { + return ColToNpy(slm.Pmax); }); + slm.def("Ppeak", [](const SLM &slm) + { + return ColToNpy(slm.Ppeak); }); - slm.def("Leq", [](const SLM &slm) { return ColToNpy(slm.Leq()); }); - slm.def("Lmax", [](const SLM &slm) { return ColToNpy(slm.Lmax()); }); - slm.def("Lpeak", [](const SLM &slm) { return ColToNpy(slm.Lpeak()); }); + slm.def("Leq", [](const SLM &slm) + { + return ColToNpy(slm.Leq()); }); + slm.def("Lmax", [](const SLM &slm) + { + return ColToNpy(slm.Lmax()); }); + slm.def("Lpeak", [](const SLM &slm) + { + return ColToNpy(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(freq); vd Xa = NpyToCol(X); - return ColToNpy(freqSmooth(freqa, Xa, w)); - }); + return ColToNpy(freqSmooth(freqa, Xa, w)); }); } /** @} */ diff --git a/python_src/lasp/lasp_measurement.py b/python_src/lasp/lasp_measurement.py index 4f7dff8..8192961 100644 --- 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)