Alternative Smoothing: transformed to matrix method, minor bug fixes, merged test scripts

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
@@ -30,14 +27,13 @@ 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
 """
 
-__all__ = ['SmoothingType', 'smoothSpectralData', 'SmoothingWidth', 'smoothSpectralDataAlt', 'intHann']
+__all__ = ['SmoothingType', 'smoothSpectralData', 'SmoothingWidth']
 
 from enum import Enum, unique
 import bisect
 import copy
 import numpy as np
-from math import ceil, floor, sin
+
-from numpy import arange, log2, log10, pi
 
 @unique
 class SmoothingWidth(Enum):
@@ -95,14 +91,17 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth):
     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"
+    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
@@ -129,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
@@ -140,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
 
 
@@ -154,7 +157,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     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
@@ -227,63 +230,10 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
 
     return Psm
 
+# %% Alternative algorithm
+from numpy import arange, log2, log10, pi, ceil, floor, sin
 
-def smoothSpectralDataAlt(freq, Mdb, sw: SmoothingWidth):
+# Integrated Hann window
-    Noct = 1/sw.value[0]
-    # M = Mdb
-    M = 10**(Mdb/20)
-
-    f0 = 0
-    if freq[0] == 0:
-        f0 += 1
-
-    Nfreq = len(freq)
-
-    test_smoothed = np.array(M)
-
-    ifreq = freq/(freq[1]-freq[0])
-    ifreq = np.array(ifreq.astype(int))
-
-    ifreqMin = ifreq[f0]
-    ifreqMax = ifreq[Nfreq-1]
-
-    sfact = 2**(Noct/2)
-
-    maxNkp = ifreqMax - floor((ifreqMax-1)/sfact**2)+1
-    W = np.zeros(maxNkp)
-
-    kpmin = np.floor(ifreq/sfact).astype(int)
-    kpmax = np.ceil(ifreq*sfact).astype(int)
-
-    for ff in range(f0, len(M)):
-        if kpmin[ff] < ifreqMin:
-            kpmin[ff] = ifreqMin
-            kpmax[ff] = ceil(ifreq[ff]**2/ifreqMin)
-            NoctAct = log2(kpmax[ff]/kpmin[ff])
-        elif kpmax[ff] > ifreqMax:
-            kpmin[ff] = floor(ifreq[ff]**2/ifreqMax)
-            kpmax[ff] = ifreqMax
-            NoctAct = log2(kpmax[ff]/kpmin[ff])
-        else:
-            NoctAct = Noct
-
-        kp = arange(kpmin[ff], kpmax[ff]+1)
-
-        if NoctAct:
-            Phi1 = log2((kp - 0.5)/ifreq[ff])/NoctAct
-            Phi2 = log2((kp + 0.5)/ifreq[ff])/NoctAct
-
-            for ii in range(len(kp)):
-                W[ii] = intHann(Phi1[ii], Phi2[ii])
-
-            test_smoothed[ff] = np.dot(M[kpmin[ff]-ifreq[0]:kpmax[ff]-ifreq[0]+1], W[:ii+1])
-
-    test_smoothed = 20*log10(test_smoothed)
-
-    return test_smoothed
-
-
-# %% Integrated Hann window
 def intHann(x1, x2):
     if (x2 <= -1/2) or (x1 >= 1/2):
         return 0
@@ -299,6 +249,242 @@ def intHann(x1, x2):
             return (sin(2*pi*x2) - sin(2*pi*x1))/(2*pi) + (x2-x1)
 
 
+def smoothSpectralDataAlt(freq, MdB, sw: SmoothingWidth,
+                          st: SmoothingType = SmoothingType.levels):
+    """
+    According to Tylka_JAES_SmoothingWeights.pdf
+    "A Generalized Method for Fractional-Octave Smoothing of Transfer Functions
+    that Preserves Log-Frequency Symmetry"
+    https://doi.org/10.17743/jaes.2016.0053
+    par 1.3
+    eq. 16
+    """
+    Noct = 1/sw.value[0]
+    # M = MdB
+    M = 10**(MdB/20)
+
+    f0 = 0
+    if freq[0] == 0:
+        f0 += 1
+
+    Nfreq = len(freq)  # Number of frequenties
+    test_smoothed = np.array(M)  # Input [Power]
+
+    ifreq = freq/(freq[1]-freq[0])  # Frequency, normalized to step=1
+    ifreq = np.array(ifreq.astype(int))
+
+    ifreqMin = ifreq[f0]        # Min. freq, normalized to step=1
+    ifreqMax = ifreq[Nfreq-1]   # Max. freq, normalized to step=1
+
+    sfact = 2**(Noct/2)  # bounds are this factor from the center freq
+
+    maxNkp = ifreqMax - floor((ifreqMax-1)/sfact**2)+1
+    # W = np.zeros(int(np.round(maxNkp)))
+
+    kpmin = np.floor(ifreq/sfact).astype(int)  # min freq of window
+    kpmax = np.ceil(ifreq*sfact).astype(int)   # max freq of window
+
+    for ff in range(f0, len(M)):  # loop over input freq
+        if kpmin[ff] < ifreqMin:
+            kpmin[ff] = ifreqMin
+            kpmax[ff] = ceil(ifreq[ff]**2/ifreqMin)  # achieved Noct
+            if np.isclose(kpmin[ff], kpmax[ff]):
+                kpmax[ff] += 1
+            NoctAct = log2(kpmax[ff]/kpmin[ff])
+        elif kpmax[ff] > ifreqMax:
+            kpmin[ff] = floor(ifreq[ff]**2/ifreqMax)  # achieved Noct
+            kpmax[ff] = ifreqMax
+            if np.isclose(kpmin[ff], kpmax[ff]):
+                kpmin[ff] -= 1
+            NoctAct = log2(kpmax[ff]/kpmin[ff])
+        else:
+            NoctAct = Noct  # Noct = smoothing width (Noct=6 --> 1/6th octave)
+
+        kp = arange(kpmin[ff], kpmax[ff]+1)  # freqs of window
+
+        Phi1 = log2((kp - 0.5)/ifreq[ff])/NoctAct  # integration bounds for hann window
+        Phi2 = log2((kp + 0.5)/ifreq[ff])/NoctAct
+
+        W = np.zeros(len(kp))
+        for ii in range(len(kp)):
+            W[ii] = intHann(Phi1[ii], Phi2[ii])  # weight = integration of hann window between Phi1 and Phi2
+
+        test_smoothed[ff] = np.dot( M[kpmin[ff]-ifreq[0]:kpmax[ff]-ifreq[0]+1], W[:ii+1] )  # eq 16
+
+    test_smoothed = 20*log10(test_smoothed)
+
+    return test_smoothed
+
+
+def smoothCalcMatrixAlt(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)
+
+    According to Tylka_JAES_SmoothingWeights.pdf
+    "A Generalized Method for Fractional-Octave Smoothing of Transfer Functions
+    that Preserves Log-Frequency Symmetry"
+    https://doi.org/10.17743/jaes.2016.0053
+    par 1.3
+    eq. 16
+    """
+    # 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"
+    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
+
+    Noct /= 1.5  # empirical correction factor: window @ -6 dB at Noct bounds
+    ifreq = freq/(freq[1]-freq[0])  # frequency, normalized to step=1
+    ifreq = np.array(ifreq.astype(int))
+
+    ifreqMin = ifreq[x0]    # min. freq, normalized to step=1
+    ifreqMax = ifreq[L-1]   # max. freq, normalized to step=1
+
+    sfact = 2**((1/Noct)/2)  # bounds are this factor from the center freq
+
+    kpmin = np.floor(ifreq/sfact).astype(int)  # min freq of window
+    kpmax = np.ceil(ifreq*sfact).astype(int)   # max freq of window
+
+    for ff in range(x0, len(M)):  # loop over input freq
+        if kpmin[ff] < ifreqMin:
+            kpmin[ff] = ifreqMin
+            kpmax[ff] = ceil(ifreq[ff]**2/ifreqMin)  # decrease smooth. width
+            if np.isclose(kpmin[ff], kpmax[ff]):
+                kpmax[ff] += 1
+            NoctAct = 1/log2(kpmax[ff]/kpmin[ff])
+        elif kpmax[ff] > ifreqMax:
+            kpmin[ff] = floor(ifreq[ff]**2/ifreqMax)  # decrease smooth. width
+            kpmax[ff] = ifreqMax
+            if np.isclose(kpmin[ff], kpmax[ff]):
+                kpmin[ff] -= 1
+            NoctAct = 1/log2(kpmax[ff]/kpmin[ff])
+        else:
+            NoctAct = Noct
+
+        kp = arange(kpmin[ff], kpmax[ff]+1)  # freqs of window
+
+        # Integration bounds for Hann window
+        Phi1 = log2((kp - 0.5)/ifreq[ff]) * NoctAct
+        Phi2 = log2((kp + 0.5)/ifreq[ff]) * NoctAct
+
+        # Weights within window = integration of hann window between Phi1, Phi2
+        W = np.zeros(len(kp))
+        for ii in range(len(kp)):
+            W[ii] = intHann(Phi1[ii], Phi2[ii])
+
+        # Insert W at input freq ii, starting at index 'kpmin[ff]-ifreq[0]'
+        Q[ff, kpmin[ff]-ifreq[0]:kpmax[ff]-ifreq[0]+1] = W
+
+    # Normalize to conserve input power
+    Qpower = np.sum(Q, axis=0)
+    Q = Q / Qpower[np.newaxis, :]
+
+    return Q
+
+
+def smoothSpectralDataAltMatrix(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).
+    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
+        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.
+
+    # Safety
+    MM = copy.deepcopy(M)
+    Noct = sw.value[0]
+    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')
+    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
+
+    if st == SmoothingType.levels:
+        P = 10**(MM/10)  # magnitude [dB] --> power
+    else:
+        P = MM  # 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
+    if freq[0] == 0:
+        Psm[0] = P[0]              # Reuse old value in case first data..
+                                   # ..point is skipped. Not plotted any way.
+
+    # # Re-use smoothing matrix Q if available. Otherwise, calculate.
+    # # Store in dict 'Qdict'
+    # nfft = int(2*(len(freq)-1))
+    # key = f"nfft{nfft}_Noct{Noct}"  # matrix name
+
+    # if 'Qdict' not in globals():  # Guarantee Qdict exists
+    #     global Qdict
+    #     Qdict = {}
+
+    # if key in Qdict:
+    #     Q = Qdict[key]
+    # else:
+    #     Q = smoothCalcMatrixAlt(freq, sw)
+    #     Qdict[key] = Q
+
+    Q = smoothCalcMatrixAlt(freq, sw)
+
+    # Apply smoothing
+    Psm = np.matmul(Q, P)
+
+    if st == SmoothingType.levels:
+        Psm = 10*np.log10(Psm)
+
+    return Psm
+
 # %% Test
 if __name__ == "__main__":
     """ Test function for evaluation and debugging
@@ -307,53 +493,76 @@ if __name__ == "__main__":
     points. They should be treated and weighted differently.
     """
     import matplotlib.pyplot as plt
+    import time
+    plt.close('all')
+
     # Initialize
-    Noct = 2  # Noct = 6 for 1/6 oct. smoothing
+    Noct = 3  # Noct = 6 for 1/6 oct. smoothing
 
-    # Create dummy data set 1: noise
+    # # Create dummy data set 1: noise
-    fmin = 1e3   # [Hz] min freq
+    # fmin = 1e3   # [Hz] min freq
-    fmax = 24e3  # [Hz] max freq
+    # fmax = 24e3  # [Hz] max freq
-    Ndata = 200  # number of data points
+    # Ndata = 200  # number of data points
+    # freq = np.linspace(fmin, fmax, Ndata)  # frequency points
+    # # 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: single tone
+    fmin = 0      # [Hz] min freq
+    fmax = 5e3    # [Hz] max freq
+    Ndata = 2501  # number of data points
     freq = np.linspace(fmin, fmax, Ndata)  # frequency points
-    M = abs(0.4*np.random.normal(size=(Ndata,)))+0.01  #
+    M = 1e-4*np.random.normal(size=(Ndata,))
-    M = 20*np.log10(M)
+    M[500] = 1
+    MdB = 20*np.log10(abs(M))
 
-#    # Create dummy data set 2: dirac delta
-#    fmin = 3e3   # [Hz] min freq
-#    fmax = 24e3  # [Hz] max freq
-#    Ndata = 200  # number of data points
-#    freq = np.linspace(fmin, fmax, Ndata)  # frequency points
-#    M = 0 * abs(1+0.4*np.random.normal(size=(Ndata,))) + 0.01  #
-#    M[int(100)] = 1
-#    M = 20*np.log10(M)
-
-    # Apply function
     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)
+    if 'Qdict' in globals():
+        del Qdict
+
+    t0 = time.time()
+    Msm = smoothSpectralData(freq, MdB, sw, st)  # current algorithm
+    t1 = time.time()
+    MsmAlt = smoothSpectralDataAlt(freq, MdB, sw, st)  # alternative algorithm
+    t2 = time.time()
+    MsmAltMatrix = smoothSpectralDataAltMatrix(freq, MdB, sw, st)  # alternative algorithm, matrix method
+    t3 = time.time()
     fsm = freq
 
+    print(f"Smoothing time: {t1-t0} s")
+    print(f"Smoothing time: {t2-t1} s (Alt)")
+    print(f"Smoothing time: {t3-t2} s (Alt Matrix)")
 
     # Plot - lin frequency
     plt.figure()
-    plt.plot(freq, M, '.b')
+    plt.plot(freq, MdB, '.b')
     plt.plot(fsm, Msm,  'r')
+    plt.plot(fsm, MsmAlt,  'g')
+    plt.plot(fsm, MsmAltMatrix,  '--k')
     plt.xlabel('f (Hz)')
     plt.ylabel('magnitude')
-    plt.xlim([100, fmax])
+    plt.xlim((0, fmax))
+    plt.ylim((-90, 1))
+    plt.grid('both')
     plt.title('lin frequency')
-    plt.legend(['Raw', 'Smooth'])
+    plt.legend(['Raw', 'Smooth', 'SmoothAlt', 'SmoothAltMatrix'])
 
     # Plot - log frequency
     plt.figure()
-    plt.semilogx(freq, M, '.b')
+    plt.semilogx(freq, MdB, '.b')
     plt.semilogx(fsm, Msm, 'r')
+    plt.semilogx(fsm, MsmAlt, 'g')
+    plt.semilogx(fsm, MsmAltMatrix,  '--k')
     plt.xlabel('f (Hz)')
     plt.ylabel('magnitude')
-    plt.xlim([100, fmax])
+    plt.xlim((100, fmax))
+    plt.ylim((-90, 1))
+    plt.grid('both')
     plt.title('log frequency')
-    plt.legend(['Raw', 'Smooth 1'])
+    plt.legend(['Raw', 'Smooth', 'SmoothAlt', 'SmoothAltMatrix'])