Made power correction in smoothing algorithm optional. Window decreases in size symmetrically around the edged of the frequency spectrum

cpp_src/dsp/lasp_freqsmooth.cpp

@@ -1,11 +1,14 @@
 // #define DEBUGTRACE_ENABLED
 #include "lasp_freqsmooth.h"
 
+#include <cassert>
+
 #include "debugtrace.hpp"
 
 using rte = std::runtime_error;
 
-vd freqSmooth(const vd& freq, const vd& X, const unsigned w) {
+vd freqSmooth(const vd& freq, const vd& X, const unsigned w,
+              bool power_correct) {
   DEBUGTRACE_ENTER;
   if (freq.size() < 2) {
     throw rte("Invalid frequency size. Should be > 2");
@@ -34,11 +37,15 @@ vd freqSmooth(const vd& freq, const vd& X, const unsigned w) {
   d ac_pwr;
   if (firstFreqEqZero) {
     freq_min = freq(1);
+    if (power_correct) {
       ac_pwr = arma::sum(X.subvec(1, Nfreq - 1));
+    }
   } else {
     freq_min = freq(0);
+    if (power_correct) {
       ac_pwr = arma::sum(X);
     }
+  }
   DEBUGTRACE_PRINT(freq_min);
   DEBUGTRACE_PRINT(freq_max);
   const vd freq_log =
@@ -72,10 +79,20 @@ vd freqSmooth(const vd& freq, const vd& X, const unsigned w) {
   // Long is at least 32 bits. So +/- 2M points length
   for (long k = 0; k < Nfreq_sm; k++) {
     // const d fcur = freq_log(k);
-    const us idx_start = std::max(k - mu, 0l);
+    long idx_start = std::max(k - mu, 0l);
-    const us idx_stop = std::min(k + mu, Nfreq_sm - 1);
+    long idx_stop = std::min(k + mu, Nfreq_sm - 1);
-    // DEBUGTRACE_PRINT(idx_start)
+
-    // DEBUGTRACE_PRINT(idx_stop);
+    // Make window smaller at the sides (close to the end of the array)
+    if (idx_start == 0 || idx_stop == Nfreq_sm - 1) {
+      const long mu_edge = std::min(k - idx_start, idx_stop - k);
+      idx_start = k - mu_edge;
+      idx_stop = k + mu_edge;
+    }
+    assert(idx_stop < Nfreq_sm);
+    assert(idx_start < Nfreq_sm);
+
+    DEBUGTRACE_PRINT(idx_start)
+    DEBUGTRACE_PRINT(idx_stop);
 
     Xsm_log(k) = arma::mean(X_log.subvec(idx_start, idx_stop));
   }
@@ -95,8 +112,10 @@ vd freqSmooth(const vd& freq, const vd& X, const unsigned w) {
     Xsm(Nfreq - 1) = Xsm_log(Nfreq_sm - 1);
 
     // Final step: power-correct smoothed spectrum
+    if (power_correct) {
       d new_acpwr = arma::sum(Xsm.subvec(1, Nfreq - 1));
       Xsm.subvec(1, Nfreq - 1) *= ac_pwr / new_acpwr;
+    }
 
   } else {
     arma::interp1(freq_log, Xsm_log, freq, Xsm, "*linear");
@@ -104,9 +123,11 @@ vd freqSmooth(const vd& freq, const vd& X, const unsigned w) {
     Xsm(Nfreq - 1) = Xsm_log(Nfreq_sm - 1);
 
     // Final step: power-correct smoothed spectrum
+    if (power_correct) {
       d new_acpwr = arma::sum(Xsm);
       Xsm *= ac_pwr / new_acpwr;
     }
+  }
 
   return Xsm;
 }

cpp_src/dsp/lasp_freqsmooth.h

@@ -1,21 +1,28 @@
 #pragma once
-#include <vector>
 #include <memory>
-#include "lasp_types.h"
+#include <vector>
-#include "lasp_mathtypes.h"
 
+#include "lasp_mathtypes.h"
+#include "lasp_types.h"
 
 /**
  * \addtogroup dsp
  * @{
  */
 
-
 /**
- * @brief Frequency domain smoothing. Returns a smoothed variant of the power
+ * @brief Apply frequency domain smoothing to a Frequency domain (single
- * spectrum given as input.
+ * sided)signal power spectrum
  * 
+ * @param freq Frequency range
+ * @param X Signal pwr
+ * @param w Parameter determining the smoothing with. 1 = 1/1 octave, 3 = 1/3th
+ * octave and so on
+ * @param power_correct Apply a correction to the whole spectrum to make the
+ * signal power equal to the unsmoothed signal power.
+ * @return vd Smoothed spectrum
  */
-vd freqSmooth(const vd& freq,const vd& X,const unsigned w);
+vd freqSmooth(const vd& freq, const vd& X, const unsigned w,
+              bool power_correct = false);
 
 /** @} */