From bdef0b45f36451e2ebbb41eddb010acdb0cfbdc1 Mon Sep 17 00:00:00 2001 From: Thijs Hekman Date: Wed, 22 Mar 2023 16:23:57 +0100 Subject: [PATCH 1/4] Added first-order HP and LP compensation filters to the biquad class --- src/lasp/filter/biquad.py | 57 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 54 insertions(+), 3 deletions(-) diff --git a/src/lasp/filter/biquad.py b/src/lasp/filter/biquad.py index 1ce16f5..209724d 100644 --- a/src/lasp/filter/biquad.py +++ b/src/lasp/filter/biquad.py @@ -23,7 +23,8 @@ y[n] = 1/ba[3] * ( ba[0] * x[n] + ba[1] * x[n-1] + ba[2] * x[n-2] + """ __all__ = ['peaking', 'biquadTF', 'notch', 'lowpass', 'highpass', - 'highshelf', 'lowshelf', 'LPcompensator', 'HPcompensator'] + 'highshelf', 'lowshelf', 'LP1compensator', 'LP2compensator', + 'HP1compensator', 'HP2compensator'] from numpy import array, cos, pi, sin, sqrt from scipy.interpolate import interp1d @@ -157,7 +158,32 @@ def lowshelf(fs, f0, Q, gain): a2 = (A+1) + (A-1)*cos(w0) - 2*sqrt(A)*alpha return array([b0/a0, b1/a0, b2/a0, a0/a0, a1/a0, a2/a0]) -def LPcompensator(fs, f0o, Qo, f0n, Qn): +def LP1compensator(fs, f0o, f0n): + """ + Shelving type filter that, when multiplied with a first-order low-pass + filter, alters the response of that filter to a different first-order + low-pass filter. + + Args: + fs: Sampling frequency [Hz] + f0o: Cut-off frequency of the original filter [Hz] + f0n: Desired cut-off frequency [Hz] + """ + + omg0o = 2*pi*f0o + omg0n = 2*pi*f0n + + z = -omg0o + p = -omg0n + k = p/z + + zd, pd, kd = bilinear_zpk(z, p, k, fs) + + sos = zpk2sos(zd,pd,kd) + + return sos[0] + +def LP2compensator(fs, f0o, Qo, f0n, Qn): """ Shelving type filter that, when multiplied with a second-order low-pass filter, alters the response of that filter to a different second-order @@ -194,7 +220,32 @@ def LPcompensator(fs, f0o, Qo, f0n, Qn): return sos[0] -def HPcompensator(fs, f0o, Qo, f0n, Qn): +def HP1compensator(fs, f0o, f0n): + """ + Shelving type filter that, when multiplied with a first-order high-pass + filter, alters the response of that filter to a different first-order + high-pass filter. + + Args: + fs: Sampling frequency [Hz] + f0o: Cut-on frequency of the original filter [Hz] + f0n: Desired cut-on frequency [Hz] + """ + + omg0o = 2*pi*f0o + omg0n = 2*pi*f0n + + z = -omg0o + p = -omg0n + k = 1 + + zd, pd, kd = bilinear_zpk(z, p, k, fs) + + sos = zpk2sos(zd,pd,kd) + + return sos[0] + +def HP2compensator(fs, f0o, Qo, f0n, Qn): """ Shelving type filter that, when multiplied with a second-order high-pass filter, alters the response of that filter to a different second-order From f1348ede80842717e235c85f2efba77ff74663db Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Mon, 3 Apr 2023 13:16:39 +0200 Subject: [PATCH 2/4] Cached time string of measurement time stamp --- src/lasp/lasp_measurement.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/lasp/lasp_measurement.py b/src/lasp/lasp_measurement.py index 8d2a012..924d466 100644 --- a/src/lasp/lasp_measurement.py +++ b/src/lasp/lasp_measurement.py @@ -52,6 +52,7 @@ import os, time, wave, logging from .lasp_common import SIQtys, Qty, getFreq from .lasp_cpp import Window, DaqChannel, LASP_VERSION_MAJOR, AvPowerSpectra from typing import List +from functools import lru_cache def getSampWidth(dtype): @@ -248,6 +249,9 @@ class Measurement: except KeyError: self._sens = np.ones(self.nchannels) + # The time is cached AND ALWAYS ASSUMED TO BE AN IMMUTABLE OBJECT. + # It is also cached. Changing the measurement timestamp should not + # be done. self._time = f.attrs['time'] # Quantity stored as channel. @@ -375,6 +379,7 @@ class Measurement: self._comment = cmt @property + @lru_cache() def recTime(self): """Returns the total recording time of the measurement, in float seconds.""" @@ -385,6 +390,19 @@ class Measurement: """Returns the measurement time in seconds since the epoch.""" return self._time + @property + @lru_cache() + def timestr(self): + """ + Return a properly formatted string of the measurement time, in order of + + year-month-day hour etc. + + """ + time_struct = time.localtime(self.time) + time_string = time.strftime('%Y-%m-%d %H:%M:%S', time_struct) + return time_string + def rms(self, channels=None, substract_average=False): """Returns the root mean square values for each channel From e09b00d801be1e343ece92a9a58d53d813ae6896 Mon Sep 17 00:00:00 2001 From: Casper Date: Fri, 14 Apr 2023 17:04:27 +0200 Subject: [PATCH 3/4] Changed argument of Measurement.exportAsWave() from newsampwidth to dtype, to allow export as float --- src/lasp/lasp_measurement.py | 76 ++++++++++++++++++++---------------- 1 file changed, 42 insertions(+), 34 deletions(-) diff --git a/src/lasp/lasp_measurement.py b/src/lasp/lasp_measurement.py index 924d466..cbe0d91 100644 --- a/src/lasp/lasp_measurement.py +++ b/src/lasp/lasp_measurement.py @@ -264,7 +264,7 @@ class Measurement: try: qtys_json = f.attrs['qtys'] # Load quantity data - self._qtys = [Qty.from_json(qty_json) for qty_json in qtys_json] + self._qtys = [Qty.from_json(qty_json) for qty_json in qtys_json] except KeyError: # If quantity data is not available, this is an 'old' # measurement file. @@ -274,7 +274,6 @@ class Measurement: self._qtys = [SIQtys.default() for i in range(self.nchannels)] logging.debug(f'Physical quantity data not available in measurement file. Assuming {SIQtys.default}') - def setAttribute(self, atrname, value): """ Set an attribute in the measurement file, and keep a local copy in @@ -395,7 +394,7 @@ class Measurement: def timestr(self): """ Return a properly formatted string of the measurement time, in order of - + year-month-day hour etc. """ @@ -426,7 +425,7 @@ class Measurement: meansquare += np.sum(block**2, axis=0) / self.N avg = sum_/N - # In fact, this is not the complete RMS, as in includes the DC + # In fact, this is not the complete RMS, as in includes the DC # If p = p_dc + p_osc, then rms(p_osc) = sqrt(ms(p)-ms(p_dc)) if substract_average: meansquare -= avg**2 @@ -437,7 +436,7 @@ class Measurement: return self.rms(substract_average=True) def rawData(self, channels=None, **kwargs): - """Returns the raw data as stored in the measurement file, + """Returns the raw data as stored in the measurement file, without any transformations applied args: @@ -490,10 +489,10 @@ class Measurement: window: Window type overlap: Overlap percentage (value between 0.0 and up to and including 100.0) - weighting: + weighting: Returns: - Cross-power-spectra. C[freq, ch_i, ch_j] = C_ij + Cross-power-spectra. C[freq, ch_i, ch_j] = C_ij """ nfft = kwargs.pop('nfft', 2048) @@ -535,14 +534,14 @@ class Measurement: # Estimate noise power in block blocks = [signal[i*N:(i+1)*N] for i in range(Nblocks)] - + if noiseCorrection: - # The difference between the measured signal in the previous block and + # The difference between the measured signal in the previous block and # the current block en = [None] + [blocks[i] - blocks[i-1] for i in range(1,Nblocks)] - + noise_est = [None] + [-np.average(en[i]*en[i+1]) for i in range(1,len(en)-1)] - + # Create weighting coefficients sum_inverse_noise = sum([1/n for n in noise_est[1:]]) c_n = [1/(ni*sum_inverse_noise) for ni in noise_est[1:]] @@ -551,7 +550,7 @@ class Measurement: assert np.isclose(sum(c_n), 1.0) assert Nblocks-2 == len(c_n) - + # Average signal over blocks avg = np.zeros((blocks[0].shape), dtype=float) for n in range(0, Nblocks-2): @@ -576,7 +575,7 @@ class Measurement: CS = ps.compute(avg) freq = getFreq(self.samplerate, N) - return freq, CS + return freq, CS @property @@ -632,8 +631,8 @@ class Measurement: return False - def exportAsWave(self, fn=None, force=False, newsampwidth=None, - normalize=True, **kwargs): + def exportAsWave(self, fn=None, force=False, dtype=None, + normalize=False, **kwargs): """Export measurement file as wave. In case the measurement data is stored as floats, the values are scaled to the proper integer (PCM) data format. @@ -645,9 +644,8 @@ class Measurement: force: If True, overwrites any existing files with the given name , otherwise a RuntimeError is raised. - newsampwidth: 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 error is thrown + dtype: if not None, convert data to this data type. + Options are 'int16', 'int32', 'float32'. normalize: If set: normalize the level to something sensible. """ @@ -661,39 +659,49 @@ class Measurement: if os.path.exists(fn) and not force: raise RuntimeError(f'File already exists: {fn}') - if not np.isclose(self.samplerate%1,0): raise RuntimeError(f'Sample rates should be approximately integer for exporting to Wave to work') # TODO: With VERY large measurment files, this is not possible! Is this # a theoretical case? + # TODO: add sensitivity? Then use self.data() instead of self.rawData() data = self.rawData(**kwargs) - if np.issubdtype(data.dtype, np.floating) and newsampwidth is None: - raise ValueError('Newsampwidth parameter should be given for floating point raw data') if normalize: # Scale back to maximum of absolute value maxabs = np.max(np.abs(data)) data = data / maxabs # "data /= maxabs" fails if dtpyes differ - if newsampwidth is not None: - # Convert to floats, then to new sample width - sensone = np.ones_like(self.sensitivity) - data = scaleBlockSens(data, sensone) + if dtype==None: + dtype = data.dtype # keep existing + logging.debug(f"dtype not passed as arg; using dtype = {dtype}") - if newsampwidth == 2: - newtype = np.int16 - elif newsampwidth == 4: - newtype = np.int32 - else: - raise ValueError('Invalid sample width, should be 2 or 4') + # dtype conversion + if dtype=='int16': + newtype = np.int16 + newsampwidth = 2 + elif dtype=='int32': + newtype = np.int32 + newsampwidth = 4 + elif dtype=='float32': + newtype = np.float32 + elif dtype=='float64': + newtype = np.float64 + else: + logging.debug(f"cannot handle this dtype {dtype}") + pass + # Convert range to [-1, 1] + # TODO: this is wrong for float data where full scale > 1 + sensone = np.ones_like(self.sensitivity) + data = scaleBlockSens(data, sensone) + + if dtype=='int16' or dtype=='int32': + # Scale data to integer range and convert to integers scalefac = 2**(8*newsampwidth-1)-1 - - # Scale data to integer range, and then convert to integers data = (data*scalefac).astype(newtype) - wavfile.write(fn, int(self.samplerate), data) + wavfile.write(fn, int(self.samplerate), data.astype(newtype)) @staticmethod def fromtxt(fn, From fb9920d00a14243cf3bbe30e0a1e55a1e7a120a9 Mon Sep 17 00:00:00 2001 From: "J.A. de Jong - Redu-Sone B.V., ASCEE V.O.F" Date: Tue, 18 Apr 2023 11:09:01 +0200 Subject: [PATCH 4/4] Bugfix: Made Uldaq sample rate check depending on actually provided sample rates --- src/lasp/device/lasp_uldaq.cpp | 5 +---- src/lasp/device/lasp_uldaq_impl.cpp | 2 +- src/lasp/device/lasp_uldaq_impl.h | 8 ++++++++ 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/lasp/device/lasp_uldaq.cpp b/src/lasp/device/lasp_uldaq.cpp index 41de7dc..34e1f8e 100644 --- a/src/lasp/device/lasp_uldaq.cpp +++ b/src/lasp/device/lasp_uldaq.cpp @@ -67,10 +67,7 @@ void fillUlDaqDeviceInfo(DeviceInfoList &devinfolist) { DataTypeDescriptor::DataType::dtype_fl64); devinfo.prefDataTypeIndex = 0; - devinfo.availableSampleRates = {8000, 10000, 11025, 16000, 20000, - 22050, 24000, 32000, 44056, 44100, - 47250, 48000, 50000, 50400, 51000}; - + devinfo.availableSampleRates = ULDAQ_SAMPLERATES; devinfo.prefSampleRateIndex = 11; devinfo.availableFramesPerBlock = {512, 1024, 2048, 4096, 8192}; diff --git a/src/lasp/device/lasp_uldaq_impl.cpp b/src/lasp/device/lasp_uldaq_impl.cpp index 9041ff9..6177b89 100644 --- a/src/lasp/device/lasp_uldaq_impl.cpp +++ b/src/lasp/device/lasp_uldaq_impl.cpp @@ -70,7 +70,7 @@ DT9837A::DT9837A(const DeviceInfo &devinfo, const DaqConfiguration &config) throw rte("Unsensible number of samples per block chosen"); } - if (samplerate() < 10000 || samplerate() > 51000) { + if (samplerate() < ULDAQ_SAMPLERATES.at(0) || samplerate() > ULDAQ_SAMPLERATES.at(ULDAQ_SAMPLERATES.size()-1)) { throw rte("Invalid sample rate"); } diff --git a/src/lasp/device/lasp_uldaq_impl.h b/src/lasp/device/lasp_uldaq_impl.h index 9fce2e5..90c308a 100644 --- a/src/lasp/device/lasp_uldaq_impl.h +++ b/src/lasp/device/lasp_uldaq_impl.h @@ -14,6 +14,14 @@ using std::cerr; using std::endl; using rte = std::runtime_error; +/** + * @brief List of available sampling frequencies for DT9837A + */ +const std::vector ULDAQ_SAMPLERATES = {8000, 10000, 11025, 16000, 20000, + 22050, 24000, 32000, 44056, 44100, + 47250, 48000, 50000, 50400, 51000}; + + /** * @brief UlDaq-specific device information. Adds a copy of the underlying * DaqDeDaqDeviceDescriptor.