ASCEE/lasp
ASCEE/lasp
5
1
Fork 0
Code Issues 5 Pull Requests Actions Releases 14 Activity

Merge remote-tracking branch 'origin/develop' into develop
All checks were successful
continuous-integration/drone/push Build is passing
Details

Browse Source
This commit is contained in:
Anne de Jong 2023-04-19 15:56:35 +02:00
commit 3844827505
1 changed files with 42 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
src/lasp/lasp_measurement.py
View File

@ -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,