tools.py

Go to the documentation of this file.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
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 uses [for loops] in Python. Speed up.
NOTE: function requires lin frequency spaced input data
TODO: accept input data that is not lin spaced in frequency
TODO: it makes more sense to output data that is log spaced in frequency
TODO: Make SmoothSpectralData() multi-dimensional array aware.
TODO: Smoothing does not work due to complex numbers. Is it reasonable to
smooth complex data? If so, the data could be split in magnitude and phase.
"""
 
__all__ = ['SmoothingType', 'smoothSpectralData', 'SmoothingWidth']
 
from enum import Enum, unique
import copy
import numpy as np
from numpy import log2, pi, sin
 
 
@unique
class SmoothingWidth(Enum):
    none = (0, 'No smoothing')
    one = (1, '1/1st octave smoothing')
    two = (2, '1/2th octave smoothing')
    three = (3, '1/3rd octave smoothing')
    six = (6, '1/6th octave smoothing')
    twelve = (12, '1/12th octave smoothing')
    twfo = (24, '1/24th octave smoothing')
    ftei = (48, '1/48th octave smoothing')
    hundred = (100, '1/100th octave smoothing')  # useful for removing 'grass'
 
    @staticmethod
    def fillComboBox(cb):
        """
        Fill Windows to a combobox
 
        Args:
            cb: QComboBox to fill
        """
        cb.clear()
        for w in list(SmoothingWidth):
            cb.addItem(w.value[1], w)
        cb.setCurrentIndex(0)
 
    @staticmethod
    def getCurrent(cb):
        return list(SmoothingWidth)[cb.currentIndex()]
 
 
class SmoothingType:
    levels = 'l', 'Levels'  # [dB]
    # tf = 'tf', 'Transfer function',
    ps = 'ps', '(Auto) powers'
 
 
# 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):
    """
    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
    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 /= 2  # corr. factor: window @ -3dB at Noct bounds (Noct/1.5 for -6dB)
    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, L):  # loop over input freq
        # Find window bounds and actual smoothing width
        # 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
 
        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
 
    # Normalize to conserve input power
    Qpower = np.sum(Q, axis=0)
    Q = Q / Qpower[np.newaxis, :]
 
    return Q
 
 
def smoothSpectralData(freq, M, sw: SmoothingWidth,
                       st: SmoothingType = SmoothingType.levels):
    """
    Apply fractional octave smoothing to data in the frequency domain.
 
    Args:
        freq: array of frequencies of data points [Hz] - equally spaced
        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
    """
    # Safety
    MM = copy.deepcopy(M)
    Noct = sw.value[0]
    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(MM), "f and M should have equal length"
 
    if st == SmoothingType.ps:
        assert np.min(MM) >= 0, 'Power spectrum values cannot be negative'
    if st == SmoothingType.levels and isinstance(MM.dtype, complex):
        raise RuntimeError('Decibel input should be real-valued')
 
    # Convert to power
    if st == SmoothingType.levels:
        P = 10**(MM/10)
    elif st == SmoothingType.ps:
        P = MM
    else:
        raise RuntimeError(f"Incorrect SmoothingType: {st}")
 
    # 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 = smoothCalcMatrix(freq, sw)
        Qdict[key] = Q
 
    # Apply smoothing
    Psm = np.matmul(Q, P)
 
    # Convert to original format
    if st == SmoothingType.levels:
        Psm = 10*np.log10(Psm)
 
    return Psm
 
 
# %% Test
if __name__ == "__main__":
    """ Test function for evaluation and debugging
 
    Note: make a distinction between lin and log spaced (in frequency) data
    points. They should be treated and weighted differently.
    """
    import matplotlib.pyplot as plt
    # Initialize
    Noct = 2  # Noct = 6 for 1/6 oct. smoothing
 
    # Create dummy data set 1: noise
    fmin = 1e3   # [Hz] min freq
    fmax = 24e3  # [Hz] max freq
    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: 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)
    fsm = freq
 
 
    # Plot - lin frequency
    plt.figure()
    plt.plot(freq, M, '.b')
    plt.plot(fsm, Msm,  'r')
    plt.xlabel('f (Hz)')
    plt.ylabel('magnitude')
    plt.xlim([100, fmax])
    plt.title('lin frequency')
    plt.legend(['Raw', 'Smooth'])
 
    # Plot - log frequency
    plt.figure()
    plt.semilogx(freq, M, '.b')
    plt.semilogx(fsm, Msm, 'r')
    plt.xlabel('f (Hz)')
    plt.ylabel('magnitude')
    plt.xlim([100, fmax])
    plt.title('log frequency')
    plt.legend(['Raw', 'Smooth 1'])