From 3ec15ec6459cf991c69a738da095ff9ea87378bd Mon Sep 17 00:00:00 2001
From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" <j.a.dejong@ascee.nl>
Date: Mon, 11 Mar 2024 16:04:24 +0100
Subject: [PATCH] New smoothing implementation, that runs a bit faster

---
 cpp_src/dsp/CMakeLists.txt           |  1 +
 cpp_src/pybind11/lasp_dsp_pybind.cpp | 15 ++++++++--
 python_src/lasp/tools/tools.py       | 42 +++++++++++++++++++++++++++-
 3 files changed, 54 insertions(+), 4 deletions(-)

diff --git a/cpp_src/dsp/CMakeLists.txt b/cpp_src/dsp/CMakeLists.txt
index c3332d6..ff2a219 100644
--- a/cpp_src/dsp/CMakeLists.txt
+++ b/cpp_src/dsp/CMakeLists.txt
@@ -16,6 +16,7 @@ set(lasp_dsp_files
     lasp_threadedindatahandler.cpp
     lasp_ppm.cpp
     lasp_clip.cpp
+    lasp_freqsmooth.cpp
     )
 
 
diff --git a/cpp_src/pybind11/lasp_dsp_pybind.cpp b/cpp_src/pybind11/lasp_dsp_pybind.cpp
index b5cedc3..91c9647 100644
--- a/cpp_src/pybind11/lasp_dsp_pybind.cpp
+++ b/cpp_src/pybind11/lasp_dsp_pybind.cpp
@@ -1,13 +1,16 @@
+#include <pybind11/pybind11.h>
+
+#include <iostream>
+
 #include "arma_npy.h"
 #include "lasp_avpowerspectra.h"
 #include "lasp_biquadbank.h"
 #include "lasp_fft.h"
 #include "lasp_filter.h"
+#include "lasp_freqsmoothing.h"
 #include "lasp_slm.h"
 #include "lasp_streammgr.h"
 #include "lasp_window.h"
-#include <iostream>
-#include <pybind11/pybind11.h>
 
 using std::cerr;
 using std::endl;
@@ -27,7 +30,6 @@ using rte = std::runtime_error;
  */
 
 void init_dsp(py::module &m) {
-
   py::class_<Fft> fft(m, "Fft");
   fft.def(py::init<us>());
   fft.def("fft", [](Fft &f, dpyarray dat) {
@@ -152,5 +154,12 @@ void init_dsp(py::module &m) {
   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) {
+    vd freqa = NpyToCol<d, false>(freq);
+    vd Xa = NpyToCol<d, false>(X);
+    return ColToNpy(freqSmooth(freqa, Xa, w));
+  });
 }
 /** @} */
diff --git a/python_src/lasp/tools/tools.py b/python_src/lasp/tools/tools.py
index 1454dc9..8caee6b 100644
--- a/python_src/lasp/tools/tools.py
+++ b/python_src/lasp/tools/tools.py
@@ -20,6 +20,7 @@ from enum import Enum, unique
 import copy
 import numpy as np
 from numpy import log2, pi, sin
+from ..lasp_cpp import freqSmooth
 
 
 @unique
@@ -152,7 +153,7 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth):
     return Q
 
 
-def smoothSpectralData(freq, M, sw: SmoothingWidth,
+def smoothSpectralData_old(freq, M, sw: SmoothingWidth,
                        st: SmoothingType = SmoothingType.levels):
     """
     Apply fractional octave smoothing to data in the frequency domain.
@@ -220,6 +221,45 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
 
     return Psm
 
+def smoothSpectralData(freq, M, sw: SmoothingWidth,
+                       st: SmoothingType = SmoothingType.levels):
+    """
+    Apply fractional octave smoothing to data in the frequency domain.
+
+    Args:
+        freq: array of frequencies of data points [Hz] - equally spaced
+        M: array of data, either power or dB
+        the smoothing type `st`, the smoothing is applied.
+        sw: smoothing width
+        st: smoothing type = data type of input data
+
+    Returns:
+        freq : array frequencies of data points [Hz]
+        Msm : float smoothed magnitude of data points
+    """
+    # Safety
+
+    if st == SmoothingType.ps:
+        assert np.min(M) >= 0, 'Power spectrum values cannot be negative'
+    if st == SmoothingType.levels and isinstance(M.dtype, complex):
+        raise RuntimeError('Decibel input should be real-valued')
+
+    # Convert to power
+    if st == SmoothingType.levels:
+        P = 10**(M/10)
+    elif st == SmoothingType.ps:
+        P = M
+    else:
+        raise RuntimeError(f"Incorrect SmoothingType: {st}")
+
+    Psm = freqSmooth(freq, P, sw.value[0])
+
+    # Convert to original format
+    if st == SmoothingType.levels:
+        Psm = 10*np.log10(Psm)
+
+    return Psm
+
 
 # %% Test
 if __name__ == "__main__":