Small interface changes for smoothing function
This commit is contained in:
parent
9a16219059
commit
f4e3688222
@ -9,3 +9,4 @@ from .lasp_slm import *
|
||||
from .lasp_record import *
|
||||
from .lasp_siggen import *
|
||||
from .lasp_weighcal import *
|
||||
from .tools import *
|
||||
|
||||
@ -18,7 +18,7 @@ Common definitions used throughout the code.
|
||||
|
||||
__all__ = [
|
||||
'P_REF', 'FreqWeighting', 'TimeWeighting', 'getTime', 'getFreq', 'Qty',
|
||||
'SIQtys',
|
||||
'SIQtys', 'Window',
|
||||
'lasp_shelve', 'this_lasp_shelve', 'W_REF', 'U_REF', 'I_REF', 'dBFS_REF',
|
||||
'AvType'
|
||||
]
|
||||
@ -263,16 +263,14 @@ class this_lasp_shelve(Shelve):
|
||||
return os.path.join(lasp_appdir, f'{node}_config.shelve')
|
||||
|
||||
|
||||
class Window:
|
||||
@unique
|
||||
class Window(Enum):
|
||||
hann = (wWindow.hann, 'Hann')
|
||||
hamming = (wWindow.hamming, 'Hamming')
|
||||
rectangular = (wWindow.rectangular, 'Rectangular')
|
||||
bartlett = (wWindow.bartlett, 'Bartlett')
|
||||
blackman = (wWindow.blackman, 'Blackman')
|
||||
|
||||
types = (hann, hamming, rectangular, bartlett, blackman)
|
||||
default = 0
|
||||
|
||||
@staticmethod
|
||||
def fillComboBox(cb):
|
||||
"""
|
||||
@ -282,12 +280,13 @@ class Window:
|
||||
cb: QComboBox to fill
|
||||
"""
|
||||
cb.clear()
|
||||
for tw in Window.types:
|
||||
cb.addItem(tw[1], tw)
|
||||
cb.setCurrentIndex(Window.default)
|
||||
for w in list(Window):
|
||||
cb.addItem(w.value[1], w)
|
||||
cb.setCurrentIndex(0)
|
||||
|
||||
@staticmethod
|
||||
def getCurrent(cb):
|
||||
return Window.types[cb.currentIndex()]
|
||||
return list(Window)[cb.currentIndex()]
|
||||
|
||||
|
||||
class TimeWeighting:
|
||||
|
||||
@ -1,6 +1,3 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
from .config import init_backend
|
||||
from lasp.lasp_common import FreqWeighting, TimeWeighting
|
||||
__all__ = ['init_backend',
|
||||
'FreqWeighting', 'TimeWeighting']
|
||||
from .tools import *
|
||||
|
||||
@ -1,52 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""!
|
||||
Author: J.A. de Jong - ASCEE
|
||||
|
||||
Description:
|
||||
"""
|
||||
__all__ = ['init_backend', 'getReportQuality']
|
||||
|
||||
_report_quality = False
|
||||
_init = False
|
||||
|
||||
|
||||
def init_matplotlib(report_quality=False):
|
||||
global _init
|
||||
if not _init:
|
||||
_init = True
|
||||
print('Initializing matplotlib...')
|
||||
preamble = [
|
||||
r'\usepackage{libertine-type1}'
|
||||
r'\usepackage[libertine]{newtxmath}'
|
||||
# r'\usepackage{fontspec}',
|
||||
# r'\setmainfont{Libertine}',
|
||||
]
|
||||
params = {
|
||||
'font.family': 'serif',
|
||||
'text.usetex': True,
|
||||
'text.latex.unicode': True,
|
||||
'pgf.rcfonts': False,
|
||||
'pgf.texsystem': 'pdflatex',
|
||||
'pgf.preamble': preamble,
|
||||
}
|
||||
import matplotlib
|
||||
matplotlib.rcParams.update(params)
|
||||
global _report_quality
|
||||
_report_quality = report_quality
|
||||
|
||||
|
||||
def init_backend(report_quality=False):
|
||||
global _init
|
||||
if not _init:
|
||||
import matplotlib
|
||||
matplotlib.use('Qt5Agg', warn=False, force=True)
|
||||
init_matplotlib(report_quality)
|
||||
_init = True
|
||||
import matplotlib.pyplot as plt
|
||||
plt.ion()
|
||||
|
||||
|
||||
def getReportQuality():
|
||||
global _report_quality
|
||||
return _report_quality
|
||||
@ -1,210 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""!
|
||||
Author: J.A. de Jong - ASCEE
|
||||
|
||||
Description:
|
||||
"""
|
||||
__all__ = ['Figure', 'Bode', 'PS', 'PSD']
|
||||
|
||||
from .config import getReportQuality
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from cycler import cycler
|
||||
from lasp.lasp_common import (PLOT_COLORS_LIST, PLOT_NOCOLORS_LIST,
|
||||
DEFAULT_FIGSIZE_H, DEFAULT_FIGSIZE_W)
|
||||
|
||||
|
||||
class Figure:
|
||||
def __init__(self, **kwargs):
|
||||
ncols = kwargs.pop('ncols', 1)
|
||||
nrows = kwargs.pop('nrows', 1)
|
||||
color = kwargs.pop('color', (PLOT_NOCOLORS_LIST if getReportQuality()
|
||||
else PLOT_COLORS_LIST))
|
||||
if isinstance(color, bool):
|
||||
if color:
|
||||
color = PLOT_COLORS_LIST
|
||||
else:
|
||||
color = PLOT_NOCOLORS_LIST
|
||||
colors = cycler('color', color)
|
||||
|
||||
figsize = kwargs.pop('figsize', (DEFAULT_FIGSIZE_W, DEFAULT_FIGSIZE_H))
|
||||
self._f = plt.figure(figsize=figsize)
|
||||
|
||||
marker = kwargs.pop('marker', False)
|
||||
if marker:
|
||||
markers = cycler(marker=['o', 's', 'D', 'X', 'v', '^', '<', '>'])
|
||||
else:
|
||||
markers = cycler(marker=[None]*8)
|
||||
|
||||
linewidths = cycler(linewidth=[1, 2, 1, 2, 2, 3, 2, 1])
|
||||
|
||||
linestyles = cycler(
|
||||
linestyle=['-', '-', '--', ':', '-', '--', ':', '-.', ])
|
||||
|
||||
self._ax = []
|
||||
self._legend = {}
|
||||
for row in range(nrows):
|
||||
self._legend[row] = {}
|
||||
for col in range(ncols):
|
||||
self._legend[row][col] = []
|
||||
ax = self._f.add_subplot(100*nrows
|
||||
+ 10*ncols
|
||||
+ (row*ncols + col)+1)
|
||||
ax.set_prop_cycle(
|
||||
colors+linestyles+markers+linewidths)
|
||||
self._ax.append(ax)
|
||||
self._ncols = ncols
|
||||
self._cur_ax = self._ax[0]
|
||||
self._cur_col = 0
|
||||
self._cur_row = 0
|
||||
|
||||
self._zorder = -1
|
||||
|
||||
def setAx(self, row, col):
|
||||
self._cur_ax = self._ax[row*self._ncols+col]
|
||||
|
||||
@property
|
||||
def fig(self):
|
||||
return self._f
|
||||
|
||||
def markup(self):
|
||||
for ax in self._ax:
|
||||
ax.grid(True, 'both')
|
||||
self._zorder -= 1
|
||||
self.fig.show()
|
||||
|
||||
def vline(self, x):
|
||||
self._ax[0].axvline(x)
|
||||
|
||||
def plot(self, *args, **kwargs):
|
||||
line = self._cur_ax.plot(*args, **kwargs, zorder=self._zorder)
|
||||
self.markup()
|
||||
return line
|
||||
|
||||
def loglog(self, *args, **kwargs):
|
||||
line = self._cur_ax.loglog(*args, **kwargs, zorder=self._zorder)
|
||||
self.markup()
|
||||
return line
|
||||
|
||||
def semilogx(self, *args, **kwargs):
|
||||
line = self._cur_ax.semilogx(*args, **kwargs, zorder=self._zorder)
|
||||
self.markup()
|
||||
return line
|
||||
|
||||
def xlabel(self, *args, **kwargs):
|
||||
all_ax = kwargs.pop('all_ax', False)
|
||||
if all_ax:
|
||||
for ax in self._ax:
|
||||
ax.set_xlabel(*args, **kwargs)
|
||||
else:
|
||||
self._cur_ax.set_xlabel(*args, **kwargs)
|
||||
|
||||
def ylabel(self, *args, **kwargs):
|
||||
all_ax = kwargs.pop('all_ax', False)
|
||||
if all_ax:
|
||||
for ax in self._ax:
|
||||
ax.set_ylabel(*args, **kwargs)
|
||||
else:
|
||||
self._cur_ax.set_ylabel(*args, **kwargs)
|
||||
|
||||
def legend(self, leg, *args, **kwargs):
|
||||
# all_ax = kwargs.pop('all_ax', False)
|
||||
if isinstance(leg, list) or isinstance(leg, tuple):
|
||||
self._legend[self._cur_col][self._cur_col] = list(leg)
|
||||
else:
|
||||
self._legend[self._cur_col][self._cur_col].append(leg)
|
||||
self._cur_ax.legend(self._legend[self._cur_col][self._cur_col])
|
||||
|
||||
def savefig(self, *args, **kwargs):
|
||||
self.fig.savefig(*args, **kwargs)
|
||||
|
||||
def xlim(self, *args, **kwargs):
|
||||
all_ax = kwargs.pop('all_ax', False)
|
||||
if all_ax:
|
||||
for ax in self._ax:
|
||||
ax.set_xlim(*args, **kwargs)
|
||||
else:
|
||||
self._cur_ax.set_xlim(*args, **kwargs)
|
||||
|
||||
def ylim(self, *args, **kwargs):
|
||||
all_ax = kwargs.pop('all_ax', False)
|
||||
if all_ax:
|
||||
for ax in self._ax:
|
||||
ax.set_ylim(*args, **kwargs)
|
||||
else:
|
||||
self._cur_ax.set_ylim(*args, **kwargs)
|
||||
|
||||
def title(self, *args, **kwargs):
|
||||
self._cur_ax.set_title(*args, **kwargs)
|
||||
|
||||
def xticks(self, ticks):
|
||||
for ax in self._ax:
|
||||
ax.set_xticks(ticks)
|
||||
|
||||
def close(self):
|
||||
plt.close(self._f)
|
||||
|
||||
def xscale(self, scale):
|
||||
for ax in self._ax:
|
||||
ax.set_xscale(scale)
|
||||
|
||||
|
||||
class Bode(Figure):
|
||||
def __init__(self, *args, **kwargs):
|
||||
super().__init__(naxes=2, *args, **kwargs)
|
||||
|
||||
def add(self, freq, phasor, qtyname='G', **kwargs):
|
||||
L = 20*np.log10(np.abs(phasor))
|
||||
phase = np.angle(phasor)*180/np.pi
|
||||
self.semilogx(freq, L, axno=0, **kwargs)
|
||||
self.semilogx(freq, phase, axno=1, **kwargs)
|
||||
self.ylabel('$L$ [%s] [dB]' % qtyname, axno=0)
|
||||
self.ylabel(fr'$\angle$ {qtyname} [$^\circ$]', axno=1)
|
||||
self.xlabel('Frequency [Hz]', axno=1)
|
||||
|
||||
|
||||
class PS(Figure):
|
||||
def __init__(self, ref, *args, **kwargs):
|
||||
super().__init__(naxes=1, *args, **kwargs)
|
||||
self.ref = ref
|
||||
|
||||
def add(self, fs, freq, ps, qtyname='C', **kwargs):
|
||||
|
||||
overall = np.sum(ps)
|
||||
print(overall)
|
||||
overall_db = 10*np.log10(overall/self.ref**2)
|
||||
L = 10*np.log10(np.abs(ps)/self.ref**2)
|
||||
|
||||
self.semilogx(freq, L, **kwargs)
|
||||
# self.plot(freq,L,**kwargs)
|
||||
self.ylabel('Level [dB re 20$\\mu$Pa]')
|
||||
self.xlabel('Frequency [Hz]')
|
||||
self.legend('%s. Overall SPL = %0.1f dB SPL' % (qtyname, overall_db))
|
||||
|
||||
|
||||
class PSD(PS):
|
||||
def __init__(self, ref, *args, **kwargs):
|
||||
"""
|
||||
Initialize a PSD plot
|
||||
|
||||
Args:
|
||||
ref: Reference value for level in dB's
|
||||
|
||||
"""
|
||||
super().__init__(ref, *args, **kwargs)
|
||||
|
||||
def add(self, fs, freq, ps, qtyname='C', **kwargs):
|
||||
df = freq[1]-freq[0]
|
||||
nfft = fs/df
|
||||
df = fs/nfft
|
||||
psd = ps / df
|
||||
|
||||
overall = np.sum(np.abs(ps), axis=0)
|
||||
overall_db = 10*np.log10(overall/self.ref**2)
|
||||
L = 10*np.log10(abs(psd)/self.ref**2)
|
||||
|
||||
self.semilogx(freq, L, **kwargs)
|
||||
self.ylabel('$L$ [%s] [dB re %0.0e]' % (qtyname, self.ref))
|
||||
self.xlabel('Frequency [Hz]')
|
||||
self.legend('%s. Overall SPL = %0.1f dB SPL' % (qtyname, overall_db))
|
||||
@ -1,198 +0,0 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Author: J.A. de Jong - ASCEE
|
||||
|
||||
Description: backend tools for easy postprocessing of measurements
|
||||
"""
|
||||
from .plot import Figure
|
||||
from lasp.wrappers import AvPowerSpectra
|
||||
from lasp.lasp_measurement import Measurement
|
||||
from lasp.lasp_common import (FreqWeighting, TimeWeighting,
|
||||
getFreq, getTime, Window, P_REF)
|
||||
from lasp.lasp_weighcal import WeighCal
|
||||
from lasp.lasp_octavefilter import OctaveFilterBank, ThirdOctaveFilterBank
|
||||
from lasp.lasp_figuredialog import FigureDialog
|
||||
from lasp.lasp_figure import Plotable, PlotOptions
|
||||
from lasp.lasp_slm import SLM
|
||||
import numpy as np
|
||||
import sys
|
||||
|
||||
|
||||
def close():
|
||||
import matplotlib.pyplot as plt
|
||||
plt.close('all')
|
||||
|
||||
|
||||
def PSPlot(fn_list, **kwargs):
|
||||
"""
|
||||
Create a power spectral density plot, ASCEE style
|
||||
|
||||
Args:
|
||||
fn_list: list of measurement filenames to plot PSD for
|
||||
fw:
|
||||
fs:
|
||||
nfft:
|
||||
xscale:
|
||||
yscale:
|
||||
"""
|
||||
|
||||
fw = kwargs.pop('fw', FreqWeighting.A)
|
||||
nfft = kwargs.pop('nfft', 2048)
|
||||
xscale = kwargs.pop('xscale', 'log')
|
||||
yscale = kwargs.pop('yscale', 'PSD')
|
||||
ylim = kwargs.pop('ylim', (0, 100))
|
||||
xlim = kwargs.pop('xlim', (100, 10000))
|
||||
f = Figure(**kwargs)
|
||||
|
||||
print(kwargs)
|
||||
if xscale == 'log':
|
||||
pltfun = f.semilogx
|
||||
else:
|
||||
pltfun = f.plot
|
||||
|
||||
for fn in fn_list:
|
||||
meas = Measurement(fn)
|
||||
fs = meas.samplerate
|
||||
data = meas.praw()
|
||||
aps = AvPowerSpectra(nfft, 1, 50.)
|
||||
weighcal = WeighCal(fw, nchannels=1,
|
||||
fs=fs)
|
||||
weighted = weighcal.filter_(data)
|
||||
ps = aps.addTimeData(weighted)
|
||||
freq = getFreq(fs, nfft)
|
||||
if yscale == 'PSD':
|
||||
df = fs/nfft
|
||||
type_str = '/$\\sqrt{\\mathrm{Hz}}$'
|
||||
elif yscale == 'PS':
|
||||
df = 1.
|
||||
type_str = ''
|
||||
else:
|
||||
raise ValueError("'type' should be either 'PS' or 'PSD'")
|
||||
|
||||
psd_log = 10*np.log10(ps[:, 0, 0].real/df/2e-5**2)
|
||||
|
||||
pltfun(freq, psd_log)
|
||||
|
||||
f.xlabel('Frequency [Hz]')
|
||||
f.ylabel(f'Level [dB({fw[0]}) re (20$\\mu$ Pa){type_str}')
|
||||
f.ylim(ylim)
|
||||
f.xlim(xlim)
|
||||
return f
|
||||
|
||||
|
||||
def PowerSpectra(fn_list, **kwargs):
|
||||
nfft = kwargs.pop('nfft', 2048)
|
||||
window = kwargs.pop('window', Window.hann)
|
||||
fw = kwargs.pop('fw', FreqWeighting.A)
|
||||
overlap = kwargs.pop('overlap', 50.)
|
||||
ptas = []
|
||||
for fn in fn_list:
|
||||
meas = Measurement(fn)
|
||||
fs = meas.samplerate
|
||||
weighcal = WeighCal(fw, nchannels=1,
|
||||
fs=fs, calfile=None)
|
||||
praw = meas.praw()
|
||||
weighted = weighcal.filter_(praw)
|
||||
aps = AvPowerSpectra(nfft, 1, overlap, window[0])
|
||||
result = aps.addTimeData(weighted)[:, 0, 0].real
|
||||
pwr = 10*np.log10(result/P_REF**2)
|
||||
|
||||
freq = getFreq(fs, nfft)
|
||||
ptas.append(Plotable(freq, pwr))
|
||||
|
||||
pto = PlotOptions.forPower()
|
||||
pto.ylabel = f'L{fw[0]} [dB({fw[0]})]'
|
||||
return ptas
|
||||
|
||||
|
||||
def Levels(fn_list, **kwargs):
|
||||
|
||||
bank = kwargs.pop('bank', 'third')
|
||||
fw = kwargs.pop('fw', FreqWeighting.A)
|
||||
tw = kwargs.pop('tw', TimeWeighting.fast)
|
||||
xmin_txt = kwargs.pop('xmin', '100')
|
||||
xmax_txt = kwargs.pop('xmax', '16k')
|
||||
|
||||
levels = []
|
||||
leveltype = 'eq'
|
||||
|
||||
for fn in fn_list:
|
||||
meas = Measurement(fn)
|
||||
fs = meas.samplerate
|
||||
weighcal = WeighCal(fw, nchannels=1,
|
||||
fs=fs, calfile=None)
|
||||
praw = meas.praw()
|
||||
weighted = weighcal.filter_(praw)
|
||||
if bank == 'third':
|
||||
filt = ThirdOctaveFilterBank(fs)
|
||||
xmin = filt.nominal_txt_tox(xmin_txt)
|
||||
xmax = filt.nominal_txt_tox(xmax_txt)
|
||||
|
||||
elif bank == 'overall':
|
||||
slm = SLM(meas.samplerate, tw)
|
||||
slm.addData(weighted)
|
||||
levels.append(
|
||||
Plotable(' ', slm.Lmax if leveltype == 'max' else slm.Leq,
|
||||
name=meas.name))
|
||||
continue
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
# Octave bands
|
||||
# filt = OctaveFilterBank(fs)
|
||||
filtered_out = filt.filter_(weighted)
|
||||
level = np.empty((xmax - xmin + 1))
|
||||
xlabels = []
|
||||
for i, x in enumerate(range(xmin, xmax+1)):
|
||||
nom = filt.nominal_txt(x)
|
||||
xlabels.append(nom)
|
||||
filt_x = filtered_out[nom]['data']
|
||||
slm = SLM(filt.fs, tw)
|
||||
slm.addData(filt_x)
|
||||
leveli = slm.Lmax if leveltype == 'max' else slm.Leq
|
||||
level[i] = leveli
|
||||
levels.append(Plotable(xlabels, level, name=meas.name))
|
||||
return levels
|
||||
|
||||
|
||||
def LevelDifference(levels):
|
||||
assert len(levels) == 2
|
||||
return Plotable(name='Difference', x=levels[0].x,
|
||||
y=levels[1].y-levels[0].y)
|
||||
|
||||
|
||||
def LevelFigure(levels, show=True, **kwargs):
|
||||
figtype = kwargs.pop('figtype', 'bar')
|
||||
from PySide.QtGui import QApplication, QFont
|
||||
app = QApplication.instance()
|
||||
if not app:
|
||||
app = QApplication(sys.argv)
|
||||
app.setFont(QFont('Linux Libertine'))
|
||||
size = kwargs.pop('size', (1200, 600))
|
||||
if figtype == 'bar':
|
||||
opts = PlotOptions.forLevelBars()
|
||||
elif figtype == 'line':
|
||||
opts = PlotOptions.forPower()
|
||||
else:
|
||||
raise RuntimeError('figtype should be either line or bar')
|
||||
opts.ylim = kwargs.pop('ylim', (0, 100))
|
||||
opts.ylabel = kwargs.pop('ylabel', 'LAeq [dB(A)]')
|
||||
opts.xlabel = kwargs.pop('xlabel', None)
|
||||
opts.legend = kwargs.pop('legend', [level.name for level in levels])
|
||||
opts.legendpos = kwargs.pop('legendpos', None)
|
||||
opts.title = kwargs.pop('title', None)
|
||||
|
||||
def Plotter(ptas, pto):
|
||||
fig = FigureDialog(None, None, pto, figtype)
|
||||
for pta in ptas:
|
||||
fig.fig.add(pta)
|
||||
return fig
|
||||
|
||||
fig = Plotter(levels, opts)
|
||||
if show:
|
||||
fig.show()
|
||||
fig.resize(*size)
|
||||
if show:
|
||||
app.exec_()
|
||||
return fig
|
||||
@ -1,24 +1,58 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Created on Thu May 6 14:49:03 2021
|
||||
Author: C. Jansen, J.A. de Jong - ASCEE V.O.F.
|
||||
|
||||
@author: Casper
|
||||
|
||||
Smooth data in the frequency domain
|
||||
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.
|
||||
"""
|
||||
from enum import Enum, unique
|
||||
|
||||
__all__ = ['SmoothingType', 'smoothSpectralData', 'SmoothingWidth']
|
||||
|
||||
|
||||
@unique
|
||||
class SmoothingWidth(Enum):
|
||||
none = (0, 'No smoothing')
|
||||
# three = (3, '1/3th octave smoothing')
|
||||
six = (6, '1/6th octave smoothing')
|
||||
twelve = (12, '1/12th octave smoothing')
|
||||
twfo = (24, '1/24th octave smoothing')
|
||||
ftei = (48, '1/48th octave smoothing')
|
||||
|
||||
@staticmethod
|
||||
def fillComboBox(cb):
|
||||
"""
|
||||
Fill Windows to a combobox
|
||||
|
||||
Args:
|
||||
cb: QComboBox to fill
|
||||
"""
|
||||
cb.clear()
|
||||
for w in list(SmoothingWidth):
|
||||
cb.addItem(w.value[1], w)
|
||||
cb.setCurrentIndex(0)
|
||||
|
||||
@staticmethod
|
||||
def getCurrent(cb):
|
||||
return list(SmoothingWidth)[cb.currentIndex()]
|
||||
|
||||
class SmoothingType:
|
||||
levels = 'l', 'Levels'
|
||||
# tf = 'tf', 'Transfer function',
|
||||
ps = 'ps', '(Auto) powers'
|
||||
|
||||
|
||||
# TO DO: check if everything is correct
|
||||
# TO DO: add possibility to insert data that is not lin spaced in frequency
|
||||
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from scipy.signal.windows import gaussian
|
||||
|
||||
|
||||
# %% Smoothing function
|
||||
def oct_smooth(f, M, Noct, dB=False):
|
||||
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
|
||||
@ -30,49 +64,57 @@ def oct_smooth(f, M, Noct, dB=False):
|
||||
side. The deviation is largest when Noct is small (e.g. coarse smoothing).
|
||||
Casper Jansen, 07-05-2021
|
||||
|
||||
Parameters
|
||||
----------
|
||||
f : float
|
||||
frequencies of data points [Hz] - equally spaced
|
||||
M : float
|
||||
magnitude of data points [- or dB, specify in paramater 'dB']
|
||||
Noct : int
|
||||
smoothing strength: Noct=12 means 1/12 octave smoothing
|
||||
dB : Bool
|
||||
True if [M]=dB, False if [M]=absolute
|
||||
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
|
||||
-------
|
||||
f : float
|
||||
frequencies of data points [Hz]
|
||||
Msm : float
|
||||
smoothed magnitude of data points
|
||||
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
|
||||
|
||||
# Safety
|
||||
assert Noct > 0, '\'Noct\' must be absolute positive'
|
||||
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(f)==len(M), 'f and M should have equal length'
|
||||
if not dB: assert np.min(M) >= 0, 'absolute magnitude M cannot be negative'
|
||||
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 = len(M) # number of data points
|
||||
P = 10**(M/10) if dB else M**2 # convert magnitude --> power
|
||||
Psm = np.zeros(L) # smoothed power - to be calculated
|
||||
x0 = 1 if f[0]==0 else 0 # skip first data point if zero frequency
|
||||
df = f[1] - f[0] # frequency step
|
||||
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
|
||||
|
||||
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
|
||||
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 = f[x] # center freq. of smoothing window
|
||||
fc = freq[x] # center freq. of smoothing window
|
||||
Df = tr * fc / Noct # freq. range of smoothing window
|
||||
xl = int(np.ceil(x - 0.5*Df/df)) # desired lower index of frequency array to be used during smoothing
|
||||
xu = int(np.floor(x + 0.5*Df/df)) + 1 # upper index + 1 (because half-open interval)
|
||||
|
||||
# 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
|
||||
@ -110,13 +152,17 @@ def oct_smooth(f, M, Noct, dB=False):
|
||||
wind_int = np.sum(wind) # integral
|
||||
Psm[x] = np.dot(wind, P[xl:xu]) / wind_int # apply window
|
||||
|
||||
Msm = 10*np.log10(Psm) if dB else Psm**0.5 # convert power --> magnitude
|
||||
if st == SmoothingType.levels:
|
||||
Psm = 10*np.log10(Psm)
|
||||
elif st == SmoothingType.tf:
|
||||
Psm = np.sqrt(Psm)
|
||||
|
||||
return Msm
|
||||
return Psm
|
||||
|
||||
|
||||
# %% Test
|
||||
if __name__ == "__main__":
|
||||
import matplotlib.pyplot as plt
|
||||
# Initialize
|
||||
Noct = 6 # 1/6 oct. smoothing
|
||||
|
||||
|
||||
@ -74,8 +74,7 @@ cdef extern from "lasp_python.h":
|
||||
__all__ = ['AvPowerSpectra', 'SosFilterBank', 'FilterBank', 'Siggen',
|
||||
'sweep_flag_forward', 'sweep_flag_backward', 'sweep_flag_linear',
|
||||
'sweep_flag_exponential',
|
||||
'load_fft_wisdom', 'store_fft_wisdom', 'Window']
|
||||
|
||||
'load_fft_wisdom', 'store_fft_wisdom']
|
||||
|
||||
setTracerLevel(15)
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user