Updated octave smoothing function; completely different approach

lasp/tools/tools.py

@@ -3,19 +3,46 @@
 """
 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
+Smooth data in the frequency domain.
+
+TODO: This function is rather slow as it
 used Python for loops. The implementations should be speed up in the near
 future.
+TODO: check if the smoothing is correct: depending on whether the data points
+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
+math --> ReduceSpectrum.c --> ReduceSpectrum::smoothLogXScale()
+    fl = fcenter / sqrt(2^(1/Noct))  # lower cutoff
+    fu = fcenter * sqrt(2^(1/Noct))  # upper cutoff
+such that:
+    sqrt(fl * fu) = fcenter
+    fu = 2^(1/Noct) * fl
+
+Smoothing window taken from
+https://www.ap.com/technical-library/deriving-fractional-octave-spectra-from-
+    the-fft-with-apx/ 15-10-2021
+g = sqrt( 1/ ([1+[(f/fm - fm/f)*(1.507*b)]^6]) )
+where b = 3 for 1/3rd octave, f = frequency, fm = mid-band frequency, g = gain.
+Gain is related to magnitude; power is related to gain^2
 """
-from enum import Enum, unique
 
 __all__ = ['SmoothingType', 'smoothSpectralData', 'SmoothingWidth']
 
+from enum import Enum, unique
+import bisect
+import numpy as np
+
 
 @unique
 class SmoothingWidth(Enum):
     none = (0, 'No smoothing')
-    # three = (3, '1/3th octave smoothing')
+    one = (1, '1/1stoctave smoothing')
+    two = (2, '1/2th octave smoothing')
+    three = (3, '1/3rd octave smoothing')
     six = (6, '1/6th octave smoothing')
     twelve = (12, '1/12th octave smoothing')
     twfo = (24, '1/24th octave smoothing')
@@ -38,6 +65,7 @@ class SmoothingWidth(Enum):
     def getCurrent(cb):
         return list(SmoothingWidth)[cb.currentIndex()]
 
+
 class SmoothingType:
     levels = 'l', 'Levels'
     # tf = 'tf', 'Transfer function',
@@ -48,9 +76,6 @@ class SmoothingType:
 # TO DO: add possibility to insert data that is not lin spaced in frequency
 
 
-import numpy as np
-from scipy.signal.windows import gaussian
-
 def smoothSpectralData(freq, M, sw: SmoothingWidth,
                        st: SmoothingType = SmoothingType.levels):
     """
@@ -77,13 +102,14 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     # TODO: Make this function multi-dimensional array aware.
 
     # Settings
-    tr = 2  # truncate window after 2x std
+    tr = 2  # truncate window after 2x std; shorter is faster and less accurate
 
     # Safety
     Noct = sw.value[0]
     assert Noct > 0, "'Noct' must be absolute positive"
-    if Noct < 1: raise Warning('Check if \'Noct\' is entered correctly')
+    if Noct < 1:
-    assert len(freq)==len(M), 'f and M should have equal length'
+        raise Warning('Check if \'Noct\' is entered correctly')
+    assert len(freq) == len(M), 'f and M should have equal length'
 
     if st == SmoothingType.ps:
         assert np.min(M) >= 0, 'absolute magnitude M cannot be negative'
@@ -92,109 +118,218 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
 
     # Initialize
     L = M.shape[0]                   # number of data points
-    
+
-    P = M
     if st == SmoothingType.levels:
-        P = 10**(P/10)
+        P = 10**(M/10)  # magnitude [dB] --> power
+    else:
+        P = M  # data already given as power
     # TODO: This does not work due to complex numbers. Should be split up in
     # magnitude and phase.
     # elif st == SmoothingType.tf:
     #     P = P**2
 
     # P is power while smoothing. x are indices of P.
-    Psm = np.zeros_like(P)           # Smoothed power - to be calculated
+    Psm = np.zeros_like(P)         # Smoothed power - to be calculated
-    x0 = 1 if freq[0]==0 else 0      # Skip first data point if zero frequency
+    x0 = 1 if freq[0] == 0 else 0  # Skip first data point if zero frequency
-    Psm[0] = P[0]                    # Reuse old value in case first data..
+    Psm[0] = P[0]                  # Reuse old value in case first data..
-                                     # ..point is skipped. Not plotted any way.
+                                   # ..point is skipped. Not plotted any way.
-    df = freq[1] - freq[0]           # Frequency step
 
     # Loop through data points
     for x in range(x0, L):
         # Find indices of data points to calculate current (smoothed) magnitude
-        fc = freq[x]            # center freq. of smoothing window
-        Df = tr * fc / Noct  # freq. range of smoothing window
-
-        # desired lower index of frequency array to be used during smoothing
-        xl = int(np.ceil(x - 0.5*Df/df))  
-
-        # upper index + 1 (because half-open interval)
-        xu = int(np.floor(x + 0.5*Df/df)) + 1 
-
-        # Create window
-        Np = xu - xl  # number of points
-        std = Np / (2 * tr)
-        wind = gaussian(Np, std)  # Gaussian window
-
-        # Clip indices to valid range
         #
         # 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.
-        if xl < 0:
+        #
-            rl = 0 - xl       # remove this number of points at the lower end
+        # Implicitly, the window is truncated in this implementation, because
-            xl = xl + rl      # .. from f
+        # of the bisect search
-            wind = wind[rl:]  # .. and from window
+        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
+        xl = bisect.bisect_left(f, fl)  # index corresponding to fl
+        xu = bisect.bisect_left(f, fu)
 
-            # rl = 0 - xl      # remove this number of points at the lower end
+        # Calculate window
-            # xl = xl + rl     # .. from f
+        g = np.zeros(xu-xl)
-            # xu = xu - rl
+        for n, xi in enumerate(range(xl, xu)):
-            # wind = wind[rl:-rl]  # .. and from window
+            fi = f[xi]  # current frequency
-
+            g[n] = np.sqrt( 1/ ((1+((fi/fc - fc/fi)*(1.507*Noct))**6)) )
-        if xu > L:
+        g /= np.sum(g)  # normalize: integral=1
-            ru = xu - L       # remove this number of points at the upper end
-            xu = xu - ru
-            wind = wind[:-ru]
-
-            # ru = xu - L     # remove this number of points at the upper end
-            # xl = xl + ru
-            # xu = xu - ru
-            # wind = wind[ru:-ru]
 
         # Apply smoothing
-        wind_int = np.sum(wind)  # integral
+        Psm[x] = np.dot(g, P[xl:xu])
-        Psm[x] = np.dot(wind, P[xl:xu]) / wind_int  # apply window
 
     if st == SmoothingType.levels:
         Psm = 10*np.log10(Psm)
 
     return Psm
 
+## OLD VERSION
+#from scipy.signal.windows import gaussian
+#def smoothSpectralData(freq, M, sw: SmoothingWidth,
+#                       st: SmoothingType = SmoothingType.levels):
+#    """
+#    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.
+#
+#    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
+#
+#    Args:
+#        freq: array of frequencies of data points [Hz] - equally spaced
+#        M: array of either power, transfer functin or dB points. Depending on
+#        the smoothing type `st`, the smoothing is applied.
+#
+#    Returns:
+#        freq : array frequencies of data points [Hz]
+#        Msm : float smoothed magnitude of data points
+#
+#    """
+#    # TODO: Make this function multi-dimensional array aware.
+#
+#    # Settings
+#    tr = 2  # truncate window after 2x std; shorter is faster and less accurate
+#
+#    # Safety
+#    Noct = sw.value[0]
+#    assert Noct > 0, "'Noct' must be absolute positive"
+#    if Noct < 1: raise Warning('Check if \'Noct\' is entered correctly')
+#    assert len(freq)==len(M), 'f and M should have equal length'
+#
+#    if st == SmoothingType.ps:
+#        assert np.min(M) >= 0, 'absolute magnitude M cannot be negative'
+#    if st == SmoothingType.levels and isinstance(M.dtype, complex):
+#        raise RuntimeError('Decibel input should be real-valued')
+#
+#    # Initialize
+#    L = M.shape[0]                   # number of data points
+#    
+#    P = M
+#    if st == SmoothingType.levels:
+#        P = 10**(P/10)
+#    # TODO: This does not work due to complex numbers. Should be split up in
+#    # magnitude and phase.
+#    # elif st == SmoothingType.tf:
+#    #     P = P**2
+#
+#    # P is power while smoothing. x are indices of P.
+#    Psm = np.zeros_like(P)           # Smoothed power - to be calculated
+#    x0 = 1 if freq[0]==0 else 0      # Skip first data point if zero frequency
+#    Psm[0] = P[0]                    # Reuse old value in case first data..
+#                                     # ..point is skipped. Not plotted any way.
+#    df = freq[1] - freq[0]           # Frequency step
+#
+#    # Loop through data points
+#    for x in range(x0, L):
+#        # Find indices of data points to calculate current (smoothed) magnitude
+#        fc = freq[x]            # center freq. of smoothing window
+#        Df = tr * fc / Noct  # freq. range of smoothing window
+#
+#        # desired lower index of frequency array to be used during smoothing
+#        xl = int(np.ceil(x - 0.5*Df/df))  
+#
+#        # upper index + 1 (because half-open interval)
+#        xu = int(np.floor(x + 0.5*Df/df)) + 1 
+#
+#        # Create window, suitable for frequency lin-spaced data points 
+#        Np = xu - xl  # number of points
+#        std = Np / (2 * tr)
+#        wind = gaussian(Np, std)  # Gaussian window
+#
+#        # Clip indices to valid range
+#        #
+#        # 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.
+#        if xl < 0:
+#            rl = 0 - xl       # remove this number of points at the lower end
+#            xl = xl + rl      # .. from f
+#            wind = wind[rl:]  # .. and from window
+#
+#            # rl = 0 - xl      # remove this number of points at the lower end
+#            # xl = xl + rl     # .. from f
+#            # xu = xu - rl
+#            # wind = wind[rl:-rl]  # .. and from window
+#
+#        if xu > L:
+#            ru = xu - L       # remove this number of points at the upper end
+#            xu = xu - ru
+#            wind = wind[:-ru]
+#
+#            # ru = xu - L     # remove this number of points at the upper end
+#            # xl = xl + ru
+#            # xu = xu - ru
+#            # wind = wind[ru:-ru]
+#
+#        # Apply smoothing
+#        wind_int = np.sum(wind)  # integral
+#        Psm[x] = np.dot(wind, P[xl:xu]) / wind_int  # apply window
+#
+#    if st == SmoothingType.levels:
+#        Psm = 10*np.log10(Psm)
+#
+#    return Psm
+
 
 # %% Test
 if __name__ == "__main__":
+    """ Test function for evaluation and debugging
+
+    Note: make a distinction between lin and log spaced (in frequency) data
+    points. They should be treated and weighted differently.
+    """
     import matplotlib.pyplot as plt
     # Initialize
-    Noct = 6  # 1/6 oct. smoothing
+    Noct = 6  # Noct = 6 for 1/6 oct. smoothing
 
-    # Create dummy data
+    # Create dummy data set 1: noise
-    fmin = 3e3  # [Hz] min freq
+    fmin = 3e3   # [Hz] min freq
     fmax = 24e3  # [Hz] max freq
-    Ndata = 200 # number of data points
+    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  #
-    dB = False
     M = 20*np.log10(M)
-    dB = True
-
-
-    # M = f+1  # magnitude
-    # dB = False  # True if M is given in dB
 
+#    # Create dummy data set 2: dirac delta
+#    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 = 0 * abs(1+0.4*np.random.normal(size=(Ndata,))) + 0.01  #
+#    M[int(Ndata/5)] = 1
+#    M = 20*np.log10(M)
 
     # Apply function
-    Msm = oct_smooth(f, M, Noct, dB)
+    class sw:
+        value = [Noct]
+    st = SmoothingType.levels  # so data is given in dB
+
+    Msm = smoothSpectralData(f, M, sw, st)
     fsm = f
 
-    # Plot
+    # Plot - lin frequency
-
     plt.figure()
-    # plt.semilogx(f, M, '.b')
-    # plt.semilogx(fsm, Msm,  'r')
     plt.plot(f, M, '.b')
     plt.plot(fsm, Msm,  'r')
     plt.xlabel('f (Hz)')
     plt.ylabel('magnitude')
     plt.xlim([100, fmax])
+    plt.title('lin frequency')
+
+    # Plot - log frequency
+    plt.figure()
+    plt.semilogx(f, M, '.b')
+    plt.semilogx(fsm, Msm,  'r')
+    plt.xlabel('f (Hz)')
+    plt.ylabel('magnitude')
+    plt.xlim([100, fmax])
+    plt.title('log frequency')