ASCEE/lasp
ASCEE/lasp
5
1
Fork 0
Code Issues 5 Pull Requests Actions Releases 14 Activity

Fixed a bug in the new smoothing algorithm

Browse Source
This commit is contained in:
Casper Jansen 2021-10-18 08:23:42 +02:00
parent 417c1bd3d8
commit 82fc38a990
1 changed files with 21 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
lasp/tools/tools.py
View File

@ -108,7 +108,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
# Safety
MM = copy.deepcopy(M)
Noct = sw.value[0]
assert M, "Smoothing function: input array is empty"
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')
@ -141,17 +141,26 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
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. 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.
#
# Implicitly, the window is truncated in this implementation, because
# of the bisect search
# 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
xl = bisect.bisect_left(f, fl) # index corresponding to fl
xu = bisect.bisect_left(f, fu)
@ -292,14 +301,14 @@ if __name__ == "__main__":
"""
import matplotlib.pyplot as plt
# Initialize
Noct = 6 # Noct = 6 for 1/6 oct. smoothing
Noct = 2 # Noct = 6 for 1/6 oct. smoothing
# Create dummy data set 1: noise
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 = abs(1+0.4*np.random.normal(size=(Ndata,)))+0.01 #
M = abs(0.4*np.random.normal(size=(Ndata,)))+0.01 #
M = 20*np.log10(M)
# # Create dummy data set 2: dirac delta
@ -308,7 +317,7 @@ if __name__ == "__main__":
# 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[int(100)] = 1
# M = 20*np.log10(M)
# Apply function