Implemented alternative smoothing algorithm in tool

src/lasp/tools/tools.py

@@ -30,13 +30,14 @@ 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']
+__all__ = ['SmoothingType', 'smoothSpectralData', 'SmoothingWidth', 'smoothSpectralDataAlt', 'intHann']
 
 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):
@@ -298,6 +299,77 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
 #    return Psm
 
 
+def smoothSpectralDataAlt(freq, Mdb, sw: SmoothingWidth):
+    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
+    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:
+            return (sin(2*pi*x2) - sin(2*pi*x1))/(2*pi) + (x2-x1)
+
+
 # %% Test
 if __name__ == "__main__":
     """ Test function for evaluation and debugging