Smoothing: minor bug fix

src/lasp/tools/tools.py

@@ -5,13 +5,10 @@ Author: C. Jansen, J.A. de Jong - ASCEE V.O.F.
 
 Smooth data in the frequency domain.
 
-TODO: This function is rather slow as it
+TODO: This function is rather slow as it uses [for loops] in Python. Speed up.
-used Python for loops. The implementations should be speed up in the near
+NOTE: function requires lin frequency spaced input data
-future.
+TODO: accept input data that is not lin spaced in frequency
-TODO: check if the smoothing is correct: depending on whether the data points
+TODO: it makes more sense to output data that is log spaced in frequency
-are spaced lin of log in frequency, they should be given different weights.
-TODO: accept data that is not equally spaced in frequency
-TODO: output data that is log spaced in frequency
 
 Cutoff frequencies of window taken from
 http://www.huennebeck-online.de/software/download/src/index.html 15-10-2021
@@ -91,17 +88,20 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth):
     Warning: this method does not work on levels (dB)
     """
     # Settings
-    tr = 2  # truncate window after 2x std; shorter is faster and less accurate
+    tr = 3  # truncate window after 2x std; shorter is faster and less accurate
     Noct = sw.value[0]
     assert Noct > 0, "'Noct' must be absolute positive"
+    assert np.isclose(freq[-1]-freq[-2], freq[1]-freq[0]), "Input data must "\
+        "have a linear frequency spacing"
     if Noct < 1:
         raise Warning('Check if \'Noct\' is entered correctly')
 
     # Initialize
     L = len(freq)
     Q = np.zeros(shape=(L, L), dtype=np.float16)  # float16: keep size small
-
+    Q[0, 0] = 1  # in case first point is skipped
     x0 = 1 if freq[0] == 0 else 0  # Skip first data point if zero frequency
+
     # Loop over indices of raw frequency vector
     for x in range(x0, L):
         # Find indices of data points to calculate current (smoothed) magnitude
@@ -128,9 +128,9 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth):
 
         # Find indices corresponding to frequencies
         xl = bisect.bisect_left(freq, fl)  # index corresponding to fl
-        xu = bisect.bisect_left(freq, fu)
+        xu = bisect.bisect_right(freq, fu)
 
-        xl = xu-1 if xu-xl <= 0 else xl  # Guarantee window length of at least 1
+        xl = xu-1 if xu-xl <= 0 else xl  # Guarantee window length >= 1
 
         # Calculate window
         xg = np.arange(xl, xu)  # indices
@@ -139,6 +139,10 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth):
         gs /= np.sum(gs)        # normalize: integral=1
         Q[x, xl:xu] = gs        # add to matrix
 
+    # Normalize to conserve input power
+    Qpower = np.sum(Q, axis=0)
+    Q = Q / Qpower[np.newaxis, :]
+
     return Q
 
 
@@ -147,13 +151,13 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     """
     Apply fractional octave smoothing to magnitude data in frequency domain.
     Smoothing is performed to power, using a sliding Gaussian window with
-    variable length. The window is truncated after 2x std at either side.
+    variable length. The window is truncated after 3x std at either side.
 
     The implementation is not exact, because f is linearly spaced and
     fractional octave smoothing is related to log spaced data. In this
     implementation, the window extends with a fixed frequency step to either
     side. The deviation is largest when Noct is small (e.g. coarse smoothing).
-    Casper Jansen, 07-05-2021
+    07-05-2021
 
     Update 16-01-2023: speed up algorithm
         - Smoothing is performed using matrix multiplication
@@ -243,7 +247,8 @@ if __name__ == "__main__":
     fmax = 24e3  # [Hz] max freq
     Ndata = 200  # number of data points
     freq = np.linspace(fmin, fmax, Ndata)  # frequency points
-    M = abs(0.4*np.random.normal(size=(Ndata,)))+0.01  #
+    # freq = np.hstack((0, freq))
+    M = abs(0.4*np.random.normal(size=(len(freq),)))+0.01  #
     M = 20*np.log10(M)
 
 #    # Create dummy data set 2: dirac delta