From 3357a0ded6c79f322f7744615f3b0f154481a6bb Mon Sep 17 00:00:00 2001 From: "J.A. de Jong" Date: Wed, 2 May 2018 16:29:53 +0200 Subject: [PATCH] Split up GUI in different Widgets, added revtime, added figure list possibilities, added Qt Resources, added About panel, added lots of comments, export and import of measurements --- .gitattributes | 1 + .gitignore | 2 + lasp/CMakeLists.txt | 28 ++-- lasp/c/CMakeLists.txt | 2 +- lasp/lasp_common.py | 69 +++++++--- lasp/lasp_config.py | 13 +- lasp/lasp_measurement.py | 269 ++++++++++++++++++++++++++++++++++----- lasp/lasp_record.py | 69 +++++----- lasp/lasp_reverb.py | 81 ++++++++++++ lasp/lasp_slm.py | 141 ++++++++++++++++++-- lasp/lasp_weighcal.py | 56 ++++---- 11 files changed, 606 insertions(+), 125 deletions(-) create mode 100644 .gitattributes create mode 100644 lasp/lasp_reverb.py diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..947f120 --- /dev/null +++ b/.gitattributes @@ -0,0 +1 @@ +*.png binary diff --git a/.gitignore b/.gitignore index 3021a61..ed3eb64 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ __pycache__ cython_debug lasp/config.pxi lasp/wrappers.c +lasp/resources_rc.py *.so test/test_bf test/test_fft @@ -19,3 +20,4 @@ doc LASP.egg-info lasp_octave_fir.* lasp/ui_* +resources_rc.py diff --git a/lasp/CMakeLists.txt b/lasp/CMakeLists.txt index 65536cd..6f4dfb0 100644 --- a/lasp/CMakeLists.txt +++ b/lasp/CMakeLists.txt @@ -18,17 +18,27 @@ include_directories( # COMMAND MakeTable ${CMAKE_CURRENT_BINARY_DIR}/Table.h # DEPENDS MakeTable # ) - + +set(ui_files ui_apsrt ui_mainwindow ui_figure ui_about ui_apswidget ui_revtime ui_slmwidget) +foreach(fn ${ui_files}) + add_custom_command( + OUTPUT "${fn}.py" + COMMAND pyside-uic ${CMAKE_CURRENT_SOURCE_DIR}/ui/${fn}.ui -o ${CMAKE_CURRENT_SOURCE_DIR}/${fn}.py + DEPENDS "ui/${fn}.ui" + ) + add_custom_target(${fn} ALL DEPENDS "${fn}.py") +endforeach(fn) + add_custom_command( - OUTPUT "aps_ui.py" - COMMAND pyside-uic ${CMAKE_CURRENT_SOURCE_DIR}/ui/aps_ui.ui -o ${CMAKE_CURRENT_SOURCE_DIR}/aps_ui.py - DEPENDS "ui/aps_ui.ui" - ) -add_custom_target(ui ALL DEPENDS "aps_ui.py") + OUTPUT "${PROJECT_SOURCE_DIR}/resources_rc.py" + COMMAND pyside-rcc -py3 ui/resources.qrc -o ${PROJECT_SOURCE_DIR}/resources_rc.py + DEPENDS "ui/resources.qrc" +) +add_custom_target(resources_rc ALL DEPENDS "${PROJECT_SOURCE_DIR}/resources_rc.py") -set_source_files_properties(wrappers.c PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} ${CYTHON_EXTRA_C_FLAGS}") + + +set_source_files_properties(wrappers.c PROPERTIES COMPILE_FLAGS "${CMAKE_C_FLAGS} ${CYTHON_EXTRA_C_FLAGS}") cython_add_module(wrappers wrappers.pyx) target_link_libraries(wrappers lasp_lib ) - - diff --git a/lasp/c/CMakeLists.txt b/lasp/c/CMakeLists.txt index ec00dd5..bf34cea 100644 --- a/lasp/c/CMakeLists.txt +++ b/lasp/c/CMakeLists.txt @@ -17,7 +17,7 @@ add_library(lasp_lib lasp_worker.c lasp_dfifo.c lasp_filterbank.c - lasp_octave_fir.c + # lasp_octave_fir.c lasp_decimation.c lasp_sp_lowpass.c ) diff --git a/lasp/lasp_common.py b/lasp/lasp_common.py index 0730b04..8b23986 100644 --- a/lasp/lasp_common.py +++ b/lasp/lasp_common.py @@ -1,21 +1,56 @@ # -*- coding: utf-8 -*- #!/usr/bin/env python3 import numpy as np +from .wrappers import Window as wWindow """ Common definitions used throughout the code. """ -__all__ = ['P_REF','FreqWeighting','TimeWeighting','getTime'] +__all__ = ['P_REF', 'FreqWeighting', 'TimeWeighting', 'getTime', 'calfile', + 'sens'] # Reference sound pressure level P_REF = 2e-5 +# Todo: fix This +calfile = '/home/anne/wip/UMIK-1/cal/7027430_90deg.txt' +sens = 0.053690387255872614 + + +class Window: + 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): + """ + Fill Windows to a combobox + + Args: + cb: QComboBox to fill + """ + cb.clear() + for tw in Window.types: + cb.addItem(tw[1], tw) + cb.setCurrentIndex(Window.default) + + def getCurrent(cb): + return Window.types[cb.currentIndex()] + + class TimeWeighting: - none = (None,'Raw (no time weighting)') - ufast = (30e-3,'30 ms') - fast = (0.125,'Fast') - slow = (1.0,'Slow') - types = (none, ufast, fast, slow) + none = (None, 'Raw (no time weighting)') + uufast = (1e-4, '0.1 ms') + ufast = (30e-3, '30 ms') + fast = (0.125, 'Fast') + slow = (1.0, 'Slow') + types = (none, uufast, ufast, fast, slow) default = 2 @staticmethod @@ -28,20 +63,21 @@ class TimeWeighting: """ cb.clear() for tw in TimeWeighting.types: - cb.addItem(tw[1],tw) + cb.addItem(tw[1], tw) cb.setCurrentIndex(TimeWeighting.default) def getCurrent(cb): return TimeWeighting.types[cb.currentIndex()] + class FreqWeighting: """ Frequency weighting types """ - A=('A','A-weighting') - C=('C','C-weighting') - Z=('Z','Z-weighting') - types = (A,C,Z) + A = ('A', 'A-weighting') + C = ('C', 'C-weighting') + Z = ('Z', 'Z-weighting') + types = (A, C, Z) default = 0 @staticmethod @@ -54,13 +90,14 @@ class FreqWeighting: """ cb.clear() for fw in FreqWeighting.types: - cb.addItem(fw[1],fw) + cb.addItem(fw[1], fw) cb.setCurrentIndex(FreqWeighting.default) def getCurrent(cb): return FreqWeighting.types[cb.currentIndex()] -def getTime(fs,N,start=0): + +def getTime(fs, N, start=0): """ Return a time array for given number of points and sampling frequency. @@ -70,7 +107,7 @@ def getTime(fs,N,start=0): start: Optional start ofset in number of samples """ - return np.linspace(start/fs,N/fs,N,endpoint=False) + return np.linspace(start/fs, N/fs, N, endpoint=False) def getFreq(fs, nfft): @@ -81,6 +118,6 @@ def getFreq(fs, nfft): fs: Sampling frequency [Hz] nfft: Fft length (int) """ - df = fs/nfft # frequency resolution - K = nfft//2+1 # number of frequency bins + df = fs/nfft # frequency resolution + K = nfft//2+1 # number of frequency bins return np.linspace(0, (K-1)*df, K) diff --git a/lasp/lasp_config.py b/lasp/lasp_config.py index 1f0edeb..49000b5 100644 --- a/lasp/lasp_config.py +++ b/lasp/lasp_config.py @@ -10,8 +10,15 @@ import numpy as np LASP_NUMPY_FLOAT_TYPE = np.float64 def zeros(shape): - return np.zeros(shape,dtype=LASP_NUMPY_FLOAT_TYPE) + return np.zeros(shape, dtype=LASP_NUMPY_FLOAT_TYPE) def ones(shape): - return np.ones(shape,dtype=LASP_NUMPY_FLOAT_TYPE) + return np.ones(shape, dtype=LASP_NUMPY_FLOAT_TYPE) def empty(shape): - return np.empty(shape,dtype=LASP_NUMPY_FLOAT_TYPE) \ No newline at end of file + return np.empty(shape, dtype=LASP_NUMPY_FLOAT_TYPE) + +class config: + default_figsize = (12,7) + default_comment_lenght_measoverview = 50 + default_xlim_freq_plot = (50, 10000) + default_ylim_freq_plot = (20, 100) + versiontxt = '0.1' diff --git a/lasp/lasp_measurement.py b/lasp/lasp_measurement.py index 6ddd40f..074ec6f 100644 --- a/lasp/lasp_measurement.py +++ b/lasp/lasp_measurement.py @@ -4,36 +4,94 @@ Author: J.A. de Jong - ASCEE Description: Measurement class + +The ASCEE hdf5 measurement file format contains the following fields: + +- Attributes: + +'samplerate': The audio data sample rate in Hz. +'nchannels': The number of audio channels in the file +'sensitivity': (Optionally) the stored sensitivity of the record channels. + This can be a single value, or an array of sensitivities for + each channel. Both representations are allowed. + +- Datasets: + +'audio': 3-dimensional array of blocks of audio data. The first axis is the +block index, the second axis the sample number and the third axis is the channel +number. The data type is either int16, int32 or float64 / float32. In case the +data is stored as integers. The raw data should be scaled with the maximum value +that can be stored for the integer bit depth to get a number between -1.0 and +1.0. + +'video': 4-dimensional array of video frames. The first index is the frame + number, the second the x-value of the pixel and the third is the + y-value of the pixel. Then, the last axis is the color. This axis has + length 3 and the colors are stored as (r,g,b). Where typically a + color depth of 256 is used (np.uint8 data format) + +The video dataset can possibly be not present in the data. + + + """ import h5py as h5 import numpy as np from .lasp_config import LASP_NUMPY_FLOAT_TYPE -import wave,os +import wave +import os + class BlockIter: - def __init__(self,nblocks,faudio): + """ + Iterate over the blocks in the audio data of a h5 file + """ + + def __init__(self, faudio): + """ + Initialize a BlockIter object + + Args: + faudio: Audio dataset in the h5 file, accessed as f['audio'] + """ self.i = 0 - self.nblocks = nblocks + self.nblocks = faudio.shape[0] self.fa = faudio + def __iter__(self): return self + def __next__(self): + """ + Return the next block + """ if self.i == self.nblocks: raise StopIteration - self.i+=1 - return self.fa[self.i-1][:,:] + self.i += 1 + return self.fa[self.i-1][:, :] def getSampWidth(dtype): + """ + Returns the width of a single sample in bytes. + + Args: + dtype: numpy dtype + + Returns: + Size of a sample in bytes (int) + """ if dtype == np.int32: return 4 elif dtype == np.int16: return 2 + elif dtype == np.float64: + return 8 else: - raise ValueError('Invalid data type: %s' %dtype) + raise ValueError('Invalid data type: %s' % dtype) -def exportAsWave(fn,fs,data,force=False): +def exportAsWave(fn, fs, data, force=False): if not '.wav' in fn[-4:]: fn += '.wav' @@ -43,21 +101,35 @@ def exportAsWave(fn,fs,data,force=False): if os.path.exists(fn) and not force: raise RuntimeError('File already exists: %s', fn) - with wave.open(fn,'w') as wf: - wf.setparams((nchannels,sampwidth,fs,0,'NONE','NONE')) + with wave.open(fn, 'w') as wf: + wf.setparams((nchannels, sampwidth, fs, 0, 'NONE', 'NONE')) wf.writeframes(np.asfortranarray(data).tobytes()) class Measurement: - def __init__(self, fn): + """ + Provides access to measurement data stored in the h5 measurement file + format. + """ - if not '.h5' in fn: + def __init__(self, fn): + """ + Initialize a Measurement object based on the filename + """ + if '.h5' not in fn: fn += '.h5' + # Full filepath self.fn = fn + # Base filename + self.fn_base = os.path.split(fn)[1] - f = h5.File(fn,'r+') + # Open the h5 file in read-plus mode, to allow for changing the + # measurement comment. + f = h5.File(fn, 'r+') self.f = f + + # Check for video data try: f['video'] self.has_video = True @@ -74,17 +146,19 @@ class Measurement: # comment = read-write thing try: - self.f.attrs['comment'] + self._comment = self.f.attrs['comment'] except KeyError: self.f.attrs['comment'] = '' + self._comment = '' @property def comment(self): - return self.f.attrs['comment'] + return self._comment @comment.setter def comment(self, cmt): self.f.attrs['comment'] = cmt + self._comment = cmt @property def recTime(self): @@ -92,23 +166,37 @@ class Measurement: @property def time(self): + """ + Returns the measurement time in seconds since the epoch. + """ return self.f.attrs['time'] @property def prms(self): + """ + Returns the root mean square of the uncalibrated rms sound pressure + level (equivalend SPL). + + Returns: + 1D array with rms values for each channel + """ + # try: return self._prms except AttributeError: pass + sens = self.sensitivity pms = 0. for block in self.iterBlocks(): - pms += np.sum(block/sens)**2/self.N + pms += np.sum(block/sens[np.newaxis, :], axis=0)**2/self.N self._prms = np.sqrt(pms) return self._prms - - def praw(self,block=None): + def praw(self, block=None): + """ + Returns the raw uncalibrated data, converted to floating point format. + """ if block is not None: blocks = self.f['audio'][block] else: @@ -118,51 +206,172 @@ class Measurement: blocks = np.asarray(blocks) blocks = blocks.reshape(self.nblocks*self.blocksize, - self.nchannels) + self.nchannels) + + # When the data is stored as integers, we assume dB full-scale scaling. + # Hence, when we convert the data to floats, we divide by the maximum + # possible value. if blocks.dtype == np.int32: fac = 2**31 elif blocks.dtype == np.int16: fac = 2**15 + elif blocks.dtype == np.float64: + fac = 1.0 else: - raise RuntimeError('Unimplemented data type from recording: %s' %str(blocks.dtype)) - - blocks = blocks.astype(LASP_NUMPY_FLOAT_TYPE)/fac/self.sensitivity + raise RuntimeError( + f'Unimplemented data type from recording: {blocks.dtype}.') + sens = self.sensitivity + blocks = blocks.astype(LASP_NUMPY_FLOAT_TYPE)/fac/sens[np.newaxis, :] return blocks def iterBlocks(self): - return BlockIter(self.nblocks,self.f['audio']) + return BlockIter(self.f['audio']) @property def sensitivity(self): + """ + Sensitivity of the data in Pa^-1, from floating point data scaled + 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: - return 1.0 + except KeyError: + return np.ones(self.nchannels) @sensitivity.setter def sensitivity(self, sens): + """ + Set the sensitivity of the measurement in the file + + Args: + sens: sensitivity data, should be a float, or an array of floats + equal to the number of channels. + """ + if isinstance(sens, float): + sens = sens*np.ones(self.nchannels) + + valid = sens.ndim == 1 + valid &= sens.shape[0] == self.nchannels + valid &= isinstance(sens.dtype, float) + if not valid: + raise ValueError('Invalid sensitivity value(s) given') + self.f.attrs['sensitivity'] = sens - def exportAsWave(self, fn=None, force=False): + def exportAsWave(self, fn=None, force=False, sampwidth=None): """ - Export measurement file as wave + Export measurement file as wave. In case the measurement data is stored + as floats, the values are scaled between 0 and 1 + + Args: + fn: If given, this will be the filename to write to. If the + filename does not end with '.wav', this extension is added. + + force: If True, overwrites any existing files with the given name + , otherwise a RuntimeError is raised. + + sampwidth: sample width in bytes with which to export the data. + This should only be given in case the measurement data is stored as + floating point values, otherwise an + """ if fn is None: fn = self.fn fn = os.path.splitext(fn)[0] - if not '.wav' in fn[-4:]: + if '.wav' not in fn[-4:]: fn += '.wav' if os.path.exists(fn) and not force: - raise RuntimeError('File already exists: %s', fn) + raise RuntimeError(f'File already exists: {fn}') - with wave.open(fn,'w') as wf: - wf.setparams((self.nchannels,self.sampwidth,self.samplerate,0,'NONE','NONE')) + audio = self.f['audio'] + + if isinstance(audio.dtype, float): + if sampwidth is None: + raise ValueError('sampwidth parameter should be given ' + 'for float data in measurement file') + elif sampwidth == 2: + itype = np.int16 + elif sampwidth == 4: + itype = np.int32 + else: + raise ValueError('Invalid sample width, should be 2 or 4') + + # Find maximum + max = 0. for block in self.iterBlocks(): + blockmax = np.max(np.abs(block)) + if blockmax > max: + max = blockmax + # Scale with maximum value only if we have a nonzero maximum value. + if max == 0.: + max = 1. + scalefac = 2**(8*sampwidth)/max + + with wave.open(fn, 'w') as wf: + wf.setparams((self.nchannels, self.sampwidth, + self.samplerate, 0, 'NONE', 'NONE')) + for block in self.iterBlocks(): + if isinstance(block.dtype, float): + # Convert block to integral data type + block = (block*scalefac).astype(itype) wf.writeframes(np.asfortranarray(block).tobytes()) + @staticmethod + def fromtxt(fn, skiprows, samplerate, sensitivity, mfn=None, + timestamp=None, + delimiter='\t', firstcoltime=True): + """ + Converts a txt file to a LASP Measurement object and returns the + measurement. + + Args: + fn: Filename of text file + skiprows: Number of header rows in text file to skip + samplerate: Sampling frequency in [Hz] + sensitivity: 1D array of channel sensitivities + mfn: Filepath where measurement file is stored. If not given, + a h5 file will be created along fn, which shares its basename + timestamp: If given, a custom timestamp for the measurement + (integer containing seconds since epoch). If not given, the + timestamp is obtained from the last modification time. + delimiter: Column delimiter + firstcoltime: If true, the first column is the treated as the + sample time. + + """ + if not os.path.exists(fn): + raise ValueError(f'File {fn} does not exist.') + if timestamp is None: + timestamp = os.path.getmtime(fn) + if mfn is None: + mfn = os.path.splitext(fn)[0] + '.h5' + + dat = np.loadtxt(fn, skiprows=skiprows, delimiter=delimiter) + if firstcoltime: + time = dat[:, 0] + if not np.isclose(time[1] - time[0], 1/samplerate): + raise ValueError('Samplerate given does not agree with ' + 'samplerate in file') + dat = dat[:, 1:] + nchannels = dat.shape[1] + + with h5.File(mfn, 'w') as hf: + hf.attrs['samplerate'] = samplerate + hf.attrs['sensitivity'] = sensitivity + hf.attrs['time'] = timestamp + hf.attrs['blocksize'] = 1 + hf.attrs['nchannels'] = nchannels + ad = hf.create_dataset('audio', + (1, dat.shape[0], dat.shape[1]), + dtype=dat.dtype, + maxshape=(1, dat.shape[0], dat.shape[1]), + compression='gzip') + ad[0] = dat + return Measurement(mfn) def __del__(self): try: diff --git a/lasp/lasp_record.py b/lasp/lasp_record.py index cdbc9ad..db9b580 100644 --- a/lasp/lasp_record.py +++ b/lasp/lasp_record.py @@ -3,20 +3,23 @@ """ Created on Sat Mar 10 08:28:03 2018 -@author: Read data from image stream and record sound at the same time +Read data from stream and record sound and video at the same time """ import numpy as np from .lasp_atomic import Atomic from threading import Condition from .lasp_avstream import AvType -import h5py, time +import h5py +import time + -# %% class Recording: - def __init__(self, fn, stream, rectime = None): - + def __init__(self, fn, stream, rectime=None): + """ + + """ ext = '.h5' - if not ext in fn: + if ext not in fn: fn += ext self._stream = stream self.blocksize = stream.blocksize @@ -25,42 +28,43 @@ class Recording: self._running_cond = Condition() self.rectime = rectime self._fn = fn - + self._video_frame_positions = [] - + self._aframeno = Atomic(0) self._vframeno = 0 def start(self): stream = self._stream - with h5py.File(self._fn,'w') as f: + with h5py.File(self._fn, 'w') as f: self._ad = f.create_dataset('audio', - (1,stream.blocksize,stream.nchannels), + (1, stream.blocksize, stream.nchannels), dtype=stream.dtype, - maxshape=(None,stream.blocksize, + maxshape=(None, stream.blocksize, stream.nchannels), compression='gzip' ) if stream.hasVideo(): - video_x,video_y = stream.video_x,stream.video_y + video_x, video_y = stream.video_x, stream.video_y self._vd = f.create_dataset('video', - (1,video_y,video_x,3), + (1, video_y, video_x, 3), dtype='uint8', - maxshape=(None,video_y,video_x,3), + maxshape=( + None, video_y, video_x, 3), compression='gzip' ) - + f.attrs['samplerate'] = stream.samplerate f.attrs['nchannels'] = stream.nchannels f.attrs['blocksize'] = stream.blocksize f.attrs['time'] = time.time() self._running <<= True # Videothread is going to start - + if not stream.isStarted(): stream.start() - + stream.addCallback(self._callback) with self._running_cond: try: @@ -70,34 +74,33 @@ class Recording: except KeyboardInterrupt: print("Keyboard interrupt on record") self._running <<= False - + stream.removeCallback(self._callback) if stream.hasVideo(): f['video_frame_positions'] = self._video_frame_positions - print('\nEnding record') def stop(self): self._running <<= False with self._running_cond: self._running_cond.notify() - + def _callback(self, _type, data, aframe, vframe): if not self._stream.isStarted(): self._running <<= False with self._running_cond: self._running_cond.notify() - + if _type == AvType.audio: self._aCallback(data, aframe) elif _type == AvType.video: self._vCallback(data) - + def _aCallback(self, frames, aframe): -# print(self._aframeno()) - print('.',end='') + # print(self._aframeno()) + print('.', end='') curT = self._aframeno()*self.blocksize/self.samplerate if self.rectime is not None and curT > self.rectime: # We are done! @@ -105,20 +108,20 @@ class Recording: with self._running_cond: self._running_cond.notify() return - - self._ad.resize(self._aframeno()+1,axis=0) - self._ad[self._aframeno(),:,:] = frames + + self._ad.resize(self._aframeno()+1, axis=0) + self._ad[self._aframeno(), :, :] = frames self._aframeno += 1 - - def _vCallback(self,frame): + def _vCallback(self, frame): self._video_frame_positions.append(self._aframeno()) vframeno = self._vframeno - self._vd.resize(vframeno+1,axis=0) - self._vd[vframeno,:,:] = frame + self._vd.resize(vframeno+1, axis=0) + self._vd[vframeno, :, :] = frame self._vframeno += 1 - + + if __name__ == '__main__': - rec = Recording('test',5) + rec = Recording('test', 5) rec.start() diff --git a/lasp/lasp_reverb.py b/lasp/lasp_reverb.py new file mode 100644 index 0000000..ada61d0 --- /dev/null +++ b/lasp/lasp_reverb.py @@ -0,0 +1,81 @@ +# -*- coding: utf-8 -*- +"""! +Author: J.A. de Jong - ASCEE + +Description: +Reverberation time estimation tool using least squares +""" +from .lasp_common import getTime +import numpy as np + + +class ReverbTime: + """ + Tool to estimate the reverberation time + """ + def __init__(self, fs, level): + """ + Initialize Reverberation time computer. + + Args: + fs: Sampling frequency [Hz] + level: (Optionally weighted) level values as a function of time, in + dB. + + """ + assert level.ndim == 1, 'Invalid number of dimensions in level' + self._level = level + # Number of time samples + self._N = self._level.shape[0] + self._t = getTime(fs, self._N, 0) + + def compute(self, tstart, tstop): + """ + Compute the reverberation time using a least-squares solver + + Args: + tstart: Start time of reverberation interval + stop: Stop time of reverberation interval + + Returns: + dictionary with result values, contains: + - istart: start index of reberberation interval + - istop: stop index of reverb. interval + - T60: Reverberation time + - const: Constant value + - derivative: rate of change of the level in dB/s. + """ + + # Find start and stop indices. Coerce if they are outside + # of the valid domain + istart = np.where(self._t >= tstart)[0][0] + istop = np.where(self._t >= tstop)[0] + + if istop.shape[0] == 0: + istop = self._level.shape[0] + else: + istop = istop[0] + + points = self._level[istart:istop] + x = self._t[istart:istop][:, np.newaxis] + + # Solve the least-squares problem, by creating a matrix of + A = np.hstack([x, np.ones((x.shape))]) + + # derivative is dB/s of increase/decrease + sol, residuals, rank, s = np.linalg.lstsq(A, points) + + # Derivative of the decay in dB/s + derivative = sol[0] + + # Start level in dB + const = sol[1] + + # Time to reach a decay of 60 dB (reverb. time) + T60 = -60./derivative + + return {'istart': istart, + 'istop': istop, + 'const': const, + 'derivative': derivative, + 'T60': T60} diff --git a/lasp/lasp_slm.py b/lasp/lasp_slm.py index 7e84e11..7a4ed4b 100644 --- a/lasp/lasp_slm.py +++ b/lasp/lasp_slm.py @@ -4,38 +4,73 @@ Sound level meter implementation @author: J.A. de Jong - ASCEE """ -from .wrappers import FilterBank, SPLowpass +from .wrappers import SPLowpass +from .lasp_computewidget import ComputeWidget import numpy as np from .lasp_config import zeros -from .lasp_common import TimeWeighting, P_REF -__all__ = ['SLM'] +from .lasp_common import (FreqWeighting, sens, calfile, + TimeWeighting, getTime, P_REF) +from .lasp_weighcal import WeighCal +from .lasp_gui_tools import wait_cursor +from .lasp_figure import FigureDialog, PlotOptions +from .ui_slmwidget import Ui_SlmWidget + +__all__ = ['SLM', 'SlmWidget'] + class Dummy: + """ + Emulate filtering, but does not filter anything at all. + """ + def filter_(self, data): - return data + return data[:, np.newaxis] + class SLM: """ Sound Level Meter, implements the single pole lowpass filter """ + def __init__(self, fs, weighcal, tw=TimeWeighting.default, - nchannels = 1, - freq=None, cal=None): + ): + """ + 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 + """ + nchannels = weighcal.nchannels + self.nchannels = nchannels self._weighcal = weighcal - if tw is not TimeWeighting.none: + if tw[0] is not TimeWeighting.none[0]: self._lps = [SPLowpass(fs, tw[0]) for i in range(nchannels)] else: - self._lpw = [Dummy() for i in range(nchannels)] + self._lps = [Dummy() for i in range(nchannels)] self._Lmax = zeros(nchannels) @property def Lmax(self): + """ + Returns the currently maximum recorded level + """ return self._Lmax def addData(self, data): - assert data.ndim == 2 + """ + Add new fresh timedata to the Sound Level Meter + + Args: + data: + """ + if data.ndim == 1: + data = data[:, np.newaxis] + data_weighted = self._weighcal.filter_(data) # Squared @@ -44,9 +79,11 @@ class SLM: return np.array([]) tw = [] + # Time-weight the signal - for chan,lp in enumerate(self._lps): - tw.append(lp.filter_(sq[:,chan])[:,0]) + for chan, lp in enumerate(self._lps): + tw.append(lp.filter_(sq[:, chan])[:, 0]) + tw = np.asarray(tw).transpose() Level = 10*np.log10(tw/P_REF**2) @@ -56,3 +93,85 @@ class SLM: self._Lmax = curmax return Level + + +class SlmWidget(ComputeWidget, Ui_SlmWidget): + def __init__(self): + """ + Initialize the SlmWidget. + """ + super().__init__() + self.setupUi(self) + FreqWeighting.fillComboBox(self.freqweighting) + + self.setMeas(None) + + def setMeas(self, meas): + """ + Set the current measurement for this widget. + + Args: + meas: if None, the Widget is disabled + """ + self.meas = meas + if meas is None: + self.setEnabled(False) + else: + self.setEnabled(True) + rt = meas.recTime + self.starttime.setRange(0, rt, 0) + self.stoptime.setRange(0, rt, rt) + + self.channel.clear() + for i in range(meas.nchannels): + self.channel.addItem(str(i)) + self.channel.setCurrentIndex(0) + + def compute(self): + """ + Compute Sound Level using settings. This method is + called whenever the Compute button is pushed in the SLM tab + """ + meas = self.meas + fs = meas.samplerate + channel = self.channel.currentIndex() + tw = TimeWeighting.getCurrent(self.timeweighting) + fw = FreqWeighting.getCurrent(self.freqweighting) + + # Downsampling factor of result + dsf = self.downsampling.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. + weighcal = WeighCal(fw, nchannels=1, + fs=fs, calfile=calfile, + sens=sens) + + slm = SLM(fs, weighcal, tw) + praw = meas.praw()[:, [channel]] + + # Filter, downsample data + filtered = slm.addData(praw)[::dsf, :] + N = filtered.shape[0] + time = getTime(float(fs)/dsf, N) + Lmax = slm.Lmax + + pto = PlotOptions() + pto.ylabel = f'L{fw[0]} [dB({fw[0]})]' + pto.xlim = (time[0], time[-1]) + fig, new = self.getFigure(pto) + fig.fig.plot(time, filtered) + 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) diff --git a/lasp/lasp_weighcal.py b/lasp/lasp_weighcal.py index b174dda..2b14d62 100644 --- a/lasp/lasp_weighcal.py +++ b/lasp/lasp_weighcal.py @@ -19,10 +19,10 @@ class WeighCal: """ Time weighting and calibration FIR filter """ - def __init__(self, fw = FreqWeighting.default, - nchannels = 1, - fs = 48000., - calfile = None, + def __init__(self, fw=FreqWeighting.default, + nchannels=1, + fs=48000., + calfile=None, sens=1.0): """ Initialize the frequency weighting and calibration FIR filters. @@ -42,7 +42,7 @@ class WeighCal: self.calfile = calfile # Frequencies used for the filter design - freq_design = np.linspace(0,17e3,3000) + freq_design = np.linspace(0, 17e3, 3000) freq_design[-1] = fs/2 # Objective function for the frequency response @@ -51,52 +51,63 @@ class WeighCal: self._firs = [] self._fbs = [] for chan in range(self.nchannels): - fir = arbitrary_fir_design(fs,2048,freq_design, - frp_obj[:,chan], - window='rectangular') + fir = arbitrary_fir_design(fs, 2048, freq_design, + frp_obj[:, chan], + window='rectangular') self._firs.append(fir) - self._fbs.append(FilterBank(fir[:,np.newaxis],4096)) + self._fbs.append(FilterBank(fir[:, np.newaxis], 4096)) self._freq_design = freq_design - def filter_(self,data): + def filter_(self, data): """ Filter data using the calibration and frequency weighting filter. + + 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. + + Retuns: + Filtered data for each channel + """ nchan = self.nchannels + if data.ndim == 1: + data = data[:, np.newaxis] 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] + filtered = filtered[:, np.newaxis] return filtered - def frpCalObj(self, freq_design): """ Computes the objective frequency response of the calibration filter """ calfile = self.calfile if calfile is not None: - cal = np.loadtxt(calfile,skiprows=2) - freq = cal[:,0] - cal = cal[:,1:] + cal = np.loadtxt(calfile, skiprows=2) + freq = cal[:, 0] + cal = cal[:, 1:] if cal.shape[1] != self.nchannels: raise ValueError('Number of channels in calibration file does' ' not equal to given number of channels') calfac = 10**(-cal/20) - filter_calfac = empty((freq_design.shape[0],self.nchannels)) + filter_calfac = empty((freq_design.shape[0], self.nchannels)) for chan in range(self.nchannels): - filter_calfac[:,chan] = np.interp(freq_design,freq,calfac[:,chan]) + filter_calfac[:, chan] = np.interp(freq_design, freq, + calfac[:, chan]) else: - filter_calfac = ones((freq_design.shape[0],self.nchannels,)) + filter_calfac = ones((freq_design.shape[0], self.nchannels,)) return filter_calfac @@ -122,7 +133,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 @@ -132,5 +143,6 @@ class WeighCal: Returns the frequency response of the designed FIR filter """ if freq is None: - freq = np.logspace(1,np.log10(self.fs/2),500) - return freq,frp(self.fs,freq,self._firs[chan]),self.frpObj(freq)[:,chan] + freq = np.logspace(1, np.log10(self.fs/2), 500) + return (freq, frp(self.fs, freq, self._firs[chan]), + self.frpObj(freq)[:, chan])