diff --git a/lasp/lasp_common.py b/lasp/lasp_common.py index d8f50b5..d6fdf3e 100644 --- a/lasp/lasp_common.py +++ b/lasp/lasp_common.py @@ -114,8 +114,8 @@ def getTime(fs, N, start=0): N: Number of time samples start: Optional start ofset in number of samples """ - - return np.linspace(start/fs, N/fs, N, endpoint=False) + assert N > 0 and fs > 0 + return np.linspace(start, start + N/fs, N, endpoint=False) def getFreq(fs, nfft): diff --git a/lasp/lasp_measurement.py b/lasp/lasp_measurement.py index ec72d4a..a11411f 100644 --- a/lasp/lasp_measurement.py +++ b/lasp/lasp_measurement.py @@ -41,13 +41,13 @@ import numpy as np from .lasp_config import LASP_NUMPY_FLOAT_TYPE import wave import os +import time class BlockIter: """ Iterate over the blocks in the audio data of a h5 file """ - def __init__(self, f): """ Initialize a BlockIter object @@ -69,7 +69,7 @@ class BlockIter: if self.i == self.nblocks: raise StopIteration self.i += 1 - return self.fa[self.i-1][:, :] + return self.fa[self.i - 1][:, :] def getSampWidth(dtype): @@ -107,10 +107,10 @@ def scaleBlockSens(block, sens): assert sens.size == block.shape[1] if np.issubdtype(block.dtype.type, np.integer): sw = getSampWidth(block.dtype) - fac = 2**(8*sw - 1) - 1 + fac = 2**(8 * sw - 1) - 1 else: fac = 1. - return block.astype(LASP_NUMPY_FLOAT_TYPE)/fac/sens[np.newaxis, :] + return block.astype(LASP_NUMPY_FLOAT_TYPE) / fac / sens[np.newaxis, :] def exportAsWave(fn, fs, data, force=False): @@ -133,7 +133,6 @@ class Measurement: Provides access to measurement data stored in the h5 measurement file format. """ - def __init__(self, fn): """ Initialize a Measurement object based on the filename @@ -162,8 +161,8 @@ class Measurement: self.sampwidth = getSampWidth(dtype) self.samplerate = f.attrs['samplerate'] - self.N = (self.nblocks*self.blocksize) - self.T = self.N/self.samplerate + self.N = (self.nblocks * self.blocksize) + self.T = self.N / self.samplerate # comment = read-write thing try: @@ -232,7 +231,7 @@ class Measurement: Returns the total recording time of the measurement, in float seconds. """ - return self.blocksize*self.nblocks/self.samplerate + return self.blocksize * self.nblocks / self.samplerate @property def time(self): @@ -273,7 +272,7 @@ class Measurement: with self.file() as f: for block in self.iterBlocks(f): block = self.scaleBlock(block) - pms += np.sum(block**2, axis=0)/self.N + pms += np.sum(block**2, axis=0) / self.N self._prms = np.sqrt(pms) return self._prms @@ -291,7 +290,7 @@ class Measurement: for block in self.iterBlocks(f): blocks.append(block) blocks = np.asarray(blocks) - blocks = blocks.reshape(self.nblocks*self.blocksize, + blocks = blocks.reshape(self.nblocks * self.blocksize, self.nchannels) # Apply scaling (sensitivity, integer -> float) @@ -326,7 +325,7 @@ class Measurement: equal to the number of channels. """ if isinstance(sens, float): - sens = sens*np.ones(self.nchannels) + sens = sens * np.ones(self.nchannels) valid = sens.ndim == 1 valid &= sens.shape[0] == self.nchannels @@ -386,24 +385,30 @@ class Measurement: # Scale with maximum value only if we have a nonzero maximum value. if max == 0.: max = 1. - scalefac = 2**(8*sampwidth)/max + scalefac = 2**(8 * sampwidth) / max with wave.open(fn, 'w') as wf: - wf.setparams((self.nchannels, self.sampwidth, - self.samplerate, 0, 'NONE', 'NONE')) + 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) + block = (block * scalefac).astype(itype) wf.writeframes(np.asfortranarray(block).tobytes()) @staticmethod - def fromtxt(fn, skiprows, samplerate, sensitivity, mfn=None, + def fromtxt(fn, + skiprows, + samplerate, + sensitivity, + mfn=None, timestamp=None, - delimiter='\t', firstcoltime=True): + delimiter='\t', + firstcoltime=True): """ - Converts a txt file to a LASP Measurement object and returns the - measurement. + Converts a txt file to a LASP Measurement file, opens the associated + Measurement object and returns it. The measurement file will have + the same file name as the txt file, except with h5 extension. Args: fn: Filename of text file @@ -430,11 +435,16 @@ class Measurement: dat = np.loadtxt(fn, skiprows=skiprows, delimiter=delimiter) if firstcoltime: time = dat[:, 0] - if not np.isclose(time[1] - time[0], 1/samplerate): + if not np.isclose(time[1] - time[0], 1 / samplerate): raise ValueError('Samplerate given does not agree with ' 'samplerate in file') + + # Chop off first column dat = dat[:, 1:] nchannels = dat.shape[1] + if nchannels != sensitivity.shape[0]: + raise ValueError( + f'Invalid sensitivity length given. Should be: {nchannels}') with h5.File(mfn, 'w') as hf: hf.attrs['samplerate'] = samplerate @@ -442,13 +452,62 @@ class Measurement: hf.attrs['time'] = timestamp hf.attrs['blocksize'] = 1 hf.attrs['nchannels'] = nchannels - ad = hf.create_dataset('audio', - (1, dat.shape[0], dat.shape[1]), + 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): - # self.f.close() + @staticmethod + def fromnpy(data, + samplerate, + sensitivity, + mfn, + timestamp=None): + """ + Converts a numpy array to a LASP Measurement file, opens the associated + Measurement object and returns it. The measurement file will have + the same file name as the txt file, except with h5 extension. + + + Args: + data: Numpy array, first column is sample, second is channel. Can + also be specified with a single column for single-channel data + samplerate: Sampling frequency in [Hz] + sensitivity: 1D array of channel sensitivities [Pa^-1] + mfn: Filepath where measurement file is stored. + 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 os.path.exists(mfn): + raise ValueError(f'File {mfn} already exist.') + if timestamp is None: + timestamp = int(time.time()) + + if data.ndim != 2: + data = data[:, np.newaxis] + + nchannels = data.shape[1] + if nchannels != sensitivity.shape[0]: + raise ValueError( + f'Invalid sensitivity length given. Should be: {nchannels}') + + 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, data.shape[0], data.shape[1]), + dtype=data.dtype, + maxshape=(1, data.shape[0], data.shape[1]), + compression='gzip') + ad[0] = data + return Measurement(mfn) diff --git a/lasp/lasp_reverb.py b/lasp/lasp_reverb.py index fecd991..dbc829d 100644 --- a/lasp/lasp_reverb.py +++ b/lasp/lasp_reverb.py @@ -14,7 +14,7 @@ class ReverbTime: Tool to estimate the reverberation time """ - def __init__(self, fs, level): + def __init__(self, fs, level, channel=0): """ Initialize Reverberation time computer. @@ -22,22 +22,25 @@ class ReverbTime: fs: Sampling frequency [Hz] level: (Optionally weighted) level values as a function of time, in dB. + channel: Channel index to compute from """ 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 + + self._level = level[:, channel][:, np.newaxis] # Number of time samples + self._channel = channel self._N = self._level.shape[0] self._t = getTime(fs, self._N, 0) + print(f't: {self._t}') - def compute(self, tstart, tstop): + def compute(self, istart, istop): """ Compute the reverberation time using a least-squares solver Args: - tstart: Start time of reverberation interval - stop: Stop time of reverberation interval + istart: Start time index reverberation interval + istop: Stop time index of reverberation interval Returns: dictionary with result values, contains: @@ -48,36 +51,33 @@ class ReverbTime: - 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))]) + A = np.hstack([x, np.ones(x.shape)]) + + print(A.shape) + print(points.shape) # derivative is dB/s of increase/decrease sol, residuals, rank, s = np.linalg.lstsq(A, points) + + print(f'sol: {sol}') # Derivative of the decay in dB/s - derivative = sol[0] + derivative = sol[0][0] # Start level in dB - const = sol[1] + const = sol[1][0] # Time to reach a decay of 60 dB (reverb. time) T60 = -60./derivative - - return {'istart': istart, + res = {'istart': istart, 'istop': istop, 'const': const, 'derivative': derivative, 'T60': T60} + print(res) + return res diff --git a/lasp/wrappers.pyx b/lasp/wrappers.pyx index 7fd4805..73243ca 100644 --- a/lasp/wrappers.pyx +++ b/lasp/wrappers.pyx @@ -10,6 +10,7 @@ cdef extern from "cblas.h": int openblas_get_num_threads() void openblas_set_num_threads(int) + # If we touch this variable: we get segfaults when running from # Spyder! # openblas_set_num_threads(8) diff --git a/scripts/lasp_record b/scripts/lasp_record index 5f5db82..d070a4c 100755 --- a/scripts/lasp_record +++ b/scripts/lasp_record @@ -1,8 +1,5 @@ #!/usr/bin/python import argparse -from lasp.lasp_record import Recording -from lasp.lasp_avstream import AvStream -from lasp.device.lasp_daqconfig import default_soundcard, roga_plugndaq, umik parser = argparse.ArgumentParser( description='Acquire data and store a measurement file' ) @@ -16,11 +13,25 @@ parser.add_argument('--comment', '-c', type=str, device_help = 'DAQ Device to record from' parser.add_argument('--input-daq', '-i', help=device_help, type=str, - choices=['roga', 'umik', 'default'], default='roga') + choices=['roga', 'umik', 'nidaq', 'default'], default='roga') args = parser.parse_args() device_str = args.input_daq + +if device_str == 'nidaq': + + # Not-integrated way to record with the NI USB 4431 DAQ device + from lasp.device.record_ni import USBDAQRecording + rec = USBDAQRecording(args.filename, [0, 1]) + rec.start() + exit(0) + + +from lasp.lasp_record import Recording +from lasp.lasp_avstream import AvStream +from lasp.device.lasp_daqconfig import default_soundcard, roga_plugndaq, umik + if 'roga' == device_str: device = roga_plugndaq elif 'default' == device_str: