From b3fb7ddb6d1c9c0a30205140311b98a256be7fc1 Mon Sep 17 00:00:00 2001 From: Casper Date: Thu, 23 Feb 2023 16:37:08 +0100 Subject: [PATCH 1/3] Smoothing: minor bug fix --- src/lasp/tools/tools.py | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/src/lasp/tools/tools.py b/src/lasp/tools/tools.py index f8a938b..c17b0af 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 @@ -91,17 +88,20 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth): Warning: this method does not work on levels (dB) """ # Settings - tr = 2 # truncate window after 2x std; shorter is faster and less accurate + tr = 3 # 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 @@ -128,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 @@ -139,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 @@ -147,13 +151,13 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, """ 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. + variable length. The window is truncated after 3x 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 + 07-05-2021 Update 16-01-2023: speed up algorithm - Smoothing is performed using matrix multiplication @@ -243,7 +247,8 @@ if __name__ == "__main__": fmax = 24e3 # [Hz] max freq Ndata = 200 # number of data points freq = np.linspace(fmin, fmax, Ndata) # frequency points - M = abs(0.4*np.random.normal(size=(Ndata,)))+0.01 # + # 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 From fa8f5e64ad83ef5bb3c2effea3c50a8908930452 Mon Sep 17 00:00:00 2001 From: Casper Date: Thu, 23 Feb 2023 18:10:06 +0100 Subject: [PATCH 2/3] Update to new smoothing algorithm. Should be made faster. --- src/lasp/tools/tools.py | 160 +++++++++++++++++++--------------------- 1 file changed, 75 insertions(+), 85 deletions(-) diff --git a/src/lasp/tools/tools.py b/src/lasp/tools/tools.py index c17b0af..1454dc9 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, C. Jansen, J.A. de Jong - ASCEE V.O.F. Smooth data in the frequency domain. @@ -9,30 +9,17 @@ 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 -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 +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 bisect import copy import numpy as np +from numpy import log2, pi, sin @unique @@ -73,6 +60,22 @@ class SmoothingType: # 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): @@ -83,12 +86,18 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth): Returns: freq: array frequencies of data points [Hz] - Q: matrix to smooth power: {fsm} = [Q] * {fraw} + 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 = 3 # 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 "\ @@ -100,44 +109,41 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth): 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 + 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 + 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 - # 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 + kpmin = np.floor(ifreq/sfact).astype(int) # min freq of window + kpmax = np.ceil(ifreq*sfact).astype(int) # max freq of window - # Find indices corresponding to frequencies - xl = bisect.bisect_left(freq, fl) # index corresponding to fl - xu = bisect.bisect_right(freq, fu) + 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 - xl = xu-1 if xu-xl <= 0 else xl # Guarantee window length >= 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 + 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) @@ -149,57 +155,40 @@ def smoothCalcMatrix(freq, sw: SmoothingWidth): 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 3x 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 + Apply fractional octave smoothing to data in the frequency domain. Args: freq: array of frequencies of data points [Hz] - equally spaced - M: array of either power, transfer functin or dB points. Depending on + 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 - """ - # 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 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(M), 'f and M should have equal length' + assert len(freq) == len(MM), "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): + 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') - # Initialize - L = M.shape[0] # number of data points - + # Convert to power if st == SmoothingType.levels: - P = 10**(MM/10) # magnitude [dB] --> power + P = 10**(MM/10) + elif st == SmoothingType.ps: + P = MM 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 + 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 @@ -225,6 +214,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth, # Apply smoothing Psm = np.matmul(Q, P) + # Convert to original format if st == SmoothingType.levels: Psm = 10*np.log10(Psm) From 8547d0915ab51c94187c3aee58d2e00e11f4e820 Mon Sep 17 00:00:00 2001 From: Thijs Hekman Date: Thu, 9 Mar 2023 10:16:44 +0100 Subject: [PATCH 3/3] Added low-pass and high-pass compensator filters --- src/lasp/filter/biquad.py | 77 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 75 insertions(+), 2 deletions(-) diff --git a/src/lasp/filter/biquad.py b/src/lasp/filter/biquad.py index dfff947..1ce16f5 100644 --- a/src/lasp/filter/biquad.py +++ b/src/lasp/filter/biquad.py @@ -23,11 +23,11 @@ y[n] = 1/ba[3] * ( ba[0] * x[n] + ba[1] * x[n-1] + ba[2] * x[n-2] + """ __all__ = ['peaking', 'biquadTF', 'notch', 'lowpass', 'highpass', - 'highshelf', 'lowshelf'] + 'highshelf', 'lowshelf', 'LPcompensator', 'HPcompensator'] from numpy import array, cos, pi, sin, sqrt from scipy.interpolate import interp1d -from scipy.signal import sosfreqz +from scipy.signal import sosfreqz, bilinear_zpk, zpk2sos def peaking(fs, f0, Q, gain): @@ -157,6 +157,79 @@ def lowshelf(fs, f0, Q, gain): a2 = (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha return array([b0/a0, b1/a0, b2/a0, a0/a0, a1/a0, a2/a0]) +def LPcompensator(fs, f0o, Qo, f0n, Qn): + """ + Shelving type filter that, when multiplied with a second-order low-pass + filter, alters the response of that filter to a different second-order + low-pass filter. + + Args: + fs: Sampling frequency [Hz] + f0o: Cut-off frequency of the original filter [Hz] + Qo: Quality factor of the original filter (~inverse of bandwidth) + f0n: Desired cut-off frequency [Hz] + Qn: Desired quality factor(~inverse of bandwidth) + """ + + omg0o = 2*pi*f0o + omg0n = 2*pi*f0n + + zRe = omg0o/(2*Qo) + zIm = omg0o*sqrt(1-1/(4*Qo**2)) + z1 = -zRe + zIm*1j + z2 = -zRe - zIm*1j + + pRe = omg0n/(2*Qn) + pIm = omg0n*sqrt(1-1/(4*Qn**2)) + p1 = -pRe + pIm*1j + p2 = -pRe - pIm*1j + + z= [z1, z2] + p = [p1, p2] + k = (pRe**2 + pIm**2)/(zRe**2 + zIm**2) + + zd, pd, kd = bilinear_zpk(z, p, k, fs) + + sos = zpk2sos(zd,pd,kd) + + return sos[0] + +def HPcompensator(fs, f0o, Qo, f0n, Qn): + """ + Shelving type filter that, when multiplied with a second-order high-pass + filter, alters the response of that filter to a different second-order + high-pass filter. + + Args: + fs: Sampling frequency [Hz] + f0o: Cut-on frequency of the original filter [Hz] + Qo: Quality factor of the original filter (~inverse of bandwidth) + f0n: Desired cut-on frequency [Hz] + Qn: Desired quality factor(~inverse of bandwidth) + """ + + omg0o = 2*pi*f0o + omg0n = 2*pi*f0n + + zRe = omg0o/(2*Qo) + zIm = omg0o*sqrt(1-1/(4*Qo**2)) + z1 = -zRe + zIm*1j + z2 = -zRe - zIm*1j + + pRe = omg0n/(2*Qn) + pIm = omg0n*sqrt(1-1/(4*Qn**2)) + p1 = -pRe + pIm*1j + p2 = -pRe - pIm*1j + + z= [z1, z2] + p = [p1, p2] + k = 1 + + zd, pd, kd = bilinear_zpk(z, p, k, fs) + + sos = zpk2sos(zd,pd,kd) + + return sos[0] def biquadTF(fs, freq, sos): """