[FIX] window the impulse response before filtering for C80, D50 and support (ISO 3382-1, A.3.4)

pytta/_plot.py

@@ -1136,7 +1136,7 @@ def waterfall(sigObjs, step=2 ** 9, n=2 ** 13, fmin=None, fmax=None, pmin=None,
         # Apply rolling window
         sliced = np.asarray(list(more_itertools.windowed(ht, n=n, step=step, fillvalue=0)))
         window = np.asarray(
-            [ss.tukey(n, alpha=winAlpha) for s in sliced])  # Tukey window (alpha=0 is a rectangular window)
+            [ss.windows.tukey(n, alpha=winAlpha) for s in sliced])  # Tukey window (alpha=0 is a rectangular window)
         windowed = sliced * window  # Apply window
         windowedFFT = abs(2 / n * fft(windowed))  # Apply FFT
         windowedFFT = windowedFFT[:, idx_min_FFT - fpad: idx_max_FFT + 1 + fpad]  # Crop the FFT

pytta/apps/roomir.py

@@ -67,7 +67,7 @@
 from os import getcwd, listdir, mkdir
 from os.path import isfile, join, exists
 import copy as cp
-import traceback
+from pytta.classes._instanceinfo import caller_name
 
 # Dict with the measurementKinds
 # TO DO FUTURE: add 'inchcalibration', 'outchcalibration'
@@ -2236,21 +2236,7 @@ def RT(self, roomirsGetDict, decay=20, IREndManualCut=None):
                 each source-receiver configuration, which is a key;
 
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(frame, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         # TR in all positions and averagßes
         # TR_avgs = {'S1R1': [Analysis_avg1, An_avg2, ..., An_avgn]}
@@ -2320,21 +2306,7 @@ def G_Lps(self, recalibirsGetDict):
                 level;
 
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(frame, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         recalibirs = recalibirsGetDict
         # Calculating Lps for each recalibir MeasuredThing
@@ -2401,21 +2373,7 @@ def G_Lpe_inSitu(self, roomirsGetDict, IREndManualCut=None):
                 key;
 
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(frame, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         #  Lpe in all positions and averagßes
         # Lpe_avgs = {'S1R1': [Analysis_avg1, An_avg2, ..., An_avgn]}
@@ -2461,21 +2419,7 @@ def G_Lpe_revCh(self, roomirsGetDict, IREndManualCut=None):
                 all the reverberation chamber's impulsive responses;
 
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(frame, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         SigObjs = []
         if isinstance(roomirsGetDict, dict):
@@ -2548,21 +2492,7 @@ def G(self, Lpe_avgs, Lpe_revCh, V_revCh, T_revCh, Lps_revCh, Lps_inSitu):
                 source-receiver configuration, which is a key;
 
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(frame, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         # G (omfg) calculation
         # G averages in all positions
@@ -2627,21 +2557,7 @@ def G_T_revCh(self, roomirsGetDict, IREndManualCut=None, T=20):
                 all the reverberation chamber's impulsive responses;
 
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(frame, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         roomirs = roomirsGetDict
         RTs = []

pytta/classes/_instanceinfo.py

@@ -13,11 +13,25 @@
 """
 
 
+import sys
 import time
 import traceback
 import threading
 
 
+def caller_name():
+    """
+    Name assigned on the module-level line that led to the current call:
+    'x' for ``x = pytta.generate.sweep()``. The text before '=' (the whole
+    line if there is none), or '' when the source line is not available.
+    """
+    frame = sys._getframe(2)  # the caller of the function asking
+    while frame.f_code.co_varnames and frame.f_back is not None:
+        frame = frame.f_back
+    line = traceback.extract_stack(frame, 1)[0].line or ''
+    return line.split('=')[0].strip()
+
+
 class InstanceCreationError(Exception):
     pass
 

pytta/classes/analysis.py

@@ -993,15 +993,29 @@ def estimate_energy_parameters(ir: SignalObj, bands: np.ndarray,
             A dict with parameters by name.
 
         """
-        listEDC, fhSignal = cumulative_integration(ir, bypassLundeby, plotLundeby, suppressWarnings, **kwargs)
+        compensation = []
+        listEDC, fhSignal = cumulative_integration(ir, bypassLundeby, plotLundeby, suppressWarnings,
+                                                   compensation=compensation, **kwargs)
         params = {}
         params['rms'] = fhSignal.rms()
         params['SPL'] = fhSignal.spl()
         params['Ts'] = central_time(fhSignal.timeSignal**2, fhSignal.timeVector)
-        params['D50'] = definition(listEDC, ir.samplingRate)
-        params['C80'] = clarity(listEDC, ir.samplingRate)
-        params['STearly'] = st_early(listEDC, ir.samplingRate)
-        params['STlate'] = st_late(listEDC, ir.samplingRate)
+        # The parameters with time windows window the broadband response and
+        # filter each window afterwards (ISO 3382-1:2009, A.3.4), so that the
+        # delay of the band filters does not move energy across the limits;
+        # the late windows keep the Lundeby truncation and tail correction.
+        early50 = windowed_band_energy(ir, 0, 0.05, **kwargs)
+        late50 = windowed_band_energy(ir, 0.05, None, compensation, **kwargs)
+        early80 = windowed_band_energy(ir, 0, 0.08, **kwargs)
+        late80 = windowed_band_energy(ir, 0.08, None, compensation, **kwargs)
+        first10 = windowed_band_energy(ir, 0, 0.01, **kwargs)
+        with np.errstate(divide='ignore', invalid='ignore'):
+            params['D50'] = np.round(100 * early50 / (early50 + late50), 2)
+            params['C80'] = np.round(10 * np.log10(early80 / late80), 2)
+            params['STearly'] = np.round(10 * np.log10(
+                windowed_band_energy(ir, 0.02, 0.1, **kwargs) / first10), 4)
+            params['STlate'] = np.round(10 * np.log10(
+                windowed_band_energy(ir, 0.1, 1.0, **kwargs) / first10), 4)
         params['EDT'] = reverberation_time('EDT', listEDC)
         params['T20'] = reverberation_time(20, listEDC)
         params['T30'] = reverberation_time(30, listEDC)
@@ -1417,6 +1431,7 @@ def energy_decay_calculation(band, timeSignal, timeVector, samplingRate,
     else:
         interIdx = 0
         lateRT = 1
+        lundebyParams = None
 
     if interIdx == 0:
         interIdx = -1
@@ -1437,15 +1452,19 @@ def energy_decay_calculation(band, timeSignal, timeVector, samplingRate,
             print(band, "[Hz] band: could not estimate C factor")
         C = 0
         energyDecay = np.zeros(truncatedTimeVector.size)
-    return (energyDecay, truncatedTimeVector, lundebyParams)
+    return (energyDecay, truncatedTimeVector, lundebyParams, (interIdx, C))
 
 
 def cumulative_integration(inputSignal,
                            bypassLundeby,
                            plotLundebyResults,
                            suppressWarnings=True,
+                           compensation=None,
                            **kwargs):
-    """Cumulative integration with proper corrections."""
+    """Cumulative integration with proper corrections.
+
+    If a list is given as compensation, the Lundeby truncation index and the
+    tail energy correction C of each band are appended to it."""
 
     def plot_lundeby():
         c0, c1, interIdx, BGL = lundebyParams
@@ -1476,7 +1495,7 @@ def plot_lundeby():
     for ch in range(hSignal.numChannels):
         signal = hSignal[ch]
         band = bands[ch]
-        energyDecay, energyVector, lundebyParams = \
+        energyDecay, energyVector, lundebyParams, bandCompensation = \
             energy_decay_calculation(band,
                                      signal.timeSignal,
                                      signal.timeVector,
@@ -1487,6 +1506,8 @@ def plot_lundeby():
                                      bypassLundeby,
                                      suppressWarnings=suppressWarnings)
         listEDC.append((energyDecay, energyVector))
+        if compensation is not None:
+            compensation.append(bandCompensation)
         if plotLundebyResults:  # Placed here because Numba can't handle plots.
             # plot_lundeby(band, timeVector, timeSignal,  samplingRate,
             #             lundebyParams)
@@ -1530,6 +1551,46 @@ def reverberation_time(decay, listEDC):
     return np.array(RT, dtype='float32')
 
 
+def windowed_band_energy(ir, start, stop=None, compensation=None, **kwargs):
+    """
+    Energy in each band of the part of a broadband impulse response between
+    start and stop seconds (stop None for the end), windowed before the band
+    filtering, as ISO 3382-1:2009, A.3.4, recommends for the early and late
+    energy of C80, D50 and the support parameters. The window is filtered
+    with the filters of _filter (kwargs nthOct, minFreq, maxFreq) after one
+    second of zeros is appended, so that the energy the filters delay past
+    the end of the response is kept.
+
+    compensation is the list of (Lundeby truncation index, tail correction
+    C) of each band that cumulative_integration fills; with it, the filtered
+    window is cut at the truncation index of each band and C is added, as
+    the energy decay curve does for the late part of the response. A third
+    element, the noise energy per sample, is subtracted from every sample
+    of the window up to the truncation, for a decay curve computed by noise
+    subtraction.
+
+    Returns an array with one energy per band.
+    """
+    fs = ir.samplingRate
+    x = np.asarray(ir.timeSignal, dtype=float)[:, 0]
+    i0 = int(round(start * fs))
+    i1 = x.size if stop is None else min(int(round(stop * fs)), x.size)
+    window = np.zeros(x.size + fs)
+    window[i0:i1] = x[i0:i1]
+    y = np.asarray(_filter(SignalObj(window, 'time', fs, signalType='energy'),
+                           **kwargs).timeSignal, dtype=float)
+    y = y.reshape(y.shape[0], -1)
+    if compensation is None:
+        return np.sum(y**2, axis=0)
+    energy = []
+    for band, comp in enumerate(compensation):
+        cut, C = comp[:2]
+        noise = comp[2] if len(comp) > 2 else 0.
+        end = x.size - 1 if cut < 0 else cut
+        energy.append(np.sum(y[:cut, band]**2) + C - noise * max(end - i0, 0))
+    return np.array(energy)
+
+
 def definition(listEDC: list, fs: int, t: int = 50) -> np.ndarray:
     """
     Room parameter.

pytta/classes/measurement.py

@@ -10,7 +10,7 @@
 from pytta.classes.streaming import Monitor, Streaming
 from pytta.classes.signal import SignalObj, ImpulsiveResponse
 from pytta import _h5utils as _h5
-import traceback
+from pytta.classes._instanceinfo import caller_name
 
 
 # Measurement class
@@ -418,21 +418,7 @@ def run(self):
         Run method: starts recording during Tmax seconds
         Outputs a signalObj with the recording content
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(frame, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         # Record
         recording = sd.rec(frames=self.numSamples,
@@ -556,21 +542,7 @@ def run(self):
         Starts reproducing the excitation signal and recording at the same time
         Outputs a signalObj with the recording content
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(framenline, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         timeStamp = time.ctime(time.time())
         recording = sd.playrec(self.excitation.timeSignal*
@@ -807,21 +779,7 @@ def run(self):
         Starts reproducing the excitation signal and recording at the same time
         Outputs the transferfunction ImpulsiveResponse
         """
-        # Code snippet to guarantee that generated object name is
-        # the declared at global scope
-        # for frame, line in traceback.walk_stack(None):
-        for framenline in traceback.walk_stack(None):
-            # varnames = frame.f_code.co_varnames
-            varnames = framenline[0].f_code.co_varnames
-            if varnames == ():
-                break
-        # creation_file, creation_line, creation_function, \
-        #     creation_text = \
-        extracted_text = \
-            traceback.extract_stack(framenline[0], 1)[0]
-            # traceback.extract_stack(frame, 1)[0]
-        # creation_name = creation_text.split("=")[0].strip()
-        creation_name = extracted_text[3].split("=")[0].strip()
+        creation_name = caller_name()
 
         recording = super().run()
         transferfunction = ImpulsiveResponse(self.excitation,

pytta/classes/signal.py

@@ -1750,7 +1750,7 @@ def _crossfade_spectruns(self, a, b, freqLims, freqVector):
         f1idx = np.where(freqVector <= f1)[0][-1]
         totalSamples = a.shape[0]
         xsamples = f1idx - f0idx
-        win = ss.hann(2*xsamples)
+        win = ss.windows.hann(2*xsamples)
 
         rightWin = win[xsamples-1:-1]
         fullRightWin = np.concatenate((np.ones(f0idx),