added fractional octave smoothing to 'lasp/tools/tools.py'

lasp/tools/tools.py

@@ -0,0 +1,152 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May  6 14:49:03 2021
+
+@author: Casper
+
+Smooth data in the frequency domain
+"""
+
+# TO DO: check if everything is correct
+# TO DO: add possibility to insert data that is not lin spaced in frequency
+
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal.windows import gaussian
+
+
+# %% Smoothing function
+def oct_smooth(f, M, Noct, dB=False):
+    """
+    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
+
+    Parameters
+    ----------
+    f : float
+        frequencies of data points [Hz] - equally spaced
+    M : float
+        magnitude of data points [- or dB, specify in paramater 'dB']
+    Noct : int
+        smoothing strength: Noct=12 means 1/12 octave smoothing
+    dB : Bool
+        True if [M]=dB, False if [M]=absolute
+
+    Returns
+    -------
+    f : float
+        frequencies of data points [Hz]
+    Msm : float
+        smoothed magnitude of data points
+
+    """
+
+    # Settings
+    tr = 2  # truncate window after 2x std
+
+    # Safety
+    assert Noct > 0, '\'Noct\' must be absolute positive'
+    if Noct < 1: raise Warning('Check if \'Noct\' is entered correctly')
+    assert len(f)==len(M), 'f and M should have equal length'
+    if not dB: assert np.min(M) >= 0, 'absolute magnitude M cannot be negative'
+
+    # Initialize
+    L = len(M)                       # number of data points
+    P = 10**(M/10) if dB else  M**2  # convert magnitude --> power
+    Psm = np.zeros(L)                # smoothed power - to be calculated
+    x0 = 1 if f[0]==0 else 0         # skip first data point if zero frequency
+    df = f[1] - f[0]                 # frequency step
+
+    # Loop through data points
+    for x in range(x0, L):
+        # Find indices of data points to calculate current (smoothed) magnitude
+        fc = f[x]            # center freq. of smoothing window
+        Df = tr * fc / Noct  # freq. range of smoothing window
+        xl = int(np.ceil(x - 0.5*Df/df))  # desired lower index of frequency array to be used during smoothing
+        xu = int(np.floor(x + 0.5*Df/df)) + 1  # upper index + 1 (because half-open interval)
+
+        # 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
+            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
+
+    Msm = 10*np.log10(Psm) if dB else  Psm**0.5  # convert power --> magnitude
+
+    return Msm
+
+
+# %% Test
+if __name__ == "__main__":
+    # Initialize
+    Noct = 6  # 1/6 oct. smoothing
+
+    # Create dummy data
+    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 = 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
+
+
+    # Apply function
+    Msm = oct_smooth(f, M, Noct, dB)
+    fsm = f
+
+    # Plot
+
+    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])