diff --git a/src/lasp/tools/tools.py b/src/lasp/tools/tools.py index 90a8be9..f8a938b 100644 --- a/src/lasp/tools/tools.py +++ b/src/lasp/tools/tools.py @@ -78,6 +78,70 @@ class SmoothingType: # TO DO: add possibility to insert data that is not lin spaced in frequency +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) + """ + # 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" + 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 + + 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 + # + # Indices beyond [0, L] point to non-existing data. Beyond 0 does not + # 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 + # there is no data beyond the Nyquist frequency, also truncate the + # 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 + xl = bisect.bisect_left(freq, fl) # index corresponding to fl + xu = bisect.bisect_left(freq, fu) + + xl = xu-1 if xu-xl <= 0 else xl # Guarantee window length of at least 1 + + # Calculate window + xg = np.arange(xl, xu) # indices + fg = freq[xg] # [Hz] corresponding freq + 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 + + def smoothSpectralData(freq, M, sw: SmoothingWidth, st: SmoothingType = SmoothingType.levels): """ @@ -91,6 +155,10 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, 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 @@ -103,9 +171,6 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, """ # TODO: Make this function multi-dimensional array aware. - # Settings - tr = 2 # truncate window after 2x std; shorter is faster and less accurate - # Safety MM = copy.deepcopy(M) Noct = sw.value[0] @@ -134,169 +199,33 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, # 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.. + if freq[0] == 0: + 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 - # - # Indices beyond [0, L] point to non-existing data. Beyond 0 does not - # 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 + # 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 the upper (frequency) side of the window is truncated because - # there is no data beyond the Nyquist frequency, also truncate the - # 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 + if 'Qdict' not in globals(): # Guarantee Qdict exists + global Qdict + Qdict = {} - # Find indices corresponding to frequencies - xl = bisect.bisect_left(freq, fl) # index corresponding to fl - xu = bisect.bisect_left(freq, fu) + if key in Qdict: + Q = Qdict[key] + else: + Q = smoothCalcMatrix(freq, sw) + Qdict[key] = Q - # Guarantee window length of at least 1 - if xu - xl <= 0: - xl = xu - 1 - - # Calculate window - g = np.zeros(xu-xl) - for n, xi in enumerate(range(xl, xu)): - fi = freq[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 - Psm[x] = np.dot(g, P[xl:xu]) + # Apply smoothing + Psm = np.matmul(Q, P) 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__": @@ -310,7 +239,7 @@ if __name__ == "__main__": Noct = 2 # Noct = 6 for 1/6 oct. smoothing # Create dummy data set 1: noise - fmin = 3e3 # [Hz] min freq + fmin = 1e3 # [Hz] min freq fmax = 24e3 # [Hz] max freq Ndata = 200 # number of data points freq = np.linspace(fmin, fmax, Ndata) # frequency points @@ -330,10 +259,13 @@ if __name__ == "__main__": 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') @@ -342,12 +274,14 @@ if __name__ == "__main__": 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.semilogx(fsm, Msm, 'r') plt.xlabel('f (Hz)') plt.ylabel('magnitude') plt.xlim([100, fmax]) plt.title('log frequency') + plt.legend(['Raw', 'Smooth 1'])