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

Lots of code and comments improvements

Browse Source
This commit is contained in:
Anne de Jong 2019-10-27 14:19:26 +01:00
parent 5af6aaf48e
commit 1f1922272d
5 changed files with 121 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
lasp/lasp_common.py
View File

@ -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):

107
lasp/lasp_measurement.py
View File

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

40
lasp/lasp_reverb.py
View File

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

1
lasp/wrappers.pyx
View File

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

19
scripts/lasp_record
View File

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