From 82fc38a990c3a091501625b03141e2bbd8da90d5 Mon Sep 17 00:00:00 2001
From: Casper <c.jansen@ascee.nl>
Date: Mon, 18 Oct 2021 08:23:42 +0200
Subject: [PATCH] Fixed a bug in the new smoothing algorithm

---
 lasp/tools/tools.py | 33 +++++++++++++++++++++------------
 1 file changed, 21 insertions(+), 12 deletions(-)

diff --git a/lasp/tools/tools.py b/lasp/tools/tools.py
index 5f145a2..41b9d91 100644
--- a/lasp/tools/tools.py
+++ b/lasp/tools/tools.py
@@ -108,7 +108,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     # Safety
     MM = copy.deepcopy(M)
     Noct = sw.value[0]
-    assert M, "Smoothing function: input array is empty"
+    assert len(M) > 0, "Smoothing function: input array is empty"  # not sure if this works
     assert Noct > 0, "'Noct' must be absolute positive"
     if Noct < 1:
         raise Warning('Check if \'Noct\' is entered correctly')
@@ -141,17 +141,26 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     for x in range(x0, L):
         # Find indices of data points to calculate current (smoothed) magnitude
         #
-        # Indices beyond [0, L] point to non-existing data. This occurs when
-        # the smoothing windows nears the beginning or end of the series.
-        # Optional: if one end of the window is clipped, the other end
-        # could be clipped as well, to prevent an error on magnitude data with
-        # a slope. It however results in unsmoothed looking data at the ends.
-        #
-        # Implicitly, the window is truncated in this implementation, because
-        # of the bisect search
+        # Indices beyond [0, L] point to non-existing data. Beyond 0 does not
+        # occur in this implementation. Beyond L occurs when the smoothing
+        # window nears the end of the series.
+        # If one end of the window is truncated, the other end
+        # could be truncated as well, to prevent an error on magnitude data
+        # with a slope. It however results in unsmoothed looking data at the
+        # end.
         fc = freq[x]            # center freq. of smoothing window
         fl = fc / np.sqrt(2**(tr/Noct))  # lower cutoff
         fu = fc * np.sqrt(2**(tr/Noct))  # upper cutoff
+
+        # If the upper (frequency) side of the window is truncated because
+        # there is no data beyond the Nyquist frequency, also truncate the
+        # other side to keep it symmetric in a log(frequency) scale.
+        # So: fu / fc = fc / fl
+        fNq = freq[-1]
+        if fu > fNq:
+            fu = fNq  # no data beyond fNq
+            fl = fc**2 / fu  # keep window symmetric
+
         xl = bisect.bisect_left(f, fl)  # index corresponding to fl
         xu = bisect.bisect_left(f, fu)
 
@@ -292,14 +301,14 @@ if __name__ == "__main__":
     """
     import matplotlib.pyplot as plt
     # Initialize
-    Noct = 6  # Noct = 6 for 1/6 oct. smoothing
+    Noct = 2  # Noct = 6 for 1/6 oct. smoothing
 
     # Create dummy data set 1: noise
     fmin = 3e3   # [Hz] min freq
     fmax = 24e3  # [Hz] max freq
     Ndata = 200  # number of data points
     f = np.linspace(fmin, fmax, Ndata)  # frequency points
-    M = abs(1+0.4*np.random.normal(size=(Ndata,)))+0.01  #
+    M = abs(0.4*np.random.normal(size=(Ndata,)))+0.01  #
     M = 20*np.log10(M)
 
 #    # Create dummy data set 2: dirac delta
@@ -308,7 +317,7 @@ if __name__ == "__main__":
 #    Ndata = 200  # number of data points
 #    f = np.linspace(fmin, fmax, Ndata)  # frequency points
 #    M = 0 * abs(1+0.4*np.random.normal(size=(Ndata,))) + 0.01  #
-#    M[int(Ndata/5)] = 1
+#    M[int(100)] = 1
 #    M = 20*np.log10(M)
 
     # Apply function