diff --git a/cpp_src/dsp/lasp_freqsmooth.cpp b/cpp_src/dsp/lasp_freqsmooth.cpp index 8381771..29ef614 100644 --- a/cpp_src/dsp/lasp_freqsmooth.cpp +++ b/cpp_src/dsp/lasp_freqsmooth.cpp @@ -1,11 +1,14 @@ // #define DEBUGTRACE_ENABLED #include "lasp_freqsmooth.h" +#include + #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,10 +37,14 @@ vd freqSmooth(const vd& freq, const vd& X, const unsigned w) { d ac_pwr; if (firstFreqEqZero) { freq_min = freq(1); - ac_pwr = arma::sum(X.subvec(1, Nfreq - 1)); + if (power_correct) { + ac_pwr = arma::sum(X.subvec(1, Nfreq - 1)); + } } else { freq_min = freq(0); - ac_pwr = arma::sum(X); + if (power_correct) { + ac_pwr = arma::sum(X); + } } DEBUGTRACE_PRINT(freq_min); DEBUGTRACE_PRINT(freq_max); @@ -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); - const us idx_stop = std::min(k + mu, Nfreq_sm - 1); - // DEBUGTRACE_PRINT(idx_start) - // DEBUGTRACE_PRINT(idx_stop); + long idx_start = std::max(k - mu, 0l); + long idx_stop = std::min(k + mu, Nfreq_sm - 1); + + // 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 - d new_acpwr = arma::sum(Xsm.subvec(1, Nfreq - 1)); - Xsm.subvec(1, Nfreq - 1) *= ac_pwr / new_acpwr; + 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,8 +123,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 - d new_acpwr = arma::sum(Xsm); - Xsm *= ac_pwr / new_acpwr; + if (power_correct) { + d new_acpwr = arma::sum(Xsm); + Xsm *= ac_pwr / new_acpwr; + } } return Xsm; diff --git a/cpp_src/dsp/lasp_freqsmooth.h b/cpp_src/dsp/lasp_freqsmooth.h index 11220ee..3575ea4 100644 --- a/cpp_src/dsp/lasp_freqsmooth.h +++ b/cpp_src/dsp/lasp_freqsmooth.h @@ -1,21 +1,28 @@ #pragma once -#include #include -#include "lasp_types.h" -#include "lasp_mathtypes.h" +#include +#include "lasp_mathtypes.h" +#include "lasp_types.h" /** * \addtogroup dsp * @{ */ - /** - * @brief Frequency domain smoothing. Returns a smoothed variant of the power - * spectrum given as input. + * @brief Apply frequency domain smoothing to a Frequency domain (single + * 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); /** @} */