Smoothing: vectorised + minor changes

src/lasp/tools/tools.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Author: C. Jansen, J.A. de Jong - ASCEE V.O.F.
+Author: T. Hekman, J.A. de Jong, C. Jansen - ASCEE V.O.F.
 
 Smooth data in the frequency domain.
 
@@ -174,6 +174,8 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
 
     """
     # TODO: Make this function multi-dimensional array aware.
+    # TODO: This does not work due to complex numbers. Should be split up in
+    # magnitude and phase.
 
     # Safety
     MM = copy.deepcopy(M)
@@ -196,8 +198,6 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
         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
 
@@ -235,86 +235,21 @@ from numpy import arange, log2, log10, pi, ceil, floor, sin
 
 # Integrated Hann window
 def intHann(x1, x2):
-    if (x2 <= -1/2) or (x1 >= 1/2):
-        return 0
-    elif x1 <= -1/2:
-        if x2 >= 1/2:
-            return 1
-        else:
-            return sin(2*pi*x2)/(2*pi) + x2 + 1/2
-    else:
-        if x2 >= 1/2:
-            return 1/2 - sin(2*pi*x1)/(2*pi) - x1
-        else:
+    """
+    Calculate integral of (part of) Hann window.
+    If the args are vectors, the return value will match those.
+
+    Args:
+        x1: lower bound [-0.5, 0.5]
+        x2: upper bound [-0.5, 0.5]
+    Return:
+        Integral of Hann window between x1 and x2
+    """
+    x1 = np.clip(x1, -0.5, 0.5)
+    x2 = np.clip(x2, -0.5, 0.5)
     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):
     """
@@ -336,7 +271,6 @@ def smoothCalcMatrixAlt(freq, sw: SmoothingWidth):
     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 "\
@@ -350,7 +284,8 @@ def smoothCalcMatrixAlt(freq, sw: SmoothingWidth):
     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
+    # Noct /= 1.5  # empirical correction factor: window @ -6 dB at Noct bounds
+    Noct /= 2    # empirical correction factor: window @ -3 dB at Noct bounds
     ifreq = freq/(freq[1]-freq[0])  # frequency, normalized to step=1
     ifreq = np.array(ifreq.astype(int))
 
@@ -362,7 +297,9 @@ def smoothCalcMatrixAlt(freq, sw: SmoothingWidth):
     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
+        # Find window bounds and actual smoothing width
         if kpmin[ff] < ifreqMin:
             kpmin[ff] = ifreqMin
             kpmax[ff] = ceil(ifreq[ff]**2/ifreqMin)  # decrease smooth. width
@@ -385,9 +322,7 @@ def smoothCalcMatrixAlt(freq, sw: SmoothingWidth):
         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])
+        W = intHann(Phi1, Phi2)
 
         # 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
@@ -497,7 +432,7 @@ if __name__ == "__main__":
     plt.close('all')
 
     # Initialize
-    Noct = 3  # Noct = 6 for 1/6 oct. smoothing
+    Noct = 1  # Noct = 6 for 1/6 oct. smoothing
 
     # # Create dummy data set 1: noise
     # fmin = 1e3   # [Hz] min freq
@@ -529,40 +464,35 @@ if __name__ == "__main__":
     t0 = time.time()
     Msm = smoothSpectralData(freq, MdB, sw, st)  # current algorithm
     t1 = time.time()
-    MsmAlt = smoothSpectralDataAlt(freq, MdB, sw, st)  # alternative algorithm
+    MsmAlt = smoothSpectralDataAltMatrix(freq, MdB, sw, st)  # alternative algorithm, matrix method
     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)")
+    print(f"Smoothing time: {t1-t0} s (Current)")
+    print(f"Smoothing time: {t2-t1} s (Alternative)")
 
     # Plot - lin frequency
     plt.figure()
     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((0, fmax))
     plt.ylim((-90, 1))
     plt.grid('both')
     plt.title('lin frequency')
-    plt.legend(['Raw', 'Smooth', 'SmoothAlt', 'SmoothAltMatrix'])
+    plt.legend(['Raw', 'Smooth', 'SmoothAlt'])
 
     # Plot - log frequency
     plt.figure()
     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.ylim((-90, 1))
     plt.grid('both')
     plt.title('log frequency')
-    plt.legend(['Raw', 'Smooth', 'SmoothAlt', 'SmoothAltMatrix'])
+    plt.legend(['Raw', 'Smooth', 'SmoothAlt'])