From e09b00d801be1e343ece92a9a58d53d813ae6896 Mon Sep 17 00:00:00 2001 From: Casper Date: Fri, 14 Apr 2023 17:04:27 +0200 Subject: [PATCH] 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,