From 5caddec58323bbb9f0cf197e5d286b536ce333d2 Mon Sep 17 00:00:00 2001 From: Casper Date: Thu, 23 Feb 2023 16:33:41 +0100 Subject: [PATCH] Alternative Smoothing: transformed to matrix method, minor bug fixes, merged test scripts --- src/lasp/tools/tools.py | 397 ++++++++++++++++++++++++++++++---------- 1 file changed, 303 insertions(+), 94 deletions(-) diff --git a/src/lasp/tools/tools.py b/src/lasp/tools/tools.py index 8329f46..20921a9 100644 --- a/src/lasp/tools/tools.py +++ b/src/lasp/tools/tools.py @@ -5,13 +5,10 @@ Author: C. Jansen, J.A. de Jong - ASCEE V.O.F. Smooth data in the frequency domain. -TODO: This function is rather slow as it -used Python for loops. The implementations should be speed up in the near -future. -TODO: check if the smoothing is correct: depending on whether the data points -are spaced lin of log in frequency, they should be given different weights. -TODO: accept data that is not equally spaced in frequency -TODO: output data that is log spaced in frequency +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 Cutoff frequencies of window taken from http://www.huennebeck-online.de/software/download/src/index.html 15-10-2021 @@ -30,14 +27,13 @@ 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', 'smoothSpectralDataAlt', 'intHann'] +__all__ = ['SmoothingType', 'smoothSpectralData', 'SmoothingWidth'] 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): @@ -95,14 +91,17 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth): 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 + # Loop over indices of raw frequency vector for x in range(x0, L): # Find indices of data points to calculate current (smoothed) magnitude @@ -129,9 +128,9 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth): # Find indices corresponding to frequencies xl = bisect.bisect_left(freq, fl) # index corresponding to fl - xu = bisect.bisect_left(freq, fu) + xu = bisect.bisect_right(freq, fu) - xl = xu-1 if xu-xl <= 0 else xl # Guarantee window length of at least 1 + xl = xu-1 if xu-xl <= 0 else xl # Guarantee window length >= 1 # Calculate window xg = np.arange(xl, xu) # indices @@ -140,6 +139,10 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth): gs /= np.sum(gs) # normalize: integral=1 Q[x, xl:xu] = gs # add to matrix + # Normalize to conserve input power + Qpower = np.sum(Q, axis=0) + Q = Q / Qpower[np.newaxis, :] + return Q @@ -154,7 +157,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, 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 + 07-05-2021 Update 16-01-2023: speed up algorithm - Smoothing is performed using matrix multiplication @@ -227,63 +230,10 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, return Psm +# %% Alternative algorithm +from numpy import arange, log2, log10, pi, ceil, floor, sin -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 +# Integrated Hann window def intHann(x1, x2): if (x2 <= -1/2) or (x1 >= 1/2): return 0 @@ -299,6 +249,242 @@ def intHann(x1, x2): 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): + """ + 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 + 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 + + Noct /= 1.5 # empirical correction factor: window @ -6 dB at Noct bounds + 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, len(M)): # loop over input freq + if kpmin[ff] < ifreqMin: + kpmin[ff] = ifreqMin + kpmax[ff] = 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] = floor(ifreq[ff]**2/ifreqMax) # decrease smooth. width + kpmax[ff] = ifreqMax + if np.isclose(kpmin[ff], kpmax[ff]): + kpmin[ff] -= 1 + NoctAct = 1/log2(kpmax[ff]/kpmin[ff]) + else: + NoctAct = Noct + + kp = arange(kpmin[ff], kpmax[ff]+1) # freqs of window + + # Integration bounds for Hann window + Phi1 = log2((kp - 0.5)/ifreq[ff]) * NoctAct + 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]) + + # 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 + + # Normalize to conserve input power + Qpower = np.sum(Q, axis=0) + Q = Q / Qpower[np.newaxis, :] + + return Q + + +def smoothSpectralDataAltMatrix(freq, M, sw: SmoothingWidth, + st: SmoothingType = SmoothingType.levels): + """ + 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). + 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 + the smoothing type `st`, the smoothing is applied. + + 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 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' + +# if st == SmoothingType.ps: +# assert np.min(M) >= 0, 'absolute magnitude M cannot be negative' + if st == SmoothingType.levels and isinstance(M.dtype, complex): + raise RuntimeError('Decibel input should be real-valued') + + # Initialize + L = M.shape[0] # number of data points + + if st == SmoothingType.levels: + 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 + + # 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 = smoothCalcMatrixAlt(freq, sw) + # Qdict[key] = Q + + Q = smoothCalcMatrixAlt(freq, sw) + + # Apply smoothing + Psm = np.matmul(Q, P) + + if st == SmoothingType.levels: + Psm = 10*np.log10(Psm) + + return Psm + # %% Test if __name__ == "__main__": """ Test function for evaluation and debugging @@ -307,53 +493,76 @@ if __name__ == "__main__": points. They should be treated and weighted differently. """ import matplotlib.pyplot as plt + import time + plt.close('all') + # Initialize - Noct = 2 # Noct = 6 for 1/6 oct. smoothing + Noct = 3 # 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 + # # 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: single tone + fmin = 0 # [Hz] min freq + fmax = 5e3 # [Hz] max freq + Ndata = 2501 # number of data points freq = np.linspace(fmin, fmax, Ndata) # frequency points - M = abs(0.4*np.random.normal(size=(Ndata,)))+0.01 # - M = 20*np.log10(M) + M = 1e-4*np.random.normal(size=(Ndata,)) + M[500] = 1 + MdB = 20*np.log10(abs(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) + if 'Qdict' in globals(): + del Qdict + + t0 = time.time() + Msm = smoothSpectralData(freq, MdB, sw, st) # current algorithm + t1 = time.time() + MsmAlt = smoothSpectralDataAlt(freq, MdB, sw, st) # alternative algorithm + 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)") # Plot - lin frequency plt.figure() - plt.plot(freq, M, '.b') + 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([100, fmax]) + plt.xlim((0, fmax)) + plt.ylim((-90, 1)) + plt.grid('both') plt.title('lin frequency') - plt.legend(['Raw', 'Smooth']) + plt.legend(['Raw', 'Smooth', 'SmoothAlt', 'SmoothAltMatrix']) # Plot - log frequency plt.figure() - plt.semilogx(freq, M, '.b') + 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.xlim((100, fmax)) + plt.ylim((-90, 1)) + plt.grid('both') plt.title('log frequency') - plt.legend(['Raw', 'Smooth 1']) + plt.legend(['Raw', 'Smooth', 'SmoothAlt', 'SmoothAltMatrix'])