diff --git a/CMakeLists.txt b/CMakeLists.txt index 8805377..b9aaf31 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -30,7 +30,7 @@ endif(LASP_FLOAT STREQUAL "double") 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) # ##################### END Cmake variables converted to a macro @@ -98,12 +98,6 @@ else() set(win32 false) 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) include_directories( fftpack @@ -112,8 +106,6 @@ include_directories( add_subdirectory(lasp) add_subdirectory(test) -find_program(PYTHON "python") - if (PYTHON) # set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup.py.in") set(SETUP_PY "${CMAKE_CURRENT_BINARY_DIR}/setup.py") @@ -134,4 +126,3 @@ endif() ############################## End compiler settings - diff --git a/img/LASP.pdf b/img/LASP.pdf new file mode 100644 index 0000000..8fd406c Binary files /dev/null and b/img/LASP.pdf differ diff --git a/lasp/CMakeLists.txt b/lasp/CMakeLists.txt index 6f4dfb0..3f53de0 100644 --- a/lasp/CMakeLists.txt +++ b/lasp/CMakeLists.txt @@ -1,16 +1,15 @@ add_subdirectory(c) - configure_file(config.pxi.in config.pxi) - set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/cmake") -find_package( PythonLibs REQUIRED ) +find_package(PythonLibs REQUIRED ) include(UseCython) - +find_package(Numpy REQUIRED ) include_directories( - . - c - ) + ${PYTHON_NUMPY_INCLUDE_DIR} + . + c + ) # add the command to generate the source code # add_custom_command ( diff --git a/lasp/cmake/FindNumpy.cmake b/lasp/cmake/FindNumpy.cmake new file mode 100644 index 0000000..a95d6ac --- /dev/null +++ b/lasp/cmake/FindNumpy.cmake @@ -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) diff --git a/lasp/filter/bandpass_fir.py b/lasp/filter/bandpass_fir.py index d608015..308b644 100644 --- a/lasp/filter/bandpass_fir.py +++ b/lasp/filter/bandpass_fir.py @@ -175,6 +175,8 @@ class ThirdOctaveBankDesigner(FilterBankDesigner): '6.3k', '8k', '10k', '12.5k', '16k', '20k'] + assert len(self.xs) == len(self.nominal_txt) + @property def b(self): # Band division factor, 3 for one-third octave bands diff --git a/lasp/lasp_measurement.py b/lasp/lasp_measurement.py index 074ec6f..2e0fd2c 100644 --- a/lasp/lasp_measurement.py +++ b/lasp/lasp_measurement.py @@ -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 numpy as np 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 """ - def __init__(self, faudio): + def __init__(self, f): """ Initialize a BlockIter object @@ -55,8 +57,8 @@ class BlockIter: faudio: Audio dataset in the h5 file, accessed as f['audio'] """ self.i = 0 - self.nblocks = faudio.shape[0] - self.fa = faudio + self.nblocks = f['audio'].shape[0] + self.fa = f['audio'] def __iter__(self): return self @@ -121,35 +123,55 @@ class Measurement: # Full filepath self.fn = fn + # Base filename self.fn_base = os.path.split(fn)[1] # Open the h5 file in read-plus mode, to allow for changing the # measurement comment. - f = h5.File(fn, 'r+') - self.f = f + with h5.File(fn, 'r+') as f: + # Check for video data + try: + f['video'] + self.has_video = True + except KeyError: + self.has_video = False - # Check for video data - try: - f['video'] - self.has_video = True - except KeyError: - self.has_video = False + self.nblocks, self.blocksize, self.nchannels = f['audio'].shape + dtype = f['audio'].dtype + self.sampwidth = getSampWidth(dtype) - self.nblocks, self.blocksize, self.nchannels = f['audio'].shape - dtype = f['audio'].dtype - self.sampwidth = getSampWidth(dtype) + self.samplerate = f.attrs['samplerate'] + self.N = (self.nblocks*self.blocksize) + self.T = self.N/self.samplerate - self.samplerate = f.attrs['samplerate'] - self.N = (self.nblocks*self.blocksize) - self.T = self.N/self.samplerate + # comment = read-write thing + try: + self._comment = f.attrs['comment'] + except KeyError: + f.attrs['comment'] = '' + self._comment = '' - # comment = read-write thing - try: - self._comment = self.f.attrs['comment'] - except KeyError: - self.f.attrs['comment'] = '' - self._comment = '' + # Sensitivity + try: + self._sens = f.attrs['sensitivity'] + except KeyError: + self._sens = np.ones(self.nchannels) + + 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 def comment(self): @@ -157,19 +179,23 @@ class Measurement: @comment.setter def comment(self, cmt): - self.f.attrs['comment'] = cmt - self._comment = cmt + with self.file('r+') as f: + f.attrs['comment'] = cmt + self._comment = cmt @property 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 def time(self): """ Returns the measurement time in seconds since the epoch. """ - return self.f.attrs['time'] + return self._time @property def prms(self): @@ -186,7 +212,7 @@ class Measurement: except AttributeError: pass - sens = self.sensitivity + sens = self._sens pms = 0. for block in self.iterBlocks(): 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. """ if block is not None: - blocks = self.f['audio'][block] + with self.file() as f: + blocks = f['audio'][block] else: blocks = [] - for block in self.iterBlocks(): - blocks.append(block) + with self.file() as f: + for block in self.iterBlocks(f): + blocks.append(block) blocks = np.asarray(blocks) - blocks = blocks.reshape(self.nblocks*self.blocksize, self.nchannels) @@ -220,13 +247,19 @@ class Measurement: else: raise RuntimeError( 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, :] return blocks - def iterBlocks(self): - return BlockIter(self.f['audio']) + def iterBlocks(self, opened_file): + """ + Iterate over all the audio blocks in the opened file + + Args: + opened_file: The h5File with the data + """ + return BlockIter(opened_file) @property def sensitivity(self): @@ -235,10 +268,7 @@ class Measurement: between -1.0 and 1.0 to Pascal. If the sensitivity is not stored in the measurement file, this function returns 1.0 """ - try: - return self.f.attrs['sensitivity'] - except KeyError: - return np.ones(self.nchannels) + return self._sens @sensitivity.setter def sensitivity(self, sens): @@ -257,8 +287,8 @@ class Measurement: valid &= isinstance(sens.dtype, float) if not valid: raise ValueError('Invalid sensitivity value(s) given') - - self.f.attrs['sensitivity'] = sens + with self.file('r+') as f: + f.attrs['sensitivity'] = sens def exportAsWave(self, fn=None, force=False, sampwidth=None): """ @@ -286,8 +316,8 @@ class Measurement: if os.path.exists(fn) and not force: raise RuntimeError(f'File already exists: {fn}') - - audio = self.f['audio'] + with self.file() as f: + audio = self.f['audio'][:] if isinstance(audio.dtype, float): if sampwidth is None: @@ -373,8 +403,5 @@ class Measurement: ad[0] = dat return Measurement(mfn) - def __del__(self): - try: - self.f.close() - except AttributeError: - pass + # def __del__(self): + # self.f.close() diff --git a/lasp/lasp_octavefilter.py b/lasp/lasp_octavefilter.py index 629b073..4a970c4 100644 --- a/lasp/lasp_octavefilter.py +++ b/lasp/lasp_octavefilter.py @@ -17,21 +17,19 @@ class FilterBank: Single channel octave filter bank implementation """ - def __init__(self, fs, designer): + def __init__(self, fs): """ Initialize a OctaveFilterBank object. Args: fs: Sampling frequency of base signal - designer: FIR Filter designer for filterbank """ assert np.isclose(fs, 48000), "Only sampling frequency" \ " available is 48 kHz" - self.fs = fs - maxdecimation = designer.decimation(designer.xs[0]) + maxdecimation = self.decimation(self.xs[0]) self.decimators = [] for dec in maxdecimation: self.decimators.append(Decimator(1, dec)) @@ -44,8 +42,8 @@ class FilterBank: self.filterbanks = [] # Sort the x values in categories according to the required decimation - for x in designer.xs: - dec = designer.decimation(x) + for x in self.xs: + dec = self.decimation(x) if len(dec) == 1 and dec[0] == 1: xs_d1.append(x) 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] for xs in xs_all: 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): # These are the filters that do not require lasp_decimation # prior to filtering - nominals.append(designer.nominal(x)) - firs[:, i] = designer.createFilter(fs, x) + nominals.append(self.nominal(x)) + firs[:, i] = self.createFilter(fs, x) filterbank = {'fb': pyxFilterBank(firs, 1024), 'xs': xs, 'nominals': nominals} @@ -78,16 +76,10 @@ class FilterBank: # Filter output counters self.dec = [1, 4, 16, 64, 256] - Pdec = 128 - Pfil = 256 - # Initial filtered number - self.Nf = [-Pfil, -(Pdec/4 + Pfil), - -(Pfil + Pdec/4 + Pdec/16), - -(Pfil + Pdec/64 + Pdec/4 + Pdec/16), - -(Pfil + Pdec/256 + Pdec/64 + Pdec/4 + Pdec/16)] - - self.delay_set = [False, False, False, False, False] + # These intial delays are found experimentally using a toneburst + # response. + self.Nf = [915, 806, 780, 582, 338] def filterd(self, dec_stage, data): """ @@ -110,10 +102,12 @@ class FilterBank: oldNf = self.Nf[dec_stage] tstart = oldNf/fd tend = tstart + Nf/fd + t = np.linspace(tstart, tend, Nf, endpoint=False) self.Nf[dec_stage] += Nf 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 def filter_(self, data): @@ -142,13 +136,13 @@ class FilterBank: return output -class OctaveFilterBank(FilterBank): +class OctaveFilterBank(FilterBank, OctaveBankDesigner): def __init__(self, fs): - designer = OctaveBankDesigner() - super().__init__(fs, designer) + OctaveBankDesigner.__init__(self) + FilterBank.__init__(self, fs) -class ThirdOctaveFilterBank(FilterBank): +class ThirdOctaveFilterBank(FilterBank, ThirdOctaveBankDesigner): def __init__(self, fs): - designer = ThirdOctaveBankDesigner() - super().__init__(fs, designer) + ThirdOctaveBankDesigner.__init__(self) + FilterBank.__init__(self, fs) diff --git a/lasp/lasp_reverb.py b/lasp/lasp_reverb.py index 6c476b0..2de03df 100644 --- a/lasp/lasp_reverb.py +++ b/lasp/lasp_reverb.py @@ -24,7 +24,8 @@ class ReverbTime: 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 # Number of time samples self._N = self._level.shape[0] diff --git a/lasp/lasp_slm.py b/lasp/lasp_slm.py index df651b7..911ff3a 100644 --- a/lasp/lasp_slm.py +++ b/lasp/lasp_slm.py @@ -14,7 +14,8 @@ from .lasp_weighcal import WeighCal from .lasp_gui_tools import wait_cursor from .lasp_figure import PlotOptions, Plotable from .ui_slmwidget import Ui_SlmWidget - +from .filter.bandpass_fir import OctaveBankDesigner, ThirdOctaveBankDesigner +from .lasp_octavefilter import OctaveFilterBank, ThirdOctaveFilterBank __all__ = ['SLM', 'SlmWidget'] @@ -29,30 +30,39 @@ class Dummy: 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, - tw=TimeWeighting.default, - ): + def __init__(self, fs, tw=TimeWeighting.default): """ - Initialize a sound level meter object. Number of channels comes from - the weighcal object + Initialize a sound level meter object. + Number of channels comes from the weighcal object. Args: fs: Sampling frequency [Hz] - weighcal: WeighCal instance used for calibration and frequency - weighting. - nchannels: Number of channels to allocate filters for + tw: Time Weighting to apply """ - nchannels = weighcal.nchannels - self.nchannels = nchannels - self._weighcal = weighcal + 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: - self._lps = [Dummy() for i in range(nchannels)] - self._Lmax = zeros(nchannels) + self._lp = Dummy() + 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 def Lmax(self): @@ -68,26 +78,26 @@ class SLM: Args: data: """ - if data.ndim == 1: - data = data[:, np.newaxis] + assert data.ndim == 2 + assert data.shape[1] == 1 - data_weighted = self._weighcal.filter_(data) + P_REFsq = P_REF**2 # Squared - sq = data_weighted**2 - if sq.shape[0] == 0: - return np.array([]) + sq = data**2 - 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 - for chan, lp in enumerate(self._lps): - tw.append(lp.filter_(sq[:, chan])[:, 0]) + tw = self._lp.filter_(sq) - tw = np.asarray(tw).transpose() - - Level = 10*np.log10(tw/P_REF**2) + Level = 10*np.log10(tw/P_REFsq) + # Update maximum level curmax = np.max(Level) if curmax > self._Lmax: self._Lmax = curmax @@ -103,8 +113,8 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget): super().__init__(parent) self.setupUi(self) - FreqWeighting.fillComboBox(self.tfreqweighting) - FreqWeighting.fillComboBox(self.eqfreqweighting) + self.eqFreqBandChanged(0) + self.tFreqBandChanged(0) self.setMeas(None) def init(self, fm): @@ -115,6 +125,9 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget): fm.registerCombo(self.tfigure) fm.registerCombo(self.eqfigure) + self.tbandstart.setEnabled(False) + self.tbandstop.setEnabled(False) + def setMeas(self, meas): """ Set the current measurement for this widget. @@ -151,16 +164,33 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget): channel = self.eqchannel.currentIndex() fw = FreqWeighting.getCurrent(self.eqfreqweighting) - startpos = self.eqstarttime.value - stoppos = self.eqstoptime.value - 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") + istart, istop = self.getStartStopIndices(meas, self.eqstarttime, + self.eqstoptime) + 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(): # This one exctracts the calfile and sensitivity from global @@ -170,8 +200,28 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget): fs=fs, calfile=calfile, sens=sens) 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.fig.add(pta) fig.show() def computeT(self): @@ -184,46 +234,95 @@ class SlmWidget(ComputeWidget, Ui_SlmWidget): tw = TimeWeighting.getCurrent(self.ttimeweighting) 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 dsf = self.tdownsampling.value() - # gb = self.slmFre with wait_cursor(): # This one exctracts the calfile and sensitivity from global # variables defined at the top. # TODO: Change this to a more # robust variant. + + praw = meas.praw()[istart:istop, [channel]] + weighcal = WeighCal(fw, nchannels=1, fs=fs, calfile=calfile, sens=sens) - slm = SLM(fs, weighcal, tw) - praw = meas.praw()[:, [channel]] + weighted = weighcal.filter_(praw) - # Filter, downsample data - filtered = slm.addData(praw)[::dsf, :] - N = filtered.shape[0] - time = getTime(float(fs)/dsf, N) - Lmax = slm.Lmax + if bands == 0: + slm = SLM(fs, tw) + level = slm.addData(weighted)[::dsf] - pto = PlotOptions() - pto.ylabel = f'L{fw[0]} [dB({fw[0]})]' - pto.xlim = (time[0], time[-1]) + # Filter, downsample data + N = level.shape[0] + 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') - fig.fig.add(pta) + else: + 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() - stats = f"""Statistical results: -============================= -Applied frequency weighting: {fw[1]} -Applied time weighting: {tw[1]} -Applied Downsampling factor: {dsf} -Maximum level (L{fw[0]} max): {Lmax:4.4} [dB({fw[0]})] - - """ - self.results.setPlainText(stats) +# stats = f"""Statistical results: +# ============================= +# Applied frequency weighting: {fw[1]} +# Applied time weighting: {tw[1]} +# Applied Downsampling factor: {dsf} +# Maximum level (L{fw[0]} max): {Lmax:4.4} [dB({fw[0]})] +# +# """ +# self.results.setPlainText(stats) def compute(self): """ @@ -234,3 +333,61 @@ Maximum level (L{fw[0]} max): {Lmax:4.4} [dB({fw[0]})] self.computeT() elif self.eqtab.isVisible(): 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) diff --git a/lasp/lasp_weighcal.py b/lasp/lasp_weighcal.py index 2b14d62..9868520 100644 --- a/lasp/lasp_weighcal.py +++ b/lasp/lasp_weighcal.py @@ -4,10 +4,10 @@ Weighting and calibration filter in one @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 .filter_design.fir_design import (arbitrary_fir_design, - freqResponse as frp) +from .filter.fir_design import (arbitrary_fir_design, + freqResponse as frp) from lasp.lasp_config import ones, empty from .wrappers import FilterBank import numpy as np @@ -17,7 +17,7 @@ __all__ = ['WeighCal'] class WeighCal: """ - Time weighting and calibration FIR filter + Frequency weighting and calibration FIR filter """ def __init__(self, fw=FreqWeighting.default, nchannels=1, @@ -48,15 +48,17 @@ class WeighCal: # Objective function for the frequency response frp_obj = self.frpObj(freq_design) - self._firs = [] + P = 2048 # Filter length (number of taps) + + self._firs = np.empty((P, self.nchannels)) self._fbs = [] 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], 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 @@ -66,22 +68,21 @@ class WeighCal: Args: 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 dimension 1 in case of a single channel. + have the same number of columns as the number of channels. First + axis is assumed to be the time axis Retuns: Filtered data for each channel """ nchan = self.nchannels - if data.ndim == 1: - data = data[:, np.newaxis] + assert data.ndim == 2 assert data.shape[1] == nchan assert data.shape[0] > 0 filtered = [] 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 if filtered.ndim == 1: filtered = filtered[:, np.newaxis] @@ -133,7 +134,7 @@ class WeighCal: """ # Objective function for the frequency response frp_objective = self.frpCalObj(freq_design) * \ - self.frpWeightingObj(freq_design)[:, np.newaxis] + self.frpWeightingObj(freq_design)[:, np.newaxis] frp_objective[-1] = 0. return frp_objective diff --git a/lasp/plot/bar.py b/lasp/plot/bar.py index 6b3a791..b8b09ad 100644 --- a/lasp/plot/bar.py +++ b/lasp/plot/bar.py @@ -55,10 +55,15 @@ class BarScene(QGraphicsScene): xvals: labels and x positions of the bars G: Number of bars per x value 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) + self.setSceneRect(QRect(0,0,*size)) # self.setBackgroundBrush(ASCEEColors.bgBrush(0, size[0])) self.ylim = ylim @@ -68,8 +73,7 @@ class BarScene(QGraphicsScene): self.bgs = [] self.size = size - xsize = size[0] - ysize = size[1] + xsize, ysize = size self.xsize = xsize self.ysize = ysize diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx index 23638c7..134af1d 100644 --- a/lasp/wrappers.pyx +++ b/lasp/wrappers.pyx @@ -377,13 +377,14 @@ cdef class SPLowpass: if 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: return np.array([],dtype=NUMPY_FLOAT_TYPE) - + 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) # # Steal the pointer from output diff --git a/setup.py b/setup.py index 754c39e..d2bb70f 100644 --- a/setup.py +++ b/setup.py @@ -16,9 +16,8 @@ setup( author_email="j.a.dejong@ascee.nl", # Project uses reStructuredText, so ensure that the docutils get # 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', - description=descr, keywords="", url="http://www.ascee.nl/lasp/", # project home page, if any