Added tools to bin narrow-band data into frequency bands

.gitignore

@@ -25,3 +25,4 @@ lasp/ui_*
 resources_rc.py
 lasp/device/lasp_daqdevice.c
 .ropeproject
+.spyproject

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

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)

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

lasp/tools/bin_narrow.py

@@ -0,0 +1,86 @@
+#!/usr/bin/python
+"""
+Bin narrow band power in octave/third octave band data
+"""
+from lasp.filter.bandpass_fir import (OctaveBankDesigner, 
+                                      ThirdOctaveBankDesigner)
+import numpy as np
+import warnings
+
+
+def binPower(freq, narrow_power, band=3, start_band=-16, stop_band=13):
+    """
+    Apply binning to narrow band frequency domain power results
+    
+    Args:
+        freq: Array of frequency indices
+        narrow_power: narrow-band power values in units of [W] or [Pa^2]
+        band: 1, or 3
+    
+    Returns:
+        ( ['25', '31.5', '40', '50', ... ],
+          [float(power_25), float(power_31p5), ...]) putting NAN values where
+           inaccurate.
+    """
+
+    if band == 3:
+        designer = ThirdOctaveBankDesigner()
+    elif band == 1:
+        designer = OctaveBankDesigner()
+    else:
+        raise ValueError("Parameter 'Band' should be either '1', or '3'")
+
+    freq = np.copy(freq)
+    narrow_power = np.copy(narrow_power)
+
+    
+    # Exact midband, lower and upper frequency limit of each band
+    fm = [designer.fm(x) for x in range(start_band, stop_band+1)]
+    fl = [designer.fl(x) for x in range(start_band, stop_band+1)]
+    fu = [designer.fu(x) for x in range(start_band, stop_band+1)]
+    fex = [designer.nominal_txt(x) for x in range(start_band, stop_band+1)]
+#    print(fl)
+    binned_power = np.zeros(len(fm), dtype=float)
+    ## Start: linear interpolation between bins while Parseval is conserved
+
+    # current frequency resolution
+    df_old = freq[1]-freq[0]
+    
+    # preferred new frequency resolution
+    df_new = .1
+
+    # ratio of resolutions
+
+    ratio = int(df_old/df_new)
+    # calculate new frequency bins
+    freq_new = np.linspace(freq[0],freq[-1],(len(freq)-1)*ratio+1)         
+    
+    # calculate the new bin data
+    interp_power = np.interp(freq_new, freq, narrow_power)/ratio
+    
+    # adapt first and last bin values so that Parseval still holds
+    interp_power[0] = binned_power[0]*(1.+1./ratio)/2.
+    interp_power[-1] = binned_power[-1]*(1.+1./ratio)/2.    
+
+    # check if Parseval still holds        
+    # print(np.sum(y, axis=0))
+    # print(np.sum(y_new, axis=0))
+    
+    ## Stop: linear interpolation between bins while Parseval is conserved        
+    binned_power = np.zeros(len(fm), dtype=float)
+    for k in range(len(fm)):
+#        print(k)
+        # find the bins which are in the corresponding band        
+        bins = (fl[k] <= freq_new) & (freq_new < fu[k])
+#        print(bins)
+        # sum the output values of these bins to obtain the band value
+        binned_power[k] = np.sum(interp_power[bins], axis=0)
+        # if no frequency bin falls in a certain band, skip previous bands 
+#        if not any(bins):
+#            binned_power[0:k+1] = np.nan
+
+    # check if data is valid
+    if(np.isnan(binned_power).all()):
+        warnings.warn('Invalid frequency array, we cannot bin these values')
+
+    return fm, fex, binned_power

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)

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: