diff --git a/src/lasp/tools/tools.py b/src/lasp/tools/tools.py index 20921a9..c0ffc8b 100644 --- a/src/lasp/tools/tools.py +++ b/src/lasp/tools/tools.py @@ -1,7 +1,7 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ -Author: C. Jansen, J.A. de Jong - ASCEE V.O.F. +Author: T. Hekman, J.A. de Jong, C. Jansen - ASCEE V.O.F. Smooth data in the frequency domain. @@ -174,6 +174,8 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, """ # TODO: Make this function multi-dimensional array aware. + # TODO: This does not work due to complex numbers. Should be split up in + # magnitude and phase. # Safety MM = copy.deepcopy(M) @@ -196,8 +198,6 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, 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 @@ -235,85 +235,20 @@ from numpy import arange, log2, log10, pi, ceil, floor, sin # 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) - - -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 + 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 """ - Noct = 1/sw.value[0] - # M = MdB - M = 10**(MdB/20) + 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) - 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): @@ -336,7 +271,6 @@ def smoothCalcMatrixAlt(freq, sw: SmoothingWidth): 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 "\ @@ -350,7 +284,8 @@ def smoothCalcMatrixAlt(freq, sw: SmoothingWidth): 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 + # Noct /= 1.5 # empirical correction factor: window @ -6 dB at Noct bounds + Noct /= 2 # empirical correction factor: window @ -3 dB at Noct bounds ifreq = freq/(freq[1]-freq[0]) # frequency, normalized to step=1 ifreq = np.array(ifreq.astype(int)) @@ -362,7 +297,9 @@ def smoothCalcMatrixAlt(freq, sw: SmoothingWidth): 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 + # Find window bounds and actual smoothing width if kpmin[ff] < ifreqMin: kpmin[ff] = ifreqMin kpmax[ff] = ceil(ifreq[ff]**2/ifreqMin) # decrease smooth. width @@ -385,9 +322,7 @@ def smoothCalcMatrixAlt(freq, sw: SmoothingWidth): 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]) + W = intHann(Phi1, Phi2) # 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 @@ -497,7 +432,7 @@ if __name__ == "__main__": plt.close('all') # Initialize - Noct = 3 # Noct = 6 for 1/6 oct. smoothing + Noct = 1 # Noct = 6 for 1/6 oct. smoothing # # Create dummy data set 1: noise # fmin = 1e3 # [Hz] min freq @@ -529,40 +464,35 @@ if __name__ == "__main__": t0 = time.time() Msm = smoothSpectralData(freq, MdB, sw, st) # current algorithm t1 = time.time() - MsmAlt = smoothSpectralDataAlt(freq, MdB, sw, st) # alternative algorithm + MsmAlt = smoothSpectralDataAltMatrix(freq, MdB, sw, st) # alternative algorithm, matrix method 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)") + print(f"Smoothing time: {t1-t0} s (Current)") + print(f"Smoothing time: {t2-t1} s (Alternative)") # Plot - lin frequency plt.figure() 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((0, fmax)) plt.ylim((-90, 1)) plt.grid('both') plt.title('lin frequency') - plt.legend(['Raw', 'Smooth', 'SmoothAlt', 'SmoothAltMatrix']) + plt.legend(['Raw', 'Smooth', 'SmoothAlt']) # Plot - log frequency plt.figure() 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.ylim((-90, 1)) plt.grid('both') plt.title('log frequency') - plt.legend(['Raw', 'Smooth', 'SmoothAlt', 'SmoothAltMatrix']) + plt.legend(['Raw', 'Smooth', 'SmoothAlt'])