Merge remote-tracking branch 'origin/develop' into develop

src/lasp/tools/tools.py

@@ -78,6 +78,70 @@ class SmoothingType:
 # TO DO: add possibility to insert data that is not lin spaced in frequency
 
 
+def smoothCalcMatrix(freq, sw: SmoothingWidth):
+    """
+    Args:
+        freq: array of frequencies of data points [Hz] - equally spaced
+        sw: SmoothingWidth
+
+    Returns:
+        freq: array frequencies of data points [Hz]
+        Q: matrix to smooth power:  {fsm} = [Q] * {fraw}
+
+    Warning: this method does not work on levels (dB)
+    """
+    # Settings
+    tr = 2  # truncate window after 2x std; shorter is faster and less accurate
+    Noct = sw.value[0]
+    assert Noct > 0, "'Noct' must be absolute positive"
+    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
+
+    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
+        #
+        # 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
+
+        # Find indices corresponding to frequencies
+        xl = bisect.bisect_left(freq, fl)  # index corresponding to fl
+        xu = bisect.bisect_left(freq, fu)
+
+        xl = xu-1 if xu-xl <= 0 else xl  # Guarantee window length of at least 1
+
+        # Calculate window
+        xg = np.arange(xl, xu)  # indices
+        fg = freq[xg]           # [Hz] corresponding freq
+        gs = np.sqrt( 1/ ((1+((fg/fc - fc/fg)*(1.507*Noct))**6)) )  # raw windw
+        gs /= np.sum(gs)        # normalize: integral=1
+        Q[x, xl:xu] = gs        # add to matrix
+
+    return Q
+
+
 def smoothSpectralData(freq, M, sw: SmoothingWidth,
                        st: SmoothingType = SmoothingType.levels):
     """
@@ -91,6 +155,10 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     side. The deviation is largest when Noct is small (e.g. coarse smoothing).
     Casper Jansen, 07-05-2021
 
+    Update 16-01-2023: speed up algorithm
+        - Smoothing is performed using matrix multiplication
+        - The smoothing matrix is not calculated if it already exists
+
     Args:
         freq: array of frequencies of data points [Hz] - equally spaced
         M: array of either power, transfer functin or dB points. Depending on
@@ -103,9 +171,6 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     """
     # TODO: Make this function multi-dimensional array aware.
 
-    # Settings
-    tr = 2  # truncate window after 2x std; shorter is faster and less accurate
-
     # Safety
     MM = copy.deepcopy(M)
     Noct = sw.value[0]
@@ -134,169 +199,33 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
 
     # 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
+    if freq[0] == 0:
-    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.
 
-    # Loop through data points
+    # Re-use smoothing matrix Q if available. Otherwise, calculate.
-    for x in range(x0, L):
+    # Store in dict 'Qdict'
-        # Find indices of data points to calculate current (smoothed) magnitude
+    nfft = int(2*(len(freq)-1))
-        #
+    key = f"nfft{nfft}_Noct{Noct}"  # matrix name
-        # 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
+    if 'Qdict' not in globals():  # Guarantee Qdict exists
-        # there is no data beyond the Nyquist frequency, also truncate the
+        global Qdict
-        # other side to keep it symmetric in a log(frequency) scale.
+        Qdict = {}
-        # So: fu / fc = fc / fl
-        fNq = freq[-1]
-        if fu > fNq:
-            fu = fNq  # no data beyond fNq
-            fl = fc**2 / fu  # keep window symmetric
 
-        # Find indices corresponding to frequencies
+    if key in Qdict:
-        xl = bisect.bisect_left(freq, fl)  # index corresponding to fl
+        Q = Qdict[key]
-        xu = bisect.bisect_left(freq, fu)
+    else:
+        Q = smoothCalcMatrix(freq, sw)
+        Qdict[key] = Q
 
-        # Guarantee window length of at least 1
+    # Apply smoothing
-        if xu - xl <= 0:
+    Psm = np.matmul(Q, P)
-            xl = xu - 1
-
-        # Calculate window
-        g = np.zeros(xu-xl)
-        for n, xi in enumerate(range(xl, xu)):
-            fi = freq[xi]  # current frequency
-            g[n] = np.sqrt( 1/ ((1+((fi/fc - fc/fi)*(1.507*Noct))**6)) )
-        g /= np.sum(g)  # normalize: integral=1
-
-        # Apply smoothing
-        Psm[x] = np.dot(g, P[xl:xu])
 
     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__":
@@ -310,7 +239,7 @@ if __name__ == "__main__":
     Noct = 2  # Noct = 6 for 1/6 oct. smoothing
 
     # Create dummy data set 1: noise
-    fmin = 3e3   # [Hz] min freq
+    fmin = 1e3   # [Hz] min freq
     fmax = 24e3  # [Hz] max freq
     Ndata = 200  # number of data points
     freq = np.linspace(fmin, fmax, Ndata)  # frequency points
@@ -330,10 +259,13 @@ if __name__ == "__main__":
     class sw:
         value = [Noct]
     st = SmoothingType.levels  # so data is given in dB
+    # st = SmoothingType.ps  # so data is given in power
 
+    # Smooth
     Msm = smoothSpectralData(freq, M, sw, st)
     fsm = freq
 
+
     # Plot - lin frequency
     plt.figure()
     plt.plot(freq, M, '.b')
@@ -342,12 +274,14 @@ if __name__ == "__main__":
     plt.ylabel('magnitude')
     plt.xlim([100, fmax])
     plt.title('lin frequency')
+    plt.legend(['Raw', 'Smooth'])
 
     # Plot - log frequency
     plt.figure()
     plt.semilogx(freq, M, '.b')
-    plt.semilogx(fsm, Msm,  'r')
+    plt.semilogx(fsm, Msm, 'r')
     plt.xlabel('f (Hz)')
     plt.ylabel('magnitude')
     plt.xlim([100, fmax])
     plt.title('log frequency')
+    plt.legend(['Raw', 'Smooth 1'])