From 8912cb145c8462a8ee46bfb9ebfb1589decfed99 Mon Sep 17 00:00:00 2001 From: Casper Date: Fri, 15 Oct 2021 15:33:25 +0200 Subject: [PATCH] Updated octave smoothing function; completely different approach --- lasp/tools/tools.py | 271 +++++++++++++++++++++++++++++++++----------- 1 file changed, 203 insertions(+), 68 deletions(-) diff --git a/lasp/tools/tools.py b/lasp/tools/tools.py index c158b3c..47173bb 100644 --- a/lasp/tools/tools.py +++ b/lasp/tools/tools.py @@ -3,19 +3,46 @@ """ 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 +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 + +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 """ -from enum import Enum, unique __all__ = ['SmoothingType', 'smoothSpectralData', 'SmoothingWidth'] +from enum import Enum, unique +import bisect +import numpy as np + @unique class SmoothingWidth(Enum): none = (0, 'No smoothing') - # three = (3, '1/3th octave smoothing') + one = (1, '1/1stoctave 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') @@ -38,6 +65,7 @@ class SmoothingWidth(Enum): def getCurrent(cb): return list(SmoothingWidth)[cb.currentIndex()] + class SmoothingType: levels = 'l', 'Levels' # tf = 'tf', 'Transfer function', @@ -48,9 +76,6 @@ class SmoothingType: # TO DO: add possibility to insert data that is not lin spaced in frequency -import numpy as np -from scipy.signal.windows import gaussian - def smoothSpectralData(freq, M, sw: SmoothingWidth, st: SmoothingType = SmoothingType.levels): """ @@ -77,13 +102,14 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, # TODO: Make this function multi-dimensional array aware. # Settings - tr = 2 # truncate window after 2x std + tr = 2 # truncate window after 2x std; shorter is faster and less accurate # Safety Noct = sw.value[0] 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 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' @@ -92,109 +118,218 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, # Initialize L = M.shape[0] # number of data points - - P = M + if st == SmoothingType.levels: - P = 10**(P/10) + P = 10**(M/10) # magnitude [dB] --> power + else: + P = M # 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 - x0 = 1 if freq[0]==0 else 0 # Skip first data point if zero frequency - Psm[0] = P[0] # Reuse old value in case first data.. - # ..point is skipped. Not plotted any way. - df = freq[1] - freq[0] # Frequency step + Psm = np.zeros_like(P) # Smoothed power - to be calculated + x0 = 1 if freq[0] == 0 else 0 # Skip first data point if zero frequency + Psm[0] = P[0] # Reuse old value in case first data.. + # ..point is skipped. Not plotted any way. # Loop through data points for x in range(x0, L): # Find indices of data points to calculate current (smoothed) magnitude - fc = freq[x] # center freq. of smoothing window - Df = tr * fc / Noct # freq. range of smoothing window - - # desired lower index of frequency array to be used during smoothing - xl = int(np.ceil(x - 0.5*Df/df)) - - # upper index + 1 (because half-open interval) - xu = int(np.floor(x + 0.5*Df/df)) + 1 - - # 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 + # + # Implicitly, the window is truncated in this implementation, because + # of the bisect search + 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 + xl = bisect.bisect_left(f, fl) # index corresponding to fl + xu = bisect.bisect_left(f, fu) - # 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] + # Calculate window + g = np.zeros(xu-xl) + for n, xi in enumerate(range(xl, xu)): + fi = f[xi] # current frequency + g[n] = np.sqrt( 1/ ((1+((fi/fc - fc/fi)*(1.507*Noct))**6)) ) + g /= np.sum(g) # normalize: integral=1 # Apply smoothing - wind_int = np.sum(wind) # integral - Psm[x] = np.dot(wind, P[xl:xu]) / wind_int # apply window + Psm[x] = np.dot(g, P[xl:xu]) if st == SmoothingType.levels: Psm = 10*np.log10(Psm) return Psm +## OLD VERSION +#from scipy.signal.windows import gaussian +#def smoothSpectralData(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). +# Casper Jansen, 07-05-2021 +# +# 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. +# +# # Settings +# tr = 2 # truncate window after 2x std; shorter is faster and less accurate +# +# # Safety +# Noct = sw.value[0] +# 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 +# +# P = M +# if st == SmoothingType.levels: +# P = 10**(P/10) +# # 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 +# x0 = 1 if freq[0]==0 else 0 # Skip first data point if zero frequency +# Psm[0] = P[0] # Reuse old value in case first data.. +# # ..point is skipped. Not plotted any way. +# df = freq[1] - freq[0] # Frequency step +# +# # Loop through data points +# for x in range(x0, L): +# # Find indices of data points to calculate current (smoothed) magnitude +# fc = freq[x] # center freq. of smoothing window +# Df = tr * fc / Noct # freq. range of smoothing window +# +# # desired lower index of frequency array to be used during smoothing +# xl = int(np.ceil(x - 0.5*Df/df)) +# +# # upper index + 1 (because half-open interval) +# xu = int(np.floor(x + 0.5*Df/df)) + 1 +# +# # Create window, suitable for frequency lin-spaced data points +# 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 +# +# 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 = 6 # 1/6 oct. smoothing + Noct = 6 # Noct = 6 for 1/6 oct. smoothing - # Create dummy data - fmin = 3e3 # [Hz] min freq + # Create dummy data set 1: noise + fmin = 3e3 # [Hz] min freq fmax = 24e3 # [Hz] max freq - Ndata = 200 # number of data points - + 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 +# # Create dummy data set 2: dirac delta +# 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 = 0 * abs(1+0.4*np.random.normal(size=(Ndata,))) + 0.01 # +# M[int(Ndata/5)] = 1 +# M = 20*np.log10(M) # Apply function - Msm = oct_smooth(f, M, Noct, dB) + class sw: + value = [Noct] + st = SmoothingType.levels # so data is given in dB + + Msm = smoothSpectralData(f, M, sw, st) fsm = f - # Plot - + # Plot - lin frequency 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]) + plt.title('lin frequency') + + # Plot - log frequency + plt.figure() + plt.semilogx(f, M, '.b') + plt.semilogx(fsm, Msm, 'r') + plt.xlabel('f (Hz)') + plt.ylabel('magnitude') + plt.xlim([100, fmax]) + plt.title('log frequency')