Lots of things improved. Forgot to count, became big commit, excuse myself

This commit is contained in:
Anne de Jong 2018-07-17 11:52:02 +02:00 committed by J.A. de Jong - ASCEE
parent 3311e43ef5
commit e95b02ae9a
13 changed files with 395 additions and 178 deletions

View File

@ -30,7 +30,7 @@ endif(LASP_FLOAT STREQUAL "double")
if(NOT DEFINED LASP_DEBUG) if(NOT DEFINED LASP_DEBUG)
message(SEND_ERROR "LASP_DEBUG flag not defined. Please set -DLASP_DEBUG=TRUE or -DLASP_DEBUG=FALSE") message(FATAL_ERROR "LASP_DEBUG flag not defined. Please set -DLASP_DEBUG=TRUE or -DLASP_DEBUG=FALSE")
endif(NOT DEFINED LASP_DEBUG) endif(NOT DEFINED LASP_DEBUG)
# ##################### END Cmake variables converted to a macro # ##################### END Cmake variables converted to a macro
@ -98,12 +98,6 @@ else()
set(win32 false) set(win32 false)
endif(CMAKE_SYSTEM_NAME STREQUAL "Windows") endif(CMAKE_SYSTEM_NAME STREQUAL "Windows")
# If numpy cannot be found in the standard include path of the Python
if(DEFINED NUMPY_INCLUDE)
include_directories(${NUMPY_INCLUDE})
endif(DEFINED NUMPY_INCLUDE)
add_subdirectory(fftpack) add_subdirectory(fftpack)
include_directories( include_directories(
fftpack fftpack
@ -112,8 +106,6 @@ include_directories(
add_subdirectory(lasp) add_subdirectory(lasp)
add_subdirectory(test) add_subdirectory(test)
find_program(PYTHON "python")
if (PYTHON) if (PYTHON)
# set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in") # set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in")
set(SETUP_PY "${CMAKE_CURRENT_BINARY_DIR}/setup.py") set(SETUP_PY "${CMAKE_CURRENT_BINARY_DIR}/setup.py")
@ -134,4 +126,3 @@ endif()
############################## End compiler settings ############################## End compiler settings

BIN
img/LASP.pdf Normal file

Binary file not shown.

View File

@ -1,16 +1,15 @@
add_subdirectory(c) add_subdirectory(c)
configure_file(config.pxi.in config.pxi) configure_file(config.pxi.in config.pxi)
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
find_package( PythonLibs REQUIRED ) find_package(PythonLibs REQUIRED )
include(UseCython) include(UseCython)
find_package(Numpy REQUIRED )
include_directories( include_directories(
. ${PYTHON_NUMPY_INCLUDE_DIR}
c .
) c
)
# add the command to generate the source code # add the command to generate the source code
# add_custom_command ( # add_custom_command (

View File

@ -0,0 +1,41 @@
# Find the Python NumPy package
# PYTHON_NUMPY_INCLUDE_DIR
# PYTHON_NUMPY_FOUND
# will be set by this script
cmake_minimum_required(VERSION 2.6)
if(NOT PYTHON_EXECUTABLE)
if(NumPy_FIND_QUIETLY)
find_package(PythonInterp QUIET)
else()
find_package(PythonInterp)
set(__numpy_out 1)
endif()
endif()
if (PYTHON_EXECUTABLE)
# Find out the include path
execute_process(
COMMAND "${PYTHON_EXECUTABLE}" -c
"from __future__ import print_function\ntry: import numpy; print(numpy.get_include(), end='')\nexcept:pass\n"
OUTPUT_VARIABLE __numpy_path)
# And the version
execute_process(
COMMAND "${PYTHON_EXECUTABLE}" -c
"from __future__ import print_function\ntry: import numpy; print(numpy.__version__, end='')\nexcept:pass\n"
OUTPUT_VARIABLE __numpy_version)
elseif(__numpy_out)
message(STATUS "Python executable not found.")
endif(PYTHON_EXECUTABLE)
find_path(PYTHON_NUMPY_INCLUDE_DIR numpy/arrayobject.h
HINTS "${__numpy_path}" "${PYTHON_INCLUDE_PATH}" NO_DEFAULT_PATH)
if(PYTHON_NUMPY_INCLUDE_DIR)
set(PYTHON_NUMPY_FOUND 1 CACHE INTERNAL "Python numpy found")
endif(PYTHON_NUMPY_INCLUDE_DIR)
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(NumPy REQUIRED_VARS PYTHON_NUMPY_INCLUDE_DIR
VERSION_VAR __numpy_version)

View File

@ -175,6 +175,8 @@ class ThirdOctaveBankDesigner(FilterBankDesigner):
'6.3k', '8k', '10k', '6.3k', '8k', '10k',
'12.5k', '16k', '20k'] '12.5k', '16k', '20k']
assert len(self.xs) == len(self.nominal_txt)
@property @property
def b(self): def b(self):
# Band division factor, 3 for one-third octave bands # Band division factor, 3 for one-third octave bands

View File

@ -35,6 +35,8 @@ The video dataset can possibly be not present in the data.
""" """
__all__ = ['Measurement']
from contextlib import contextmanager
import h5py as h5 import h5py as h5
import numpy as np import numpy as np
from .lasp_config import LASP_NUMPY_FLOAT_TYPE from .lasp_config import LASP_NUMPY_FLOAT_TYPE
@ -47,7 +49,7 @@ class BlockIter:
Iterate over the blocks in the audio data of a h5 file Iterate over the blocks in the audio data of a h5 file
""" """
def __init__(self, faudio): def __init__(self, f):
""" """
Initialize a BlockIter object Initialize a BlockIter object
@ -55,8 +57,8 @@ class BlockIter:
faudio: Audio dataset in the h5 file, accessed as f['audio'] faudio: Audio dataset in the h5 file, accessed as f['audio']
""" """
self.i = 0 self.i = 0
self.nblocks = faudio.shape[0] self.nblocks = f['audio'].shape[0]
self.fa = faudio self.fa = f['audio']
def __iter__(self): def __iter__(self):
return self return self
@ -121,35 +123,55 @@ class Measurement:
# Full filepath # Full filepath
self.fn = fn self.fn = fn
# Base filename # Base filename
self.fn_base = os.path.split(fn)[1] self.fn_base = os.path.split(fn)[1]
# Open the h5 file in read-plus mode, to allow for changing the # Open the h5 file in read-plus mode, to allow for changing the
# measurement comment. # measurement comment.
f = h5.File(fn, 'r+') with h5.File(fn, 'r+') as f:
self.f = f # Check for video data
try:
f['video']
self.has_video = True
except KeyError:
self.has_video = False
# Check for video data self.nblocks, self.blocksize, self.nchannels = f['audio'].shape
try: dtype = f['audio'].dtype
f['video'] self.sampwidth = getSampWidth(dtype)
self.has_video = True
except KeyError:
self.has_video = False
self.nblocks, self.blocksize, self.nchannels = f['audio'].shape self.samplerate = f.attrs['samplerate']
dtype = f['audio'].dtype self.N = (self.nblocks*self.blocksize)
self.sampwidth = getSampWidth(dtype) self.T = self.N/self.samplerate
self.samplerate = f.attrs['samplerate'] # comment = read-write thing
self.N = (self.nblocks*self.blocksize) try:
self.T = self.N/self.samplerate self._comment = f.attrs['comment']
except KeyError:
f.attrs['comment'] = ''
self._comment = ''
# comment = read-write thing # Sensitivity
try: try:
self._comment = self.f.attrs['comment'] self._sens = f.attrs['sensitivity']
except KeyError: except KeyError:
self.f.attrs['comment'] = '' self._sens = np.ones(self.nchannels)
self._comment = ''
self._time = f.attrs['time']
@contextmanager
def file(self, mode='r'):
"""
Contextmanager which opens the storage file and yields the file.
Args:
mode: Opening mode for the file. Should either be 'r', or 'r+'
"""
if mode not in ('r','r+'):
raise ValueError('Invalid file opening mode.')
with h5.File(self.fn, mode) as f:
yield f
@property @property
def comment(self): def comment(self):
@ -157,19 +179,23 @@ class Measurement:
@comment.setter @comment.setter
def comment(self, cmt): def comment(self, cmt):
self.f.attrs['comment'] = cmt with self.file('r+') as f:
self._comment = cmt f.attrs['comment'] = cmt
self._comment = cmt
@property @property
def recTime(self): def recTime(self):
return (self.blocksize*self.nblocks)/self.samplerate """
Returns the total recording time of the measurement, in float seconds.
"""
return self.blocksize*self.nblocks/self.samplerate
@property @property
def time(self): def time(self):
""" """
Returns the measurement time in seconds since the epoch. Returns the measurement time in seconds since the epoch.
""" """
return self.f.attrs['time'] return self._time
@property @property
def prms(self): def prms(self):
@ -186,7 +212,7 @@ class Measurement:
except AttributeError: except AttributeError:
pass pass
sens = self.sensitivity sens = self._sens
pms = 0. pms = 0.
for block in self.iterBlocks(): for block in self.iterBlocks():
pms += np.sum(block/sens[np.newaxis, :], axis=0)**2/self.N pms += np.sum(block/sens[np.newaxis, :], axis=0)**2/self.N
@ -198,13 +224,14 @@ class Measurement:
Returns the raw uncalibrated data, converted to floating point format. Returns the raw uncalibrated data, converted to floating point format.
""" """
if block is not None: if block is not None:
blocks = self.f['audio'][block] with self.file() as f:
blocks = f['audio'][block]
else: else:
blocks = [] blocks = []
for block in self.iterBlocks(): with self.file() as f:
blocks.append(block) for block in self.iterBlocks(f):
blocks.append(block)
blocks = np.asarray(blocks) blocks = np.asarray(blocks)
blocks = blocks.reshape(self.nblocks*self.blocksize, blocks = blocks.reshape(self.nblocks*self.blocksize,
self.nchannels) self.nchannels)
@ -220,13 +247,19 @@ class Measurement:
else: else:
raise RuntimeError( raise RuntimeError(
f'Unimplemented data type from recording: {blocks.dtype}.') f'Unimplemented data type from recording: {blocks.dtype}.')
sens = self.sensitivity sens = self._sens
blocks = blocks.astype(LASP_NUMPY_FLOAT_TYPE)/fac/sens[np.newaxis, :] blocks = blocks.astype(LASP_NUMPY_FLOAT_TYPE)/fac/sens[np.newaxis, :]
return blocks return blocks
def iterBlocks(self): def iterBlocks(self, opened_file):
return BlockIter(self.f['audio']) """
Iterate over all the audio blocks in the opened file
Args:
opened_file: The h5File with the data
"""
return BlockIter(opened_file)
@property @property
def sensitivity(self): def sensitivity(self):
@ -235,10 +268,7 @@ class Measurement:
between -1.0 and 1.0 to Pascal. If the sensitivity is not stored in between -1.0 and 1.0 to Pascal. If the sensitivity is not stored in
the measurement file, this function returns 1.0 the measurement file, this function returns 1.0
""" """
try: return self._sens
return self.f.attrs['sensitivity']
except KeyError:
return np.ones(self.nchannels)
@sensitivity.setter @sensitivity.setter
def sensitivity(self, sens): def sensitivity(self, sens):
@ -257,8 +287,8 @@ class Measurement:
valid &= isinstance(sens.dtype, float) valid &= isinstance(sens.dtype, float)
if not valid: if not valid:
raise ValueError('Invalid sensitivity value(s) given') raise ValueError('Invalid sensitivity value(s) given')
with self.file('r+') as f:
self.f.attrs['sensitivity'] = sens f.attrs['sensitivity'] = sens
def exportAsWave(self, fn=None, force=False, sampwidth=None): def exportAsWave(self, fn=None, force=False, sampwidth=None):
""" """
@ -286,8 +316,8 @@ class Measurement:
if os.path.exists(fn) and not force: if os.path.exists(fn) and not force:
raise RuntimeError(f'File already exists: {fn}') raise RuntimeError(f'File already exists: {fn}')
with self.file() as f:
audio = self.f['audio'] audio = self.f['audio'][:]
if isinstance(audio.dtype, float): if isinstance(audio.dtype, float):
if sampwidth is None: if sampwidth is None:
@ -373,8 +403,5 @@ class Measurement:
ad[0] = dat ad[0] = dat
return Measurement(mfn) return Measurement(mfn)
def __del__(self): # def __del__(self):
try: # self.f.close()
self.f.close()
except AttributeError:
pass

View File

@ -17,21 +17,19 @@ class FilterBank:
Single channel octave filter bank implementation Single channel octave filter bank implementation
""" """
def __init__(self, fs, designer): def __init__(self, fs):
""" """
Initialize a OctaveFilterBank object. Initialize a OctaveFilterBank object.
Args: Args:
fs: Sampling frequency of base signal fs: Sampling frequency of base signal
designer: FIR Filter designer for filterbank
""" """
assert np.isclose(fs, 48000), "Only sampling frequency" \ assert np.isclose(fs, 48000), "Only sampling frequency" \
" available is 48 kHz" " available is 48 kHz"
self.fs = fs
maxdecimation = designer.decimation(designer.xs[0]) maxdecimation = self.decimation(self.xs[0])
self.decimators = [] self.decimators = []
for dec in maxdecimation: for dec in maxdecimation:
self.decimators.append(Decimator(1, dec)) self.decimators.append(Decimator(1, dec))
@ -44,8 +42,8 @@ class FilterBank:
self.filterbanks = [] self.filterbanks = []
# Sort the x values in categories according to the required decimation # Sort the x values in categories according to the required decimation
for x in designer.xs: for x in self.xs:
dec = designer.decimation(x) dec = self.decimation(x)
if len(dec) == 1 and dec[0] == 1: if len(dec) == 1 and dec[0] == 1:
xs_d1.append(x) xs_d1.append(x)
elif len(dec) == 1 and dec[0] == 4: elif len(dec) == 1 and dec[0] == 4:
@ -62,12 +60,12 @@ class FilterBank:
xs_all = [xs_d1, xs_d4, xs_d16, xs_d64, xs_d256] xs_all = [xs_d1, xs_d4, xs_d16, xs_d64, xs_d256]
for xs in xs_all: for xs in xs_all:
nominals = [] nominals = []
firs = np.empty((designer.L, len(xs)), order='F') firs = np.empty((self.L, len(xs)), order='F')
for i, x in enumerate(xs): for i, x in enumerate(xs):
# These are the filters that do not require lasp_decimation # These are the filters that do not require lasp_decimation
# prior to filtering # prior to filtering
nominals.append(designer.nominal(x)) nominals.append(self.nominal(x))
firs[:, i] = designer.createFilter(fs, x) firs[:, i] = self.createFilter(fs, x)
filterbank = {'fb': pyxFilterBank(firs, 1024), filterbank = {'fb': pyxFilterBank(firs, 1024),
'xs': xs, 'xs': xs,
'nominals': nominals} 'nominals': nominals}
@ -78,16 +76,10 @@ class FilterBank:
# Filter output counters # Filter output counters
self.dec = [1, 4, 16, 64, 256] self.dec = [1, 4, 16, 64, 256]
Pdec = 128
Pfil = 256
# Initial filtered number # These intial delays are found experimentally using a toneburst
self.Nf = [-Pfil, -(Pdec/4 + Pfil), # response.
-(Pfil + Pdec/4 + Pdec/16), self.Nf = [915, 806, 780, 582, 338]
-(Pfil + Pdec/64 + Pdec/4 + Pdec/16),
-(Pfil + Pdec/256 + Pdec/64 + Pdec/4 + Pdec/16)]
self.delay_set = [False, False, False, False, False]
def filterd(self, dec_stage, data): def filterd(self, dec_stage, data):
""" """
@ -110,10 +102,12 @@ class FilterBank:
oldNf = self.Nf[dec_stage] oldNf = self.Nf[dec_stage]
tstart = oldNf/fd tstart = oldNf/fd
tend = tstart + Nf/fd tend = tstart + Nf/fd
t = np.linspace(tstart, tend, Nf, endpoint=False) t = np.linspace(tstart, tend, Nf, endpoint=False)
self.Nf[dec_stage] += Nf self.Nf[dec_stage] += Nf
for i, nom in enumerate(self.filterbanks[dec_stage]['nominals']): for i, nom in enumerate(self.filterbanks[dec_stage]['nominals']):
output[nom] = {'t': t, 'data': filtered[:, i]} output[nom] = {'t': t, 'data': filtered[:, [i]]}
return output return output
def filter_(self, data): def filter_(self, data):
@ -142,13 +136,13 @@ class FilterBank:
return output return output
class OctaveFilterBank(FilterBank): class OctaveFilterBank(FilterBank, OctaveBankDesigner):
def __init__(self, fs): def __init__(self, fs):
designer = OctaveBankDesigner() OctaveBankDesigner.__init__(self)
super().__init__(fs, designer) FilterBank.__init__(self, fs)
class ThirdOctaveFilterBank(FilterBank): class ThirdOctaveFilterBank(FilterBank, ThirdOctaveBankDesigner):
def __init__(self, fs): def __init__(self, fs):
designer = ThirdOctaveBankDesigner() ThirdOctaveBankDesigner.__init__(self)
super().__init__(fs, designer) FilterBank.__init__(self, fs)

View File

@ -24,7 +24,8 @@ class ReverbTime:
dB. dB.
""" """
assert level.ndim == 1, 'Invalid number of dimensions in level' assert level.ndim == 2, 'Invalid number of dimensions in level'
assert level.shape[1] == 1,'Invalid number of channels, should be 1'
self._level = level self._level = level
# Number of time samples # Number of time samples
self._N = self._level.shape[0] self._N = self._level.shape[0]

View File

@ -14,7 +14,8 @@ from .lasp_weighcal import WeighCal
from .lasp_gui_tools import wait_cursor from .lasp_gui_tools import wait_cursor
from .lasp_figure import PlotOptions, Plotable from .lasp_figure import PlotOptions, Plotable
from .ui_slmwidget import Ui_SlmWidget from .ui_slmwidget import Ui_SlmWidget
from .filter.bandpass_fir import OctaveBankDesigner, ThirdOctaveBankDesigner
from .lasp_octavefilter import OctaveFilterBank, ThirdOctaveFilterBank
__all__ = ['SLM', 'SlmWidget'] __all__ = ['SLM', 'SlmWidget']
@ -29,30 +30,39 @@ class Dummy:
class SLM: class SLM:
""" """
Sound Level Meter, implements the single pole lowpass filter Multi-channel sound Level Meter. Input data: time data with a certain
sampling frequency. Output: time-weighted (fast/slow) sound pressure
levels in dB(A/C/Z).
""" """
def __init__(self, fs, weighcal, def __init__(self, fs, tw=TimeWeighting.default):
tw=TimeWeighting.default,
):
""" """
Initialize a sound level meter object. Number of channels comes from Initialize a sound level meter object.
the weighcal object Number of channels comes from the weighcal object.
Args: Args:
fs: Sampling frequency [Hz] fs: Sampling frequency [Hz]
weighcal: WeighCal instance used for calibration and frequency tw: Time Weighting to apply
weighting.
nchannels: Number of channels to allocate filters for
""" """
nchannels = weighcal.nchannels
self.nchannels = nchannels
self._weighcal = weighcal
if tw[0] is not TimeWeighting.none[0]: if tw[0] is not TimeWeighting.none[0]:
self._lps = [SPLowpass(fs, tw[0]) for i in range(nchannels)] self._lp = SPLowpass(fs, tw[0])
else: else:
self._lps = [Dummy() for i in range(nchannels)] self._lp = Dummy()
self._Lmax = zeros(nchannels) self._Lmax = 0.
# Storage for computing the equivalent level
self._sq = 0.
self._N = 0
self._Leq = 0.
@property
def Leq(self):
"""
Returns the equivalent level of the recorded signal so far
"""
return self._Leq
@property @property
def Lmax(self): def Lmax(self):
@ -68,26 +78,26 @@ class SLM:
Args: Args:
data: data:
""" """
if data.ndim == 1: assert data.ndim == 2
data = data[:, np.newaxis] assert data.shape[1] == 1
data_weighted = self._weighcal.filter_(data) P_REFsq = P_REF**2
# Squared # Squared
sq = data_weighted**2 sq = data**2
if sq.shape[0] == 0:
return np.array([])
tw = [] # Update equivalent level
N1 = sq.shape[0]
self._sq = (np.sum(sq) + self._sq*self._N)/(N1+self._N)
self._N += N1
self._Leq = 10*np.log10(self._sq/P_REFsq)
# Time-weight the signal # Time-weight the signal
for chan, lp in enumerate(self._lps): tw = self._lp.filter_(sq)
tw.append(lp.filter_(sq[:, chan])[:, 0])
tw = np.asarray(tw).transpose() Level = 10*np.log10(tw/P_REFsq)
Level = 10*np.log10(tw/P_REF**2)
# Update maximum level
curmax = np.max(Level) curmax = np.max(Level)
if curmax > self._Lmax: if curmax > self._Lmax:
self._Lmax = curmax self._Lmax = curmax
@ -103,8 +113,8 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget):
super().__init__(parent) super().__init__(parent)
self.setupUi(self) self.setupUi(self)
FreqWeighting.fillComboBox(self.tfreqweighting) self.eqFreqBandChanged(0)
FreqWeighting.fillComboBox(self.eqfreqweighting) self.tFreqBandChanged(0)
self.setMeas(None) self.setMeas(None)
def init(self, fm): def init(self, fm):
@ -115,6 +125,9 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget):
fm.registerCombo(self.tfigure) fm.registerCombo(self.tfigure)
fm.registerCombo(self.eqfigure) fm.registerCombo(self.eqfigure)
self.tbandstart.setEnabled(False)
self.tbandstop.setEnabled(False)
def setMeas(self, meas): def setMeas(self, meas):
""" """
Set the current measurement for this widget. Set the current measurement for this widget.
@ -151,16 +164,33 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget):
channel = self.eqchannel.currentIndex() channel = self.eqchannel.currentIndex()
fw = FreqWeighting.getCurrent(self.eqfreqweighting) fw = FreqWeighting.getCurrent(self.eqfreqweighting)
startpos = self.eqstarttime.value istart, istop = self.getStartStopIndices(meas, self.eqstarttime,
stoppos = self.eqstoptime.value self.eqstoptime)
N = meas.N
istart = int(startpos*fs)
if istart >= N:
raise ValueError("Invalid start position")
istop = int(stoppos*fs)
if istart > N:
raise ValueError("Invalid stop position")
bands = self.eqfreqband.currentIndex()
if bands == 0:
# 1/3 Octave bands
filt = ThirdOctaveFilterBank(fs)
xs = filt.xs
xmin = xs[0] + self.eqbandstart.currentIndex()
xmax = xs[0] + self.eqbandstop.currentIndex()
if bands == 1:
# Octave bands
filt = OctaveFilterBank(fs)
xs = filt.xs
xmin = xs[0] + self.eqbandstart.currentIndex()
xmax = xs[0] + self.eqbandstop.currentIndex()
leveltype = self.eqleveltype.currentIndex()
if leveltype == 0:
# equivalent levels
tw = TimeWeighting.fast
elif leveltype == 1:
# fast time weighting
tw = TimeWeighting.fast
elif leveltype == 2:
# slow time weighting
tw = TimeWeighting.slow
with wait_cursor(): with wait_cursor():
# This one exctracts the calfile and sensitivity from global # This one exctracts the calfile and sensitivity from global
@ -170,8 +200,28 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget):
fs=fs, calfile=calfile, fs=fs, calfile=calfile,
sens=sens) sens=sens)
praw = meas.praw()[istart:istop, [channel]] praw = meas.praw()[istart:istop, [channel]]
pto = PlotOptions()
weighted = weighcal.filter_(praw)
filtered_out = filt.filter_(weighted)
levels = np.empty((xmax - xmin + 1))
xlabels = []
for i, x in enumerate(range(xmin, xmax+1)):
nom = filt.nominal(x)
xlabels.append(nom)
filt_x = filtered_out[nom]['data']
slm = SLM(filt.fs, tw)
slm.addData(filt_x)
if leveltype > 0:
level = slm.Lmax
else:
level = slm.Leq
levels[i] = level
pto = PlotOptions.forLevelBars()
pta = Plotable(xlabels, levels)
fig, new = self.getFigure(self.eqfigure, pto, 'bar') fig, new = self.getFigure(self.eqfigure, pto, 'bar')
fig.fig.add(pta)
fig.show() fig.show()
def computeT(self): def computeT(self):
@ -184,46 +234,95 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget):
tw = TimeWeighting.getCurrent(self.ttimeweighting) tw = TimeWeighting.getCurrent(self.ttimeweighting)
fw = FreqWeighting.getCurrent(self.tfreqweighting) fw = FreqWeighting.getCurrent(self.tfreqweighting)
istart, istop = self.getStartStopIndices(meas, self.tstarttime,
self.tstoptime)
bands = self.tfreqband.currentIndex()
if bands == 0:
# Overall
filt = Dummy()
else:
# Octave bands
filt = OctaveFilterBank(
fs) if bands == 1 else ThirdOctaveFilterBank(fs)
xs = filt.xs
xmin = xs[0] + self.tbandstart.currentIndex()
xmax = xs[0] + self.tbandstop.currentIndex()
# Downsampling factor of result # Downsampling factor of result
dsf = self.tdownsampling.value() dsf = self.tdownsampling.value()
# gb = self.slmFre
with wait_cursor(): with wait_cursor():
# This one exctracts the calfile and sensitivity from global # This one exctracts the calfile and sensitivity from global
# variables defined at the top. # TODO: Change this to a more # variables defined at the top. # TODO: Change this to a more
# robust variant. # robust variant.
praw = meas.praw()[istart:istop, [channel]]
weighcal = WeighCal(fw, nchannels=1, weighcal = WeighCal(fw, nchannels=1,
fs=fs, calfile=calfile, fs=fs, calfile=calfile,
sens=sens) sens=sens)
slm = SLM(fs, weighcal, tw) weighted = weighcal.filter_(praw)
praw = meas.praw()[:, [channel]]
# Filter, downsample data if bands == 0:
filtered = slm.addData(praw)[::dsf, :] slm = SLM(fs, tw)
N = filtered.shape[0] level = slm.addData(weighted)[::dsf]
time = getTime(float(fs)/dsf, N)
Lmax = slm.Lmax
pto = PlotOptions() # Filter, downsample data
pto.ylabel = f'L{fw[0]} [dB({fw[0]})]' N = level.shape[0]
pto.xlim = (time[0], time[-1]) time = getTime(float(fs)/dsf, N)
Lmax = slm.Lmax
pta = Plotable(time, filtered) pta = Plotable(time, level,
name=f'Overall level [dB([fw[0]])]')
pto = PlotOptions()
pto.ylabel = f'L{fw[0]} [dB({fw[0]})]'
pto.xlim = (time[0], time[-1])
fig, new = self.getFigure(self.tfigure, pto, 'line')
fig.fig.add(pta)
fig, new = self.getFigure(self.tfigure, pto, 'line') else:
fig.fig.add(pta) pto = PlotOptions()
fig, new = self.getFigure(self.tfigure, pto, 'line')
pto.ylabel = f'L{fw[0]} [dB({fw[0]})]'
out = filt.filter_(weighted)
tmin = 0
tmax = 0
for x in range(xmin, xmax+1):
dec = np.prod(filt.decimation(x))
fd = filt.fs/dec
# Nominal frequency text
nom = filt.nominal(x)
leg = f'{nom} Hz - [dB({fw[0]})]'
# Find global tmin and tmax, used for xlim
time = out[nom]['t']
tmin = min(tmin, time[0])
tmax = max(tmax, time[-1])
slm = SLM(fd, tw)
level = slm.addData(out[nom]['data'])
plotable = Plotable(time[::dsf//dec],
level[::dsf//dec],
name=leg)
fig.fig.add(plotable)
pto.xlim = (tmin, tmax)
fig.fig.setPlotOptions(pto)
fig.show() fig.show()
stats = f"""Statistical results: # stats = f"""Statistical results:
============================= # =============================
Applied frequency weighting: {fw[1]} # Applied frequency weighting: {fw[1]}
Applied time weighting: {tw[1]} # Applied time weighting: {tw[1]}
Applied Downsampling factor: {dsf} # Applied Downsampling factor: {dsf}
Maximum level (L{fw[0]} max): {Lmax:4.4} [dB({fw[0]})] # Maximum level (L{fw[0]} max): {Lmax:4.4} [dB({fw[0]})]
#
""" # """
self.results.setPlainText(stats) # self.results.setPlainText(stats)
def compute(self): def compute(self):
""" """
@ -234,3 +333,61 @@ Maximum level (L{fw[0]} max): {Lmax:4.4} [dB({fw[0]})]
self.computeT() self.computeT()
elif self.eqtab.isVisible(): elif self.eqtab.isVisible():
self.computeEq() self.computeEq()
def eqFreqBandChanged(self, idx):
"""
User changes frequency bands to plot time-dependent values for
"""
self.eqbandstart.clear()
self.eqbandstop.clear()
if idx == 1:
# 1/3 Octave bands
o = OctaveBankDesigner()
for x in o.xs:
nom = o.nominal(x)
self.eqbandstart.addItem(nom)
self.eqbandstop.addItem(nom)
self.eqbandstart.setCurrentIndex(0)
self.eqbandstop.setCurrentIndex(len(o.xs)-1)
elif idx == 0:
# Octave bands
o = ThirdOctaveBankDesigner()
for x in o.xs:
nom = o.nominal(x)
self.eqbandstart.addItem(nom)
self.eqbandstop.addItem(nom)
self.eqbandstart.setCurrentIndex(2)
self.eqbandstop.setCurrentIndex(len(o.xs) - 3)
def tFreqBandChanged(self, idx):
"""
User changes frequency bands to plot time-dependent values for
"""
self.tbandstart.clear()
self.tbandstop.clear()
enabled = False
if idx == 1:
# Octave bands
enabled = True
o = OctaveBankDesigner()
for x in o.xs:
nom = o.nominal(x)
self.tbandstart.addItem(nom)
self.tbandstop.addItem(nom)
self.tbandstart.setCurrentIndex(2)
self.tbandstop.setCurrentIndex(len(o.xs)-1)
elif idx == 2:
# Octave bands
enabled = True
o = ThirdOctaveBankDesigner()
for x in o.xs:
nom = o.nominal(x)
self.tbandstart.addItem(nom)
self.tbandstop.addItem(nom)
self.tbandstart.setCurrentIndex(2)
self.tbandstop.setCurrentIndex(len(o.xs) - 3)
self.tbandstart.setEnabled(enabled)
self.tbandstop.setEnabled(enabled)

View File

@ -4,10 +4,10 @@
Weighting and calibration filter in one Weighting and calibration filter in one
@author: J.A. de Jong - ASCEE @author: J.A. de Jong - ASCEE
""" """
from .filter_design.freqweighting_fir import A,C from .filter.freqweighting_fir import A, C
from .lasp_common import FreqWeighting from .lasp_common import FreqWeighting
from .filter_design.fir_design import (arbitrary_fir_design, from .filter.fir_design import (arbitrary_fir_design,
freqResponse as frp) freqResponse as frp)
from lasp.lasp_config import ones, empty from lasp.lasp_config import ones, empty
from .wrappers import FilterBank from .wrappers import FilterBank
import numpy as np import numpy as np
@ -17,7 +17,7 @@ __all__ = ['WeighCal']
class WeighCal: class WeighCal:
""" """
Time weighting and calibration FIR filter Frequency weighting and calibration FIR filter
""" """
def __init__(self, fw=FreqWeighting.default, def __init__(self, fw=FreqWeighting.default,
nchannels=1, nchannels=1,
@ -48,15 +48,17 @@ class WeighCal:
# Objective function for the frequency response # Objective function for the frequency response
frp_obj = self.frpObj(freq_design) frp_obj = self.frpObj(freq_design)
self._firs = [] P = 2048 # Filter length (number of taps)
self._firs = np.empty((P, self.nchannels))
self._fbs = [] self._fbs = []
for chan in range(self.nchannels): for chan in range(self.nchannels):
fir = arbitrary_fir_design(fs, 2048, freq_design, fir = arbitrary_fir_design(fs, P, freq_design,
frp_obj[:, chan], frp_obj[:, chan],
window='rectangular') window='rectangular')
self._firs.append(fir) self._firs[:, chan] = fir
self._fbs.append(FilterBank(fir[:, np.newaxis], 4096)) self._fbs.append(FilterBank(fir[:, np.newaxis], 2*P))
self._freq_design = freq_design self._freq_design = freq_design
@ -66,22 +68,21 @@ class WeighCal:
Args: Args:
data: (Weighted) raw time data that needs to be filtered, should data: (Weighted) raw time data that needs to be filtered, should
have the same number of columns as the number of channels, or should have the same number of columns as the number of channels. First
have dimension 1 in case of a single channel. axis is assumed to be the time axis
Retuns: Retuns:
Filtered data for each channel Filtered data for each channel
""" """
nchan = self.nchannels nchan = self.nchannels
if data.ndim == 1: assert data.ndim == 2
data = data[:, np.newaxis]
assert data.shape[1] == nchan assert data.shape[1] == nchan
assert data.shape[0] > 0 assert data.shape[0] > 0
filtered = [] filtered = []
for chan in range(nchan): for chan in range(nchan):
filtered.append(self._fbs[chan].filter_(data[:, chan])[:, 0]) filtered.append(self._fbs[chan].filter_(data[:, [chan]])[:, 0])
filtered = np.asarray(filtered).transpose()/self.sens filtered = np.asarray(filtered).transpose()/self.sens
if filtered.ndim == 1: if filtered.ndim == 1:
filtered = filtered[:, np.newaxis] filtered = filtered[:, np.newaxis]
@ -133,7 +134,7 @@ class WeighCal:
""" """
# Objective function for the frequency response # Objective function for the frequency response
frp_objective = self.frpCalObj(freq_design) * \ frp_objective = self.frpCalObj(freq_design) * \
self.frpWeightingObj(freq_design)[:, np.newaxis] self.frpWeightingObj(freq_design)[:, np.newaxis]
frp_objective[-1] = 0. frp_objective[-1] = 0.
return frp_objective return frp_objective

View File

@ -55,10 +55,15 @@ class BarScene(QGraphicsScene):
xvals: labels and x positions of the bars xvals: labels and x positions of the bars
G: Number of bars per x value G: Number of bars per x value
ylim: y limits of the figure ylim: y limits of the figure
xlabel: label below x-axis
ylabel: label on left side of the y-axis
title: figure title
colors: color cycler
size: size of the plot in pixels
legend: list of legend strings to show.
""" """
super().__init__(parent=parent) super().__init__(parent=parent)
self.setSceneRect(QRect(0,0,*size))
# self.setBackgroundBrush(ASCEEColors.bgBrush(0, size[0])) # self.setBackgroundBrush(ASCEEColors.bgBrush(0, size[0]))
self.ylim = ylim self.ylim = ylim
@ -68,8 +73,7 @@ class BarScene(QGraphicsScene):
self.bgs = [] self.bgs = []
self.size = size self.size = size
xsize = size[0] xsize, ysize = size
ysize = size[1]
self.xsize = xsize self.xsize = xsize
self.ysize = ysize self.ysize = ysize

View File

@ -377,12 +377,13 @@ cdef class SPLowpass:
if self.lp: if self.lp:
SPLowpass_free(self.lp) SPLowpass_free(self.lp)
def filter_(self,d[:] input_): def filter_(self,d[::1,:] input_):
assert input_.shape[1] == 1
if input_.shape[0] == 0: if input_.shape[0] == 0:
return np.array([],dtype=NUMPY_FLOAT_TYPE) return np.array([],dtype=NUMPY_FLOAT_TYPE)
cdef vd input_vd = dmat_foreign_data(input_.shape[0],1, cdef vd input_vd = dmat_foreign_data(input_.shape[0],1,
&input_[0],False) &input_[0,0],False)
cdef dmat output = SPLowpass_filter(self.lp,&input_vd) cdef dmat output = SPLowpass_filter(self.lp,&input_vd)

View File

@ -16,9 +16,8 @@ setup(
author_email="j.a.dejong@ascee.nl", author_email="j.a.dejong@ascee.nl",
# Project uses reStructuredText, so ensure that the docutils get # Project uses reStructuredText, so ensure that the docutils get
# installed or upgraded on the target machine # installed or upgraded on the target machine
install_requires=['matplotlib>=1.0', 'scipy>=0.17', 'numpy>=1.0'], install_requires=['matplotlib>=1.0', 'scipy>=1.0', 'numpy>=1.0'],
license='MIT', license='MIT',
description=descr, description=descr,
keywords="", keywords="",
url="http://www.ascee.nl/lasp/", # project home page, if any url="http://www.ascee.nl/lasp/", # project home page, if any