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

src/lasp/filter/biquad.py

@@ -23,11 +23,11 @@ y[n] = 1/ba[3] * (  ba[0] * x[n]   + ba[1] * x[n-1] + ba[2] * x[n-2] +
 
 """
 __all__ = ['peaking', 'biquadTF', 'notch', 'lowpass', 'highpass',
-           'highshelf', 'lowshelf']
+           'highshelf', 'lowshelf', 'LPcompensator', 'HPcompensator']
 
 from numpy import array, cos, pi, sin, sqrt
 from scipy.interpolate import interp1d
-from scipy.signal import sosfreqz
+from scipy.signal import sosfreqz, bilinear_zpk, zpk2sos
 
 
 def peaking(fs, f0, Q, gain):
@@ -157,6 +157,79 @@ def lowshelf(fs, f0, Q, gain):
     a2 = (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha
     return array([b0/a0, b1/a0, b2/a0, a0/a0, a1/a0, a2/a0])
 
+def LPcompensator(fs, f0o, Qo, f0n, Qn):
+    """
+    Shelving type filter that, when multiplied with a second-order low-pass 
+    filter, alters the response of that filter to a different second-order
+    low-pass filter.
+
+    Args:
+        fs: Sampling frequency [Hz]
+        f0o: Cut-off frequency of the original filter [Hz]
+        Qo: Quality factor of the original filter (~inverse of bandwidth)
+        f0n: Desired cut-off frequency [Hz]
+        Qn: Desired quality factor(~inverse of bandwidth)  
+    """
+    
+    omg0o = 2*pi*f0o
+    omg0n = 2*pi*f0n
+    
+    zRe = omg0o/(2*Qo)
+    zIm = omg0o*sqrt(1-1/(4*Qo**2))
+    z1 = -zRe + zIm*1j
+    z2 = -zRe - zIm*1j
+
+    pRe = omg0n/(2*Qn)
+    pIm = omg0n*sqrt(1-1/(4*Qn**2))
+    p1 = -pRe + pIm*1j
+    p2 = -pRe - pIm*1j
+
+    z= [z1, z2]
+    p = [p1, p2]
+    k = (pRe**2 + pIm**2)/(zRe**2 + zIm**2)
+    
+    zd, pd, kd = bilinear_zpk(z, p, k, fs)
+
+    sos = zpk2sos(zd,pd,kd)
+    
+    return sos[0]
+
+def HPcompensator(fs, f0o, Qo, f0n, Qn):
+    """
+    Shelving type filter that, when multiplied with a second-order high-pass
+    filter, alters the response of that filter to a different second-order 
+    high-pass filter.
+
+    Args:
+        fs: Sampling frequency [Hz]
+        f0o: Cut-on frequency of the original filter [Hz]
+        Qo: Quality factor of the original filter (~inverse of bandwidth)
+        f0n: Desired cut-on frequency [Hz]
+        Qn: Desired quality factor(~inverse of bandwidth)  
+    """
+    
+    omg0o = 2*pi*f0o
+    omg0n = 2*pi*f0n
+    
+    zRe = omg0o/(2*Qo)
+    zIm = omg0o*sqrt(1-1/(4*Qo**2))
+    z1 = -zRe + zIm*1j
+    z2 = -zRe - zIm*1j
+    
+    pRe = omg0n/(2*Qn)
+    pIm = omg0n*sqrt(1-1/(4*Qn**2))
+    p1 = -pRe + pIm*1j
+    p2 = -pRe - pIm*1j
+
+    z= [z1, z2]
+    p = [p1, p2]
+    k = 1
+    
+    zd, pd, kd = bilinear_zpk(z, p, k, fs)
+
+    sos = zpk2sos(zd,pd,kd)
+    
+    return sos[0]
 
 def biquadTF(fs, freq, sos):
     """

src/lasp/tools/tools.py

@@ -1,41 +1,25 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Author: C. Jansen, J.A. de Jong - ASCEE V.O.F.
+Author: T. Hekman, 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: Make SmoothSpectralData() multi-dimensional array aware.
-TODO: accept data that is not equally spaced in frequency
+TODO: Smoothing does not work due to complex numbers. Is it reasonable to
-TODO: output data that is log spaced in frequency
+smooth complex data? If so, the data could be split in magnitude and phase.
-
-Cutoff frequencies of window taken from
-http://www.huennebeck-online.de/software/download/src/index.html 15-10-2021
-math --> ReduceSpectrum.c --> ReduceSpectrum::smoothLogXScale()
-    fl = fcenter / sqrt(2^(1/Noct))  # lower cutoff
-    fu = fcenter * sqrt(2^(1/Noct))  # upper cutoff
-such that:
-    sqrt(fl * fu) = fcenter
-    fu = 2^(1/Noct) * fl
-
-Smoothing window taken from
-https://www.ap.com/technical-library/deriving-fractional-octave-spectra-from-
-    the-fft-with-apx/ 15-10-2021
-g = sqrt( 1/ ([1+[(f/fm - fm/f)*(1.507*b)]^6]) )
-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']
 
 from enum import Enum, unique
-import bisect
 import copy
 import numpy as np
+from numpy import log2, pi, sin
 
 
 @unique
@@ -76,6 +60,22 @@ class SmoothingType:
 
 # TO DO: check if everything is correct
 # TO DO: add possibility to insert data that is not lin spaced in frequency
+# Integrated Hann window
+
+def intHann(x1, x2):
+    """
+    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 smoothCalcMatrix(freq, sw: SmoothingWidth):
@@ -86,58 +86,68 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth):
 
     Returns:
         freq: array frequencies of data points [Hz]
-        Q: matrix to smooth power:  {fsm} = [Q] * {fraw}
+        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
 
-    x0 = 1 if freq[0] == 0 else 0  # Skip first data point if zero frequency
+    Noct /= 2  # corr. factor: window @ -3dB at Noct bounds (Noct/1.5 for -6dB)
-    # Loop over indices of raw frequency vector
+    ifreq = freq/(freq[1]-freq[0])  # frequency, normalized to step=1
-    for x in range(x0, L):
+    ifreq = np.array(ifreq.astype(int))
-        # Find indices of data points to calculate current (smoothed) magnitude
+    ifreqMin = ifreq[x0]     # min. freq, normalized to step=1
-        #
+    ifreqMax = ifreq[L-1]    # max. freq, normalized to step=1
-        # Indices beyond [0, L] point to non-existing data. Beyond 0 does not
+    sfact = 2**((1/Noct)/2)  # bounds are this factor from the center freq
-        # 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
+    kpmin = np.floor(ifreq/sfact).astype(int)  # min freq of window
-        # there is no data beyond the Nyquist frequency, also truncate the
+    kpmax = np.ceil(ifreq*sfact).astype(int)   # max freq of window
-        # 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
+    for ff in range(x0, L):  # loop over input freq
-        xl = bisect.bisect_left(freq, fl)  # index corresponding to fl
+        # Find window bounds and actual smoothing width
-        xu = bisect.bisect_left(freq, fu)
+        # Keep window symmetrical if one side is truncated
+        if kpmin[ff] < ifreqMin:
+            kpmin[ff] = ifreqMin
+            kpmax[ff] = np.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] = np.floor(ifreq[ff]**2/ifreqMax)  # decrease smoothwidth
+            kpmax[ff] = ifreqMax
+            if np.isclose(kpmin[ff], kpmax[ff]):
+                kpmin[ff] -= 1
+            NoctAct = 1/log2(kpmax[ff]/kpmin[ff])
+        else:
+            NoctAct = Noct
 
-        xl = xu-1 if xu-xl <= 0 else xl  # Guarantee window length of at least 1
+        kp = np.arange(kpmin[ff], kpmax[ff]+1)       # freqs of window
+        Phi1 = log2((kp - 0.5)/ifreq[ff]) * NoctAct  # integr. bounds Hann wind
+        Phi2 = log2((kp + 0.5)/ifreq[ff]) * NoctAct
+        W = intHann(Phi1, Phi2)                      # smoothing window
+        Q[ff, kpmin[ff]-ifreq[0]:kpmax[ff]-ifreq[0]+1] = W
 
-        # Calculate window
+    # Normalize to conserve input power
-        xg = np.arange(xl, xu)  # indices
+    Qpower = np.sum(Q, axis=0)
-        fg = freq[xg]           # [Hz] corresponding freq
+    Q = Q / Qpower[np.newaxis, :]
-        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
 
@@ -145,57 +155,40 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth):
 def smoothSpectralData(freq, M, sw: SmoothingWidth,
                        st: SmoothingType = SmoothingType.levels):
     """
-    Apply fractional octave smoothing to magnitude data in frequency domain.
+    Apply fractional octave smoothing to data in the 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
-
-    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
+        M: array of data, either power or dB
         the smoothing type `st`, the smoothing is applied.
+        sw: smoothing width
+        st: smoothing type = data type of input data
 
     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 len(MM) > 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'
+    assert len(freq) == len(MM), "f and M should have equal length"
 
-#    if st == SmoothingType.ps:
+    if st == SmoothingType.ps:
-#        assert np.min(M) >= 0, 'absolute magnitude M cannot be negative'
+        assert np.min(MM) >= 0, 'Power spectrum values cannot be negative'
-    if st == SmoothingType.levels and isinstance(M.dtype, complex):
+    if st == SmoothingType.levels and isinstance(MM.dtype, complex):
         raise RuntimeError('Decibel input should be real-valued')
 
-    # Initialize
+    # Convert to power
-    L = M.shape[0]                   # number of data points
-
     if st == SmoothingType.levels:
-        P = 10**(MM/10)  # magnitude [dB] --> power
+        P = 10**(MM/10)
+    elif st == SmoothingType.ps:
+        P = MM
     else:
-        P = MM  # data already given as power
+        raise RuntimeError(f"Incorrect SmoothingType: {st}")
-    # 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
@@ -221,6 +214,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     # Apply smoothing
     Psm = np.matmul(Q, P)
 
+    # Convert to original format
     if st == SmoothingType.levels:
         Psm = 10*np.log10(Psm)
 
@@ -243,7 +237,8 @@ if __name__ == "__main__":
     fmax = 24e3  # [Hz] max freq
     Ndata = 200  # number of data points
     freq = np.linspace(fmin, fmax, Ndata)  # frequency points
-    M = abs(0.4*np.random.normal(size=(Ndata,)))+0.01  #
+    # 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: dirac delta