[FEATURE] EDC noise compensation methods and T20/T30 uncertainty

pytta/classes/analysis.py

@@ -957,7 +957,8 @@ def __init__(self, ir: SignalObj, nthOct: int = 1,
                  plotLundeby: bool = False,
                  bypassLundeby: bool = False,
                  suppressWarnings: bool = True,
-                 ircut: float = None, **kwargs):
+                 ircut: float = None,
+                 edcMethod: str = 'cutWithCorrection', **kwargs):
         _ir = ir.IR if type(ir) == ImpulsiveResponse else ir
         minBand = freq_to_band(minFreq, nthOct, 1000, 10)
         maxBand = freq_to_band(maxFreq, nthOct, 1000, 10)
@@ -966,6 +967,7 @@ def __init__(self, ir: SignalObj, nthOct: int = 1,
         self.ir = crop_IR(_ir, ircut)
         self._params, self.listEDC, self.fhSignal = self.estimate_energy_parameters(self.ir, self.bands, plotLundeby,
                                                        bypassLundeby, suppressWarnings,
+                                                       edcMethod=edcMethod,
                                                        nthOct=nthOct, minFreq=minFreq,
                                                        maxFreq=maxFreq)
         return
@@ -974,7 +976,9 @@ def __init__(self, ir: SignalObj, nthOct: int = 1,
     def estimate_energy_parameters(ir: SignalObj, bands: np.ndarray,
                                    plotLundeby: bool = False,
                                    bypassLundeby: bool = False,
-                                   suppressWarnings: bool = False, **kwargs):
+                                   suppressWarnings: bool = False,
+                                   edcMethod: str = 'cutWithCorrection',
+                                   **kwargs):
         """
         Estimate the Impulse Response energy parameters.
 
@@ -993,7 +997,9 @@ def estimate_energy_parameters(ir: SignalObj, bands: np.ndarray,
             A dict with parameters by name.
 
         """
-        listEDC, fhSignal = cumulative_integration(ir, bypassLundeby, plotLundeby, suppressWarnings, **kwargs)
+        listEDC, fhSignal, listLundeby = cumulative_integration(
+            ir, bypassLundeby, plotLundeby, suppressWarnings,
+            edcMethod=edcMethod, **kwargs)
         params = {}
         params['rms'] = fhSignal.rms()
         params['SPL'] = fhSignal.spl()
@@ -1005,6 +1011,13 @@ def estimate_energy_parameters(ir: SignalObj, bands: np.ndarray,
         params['EDT'] = reverberation_time('EDT', listEDC)
         params['T20'] = reverberation_time(20, listEDC)
         params['T30'] = reverberation_time(30, listEDC)
+        # peak-to-noise ratio: peak of the squared band response over the
+        # background noise level found by Lundeby (nan without it)
+        noise = np.array([p[3] if p is not None and p[3] > 0 else np.nan
+                          for p in listLundeby], dtype=float)
+        with np.errstate(divide='ignore', invalid='ignore'):
+            params['PSNR'] = 10 * np.log10(
+                np.max(fhSignal.timeSignal ** 2, axis=0) / noise)
         # self._params['BR'], self._params['TR'] = timbre_ratios(self.T20)
         return params, listEDC, fhSignal
 
@@ -1063,6 +1076,53 @@ def T30(self):
         """Reverberation time with 30 dB decay, by frequency `band`."""
         return self._params['T30']
 
+    @property
+    def PSNR(self):
+        """Peak-to-noise ratio by band, in dB: the peak of the squared band
+        response over the background noise level found by Lundeby; nan
+        with edcMethod='noCut'."""
+        return self._params['PSNR']
+
+    def timbre_ratios(self, param: str = 'T30'):
+        """
+        Bass ratio and treble ratio (Beranek) from the reverberation time
+        `param` of the octave bands: BR = (T125 + T250) / (T500 + T1000),
+        TR = (T2000 + T4000) / (T500 + T1000). nan when a band is missing.
+        """
+        bands, rt = np.asarray(self.bands), self._params[param]
+
+        def get(f):                       # band within 6 % of the nominal f
+            near = np.flatnonzero(np.isclose(bands, f, rtol=0.06))
+            return rt[near[0]] if near.size else np.nan
+        mid = get(500) + get(1000)
+        return (get(125) + get(250)) / mid, (get(2000) + get(4000)) / mid
+
+    def uncertainty(self, param: str = 'T30', numPositions: int = 1,
+                    numDecays: int = 10) -> np.ndarray:
+        """
+        Standard deviation of T20 or T30 by band, from ISO 3382-1:2009,
+        Equations (4) and (5): sigma(T) = k T sqrt((1 + c / n) / (N B T)),
+        with k, c = 0.88, 1.90 for T20 and 0.55, 1.52 for T30.
+
+        * param ('T30'), (str): 'T20' or 'T30';
+        * numPositions (1), (int): N, independent source and receiver
+          positions averaged;
+        * numDecays (10), (int): n, decays measured in each position; an
+          impulse response counts as n = 10 (clause 7.2).
+
+        B is the bandwidth of the band filter, fc (2**(1/(2 nthOct)) -
+        2**(-1/(2 nthOct))): 0.71 fc for octaves and 0.23 fc for thirds.
+        """
+        if param not in ('T20', 'T30'):
+            raise ValueError("param must be 'T20' or 'T30'.")
+        k, c = (0.88, 1.90) if param == 'T20' else (0.55, 1.52)
+        rt = self._params[param]
+        half = 2 ** (1 / (2 * self.nthOct))
+        B = np.asarray(self.bands) * (half - 1 / half)
+        with np.errstate(divide='ignore', invalid='ignore'):
+            return k * rt * np.sqrt((1 + c / numDecays)
+                                    / (numPositions * B * rt))
+
     # @property
     # def BR(self):
     #     """Reverberation time with 30 dB decay, by frequency `band`."""
@@ -1401,9 +1461,31 @@ def _Lundeby_correction(band, timeSignal, samplingRate, numSamples,
 # @njit
 def energy_decay_calculation(band, timeSignal, timeVector, samplingRate,
                              numSamples, numChannels, timeLength, bypassLundeby,
-                             suppressWarnings=True):
-    """Calculate the Energy Decay Curve."""
-    if bypassLundeby is False:
+                             suppressWarnings=True,
+                             edcMethod='cutWithCorrection'):
+    """
+    Calculate the Energy Decay Curve by backward integration, with one of
+    the noise compensations of ISO 3382-1 (Guski and Vorlaender, 2013):
+
+        * 'cutWithCorrection': truncate at the intersection of the decay
+          with the background noise (Lundeby et al., 1995) and add the
+          energy of the decay beyond it;
+        * 'justCut': truncate at the intersection, with no correction;
+        * 'noCut': integrate the whole response, noise included;
+        * 'subtractNoise': subtract the background noise level from the
+          squared response and integrate the whole of it (Chu, 1978).
+
+    bypassLundeby=True is the 'noCut' method.
+    """
+    methods = ('cutWithCorrection', 'justCut', 'noCut', 'subtractNoise')
+    if edcMethod not in methods:
+        raise ValueError(f'edcMethod must be one of {methods}.')
+    if bypassLundeby:
+        edcMethod = 'noCut'
+    if edcMethod == 'noCut':
+        lundebyParams = None
+        interIdx, lateRT, BGL = 0, 1, 0
+    else:
         lundebyParams = \
             _Lundeby_correction(band,
                                 timeSignal,
@@ -1414,22 +1496,18 @@ def energy_decay_calculation(band, timeSignal, timeVector, samplingRate,
                                 suppressWarnings=suppressWarnings)
         _, c1, interIdx, BGL = lundebyParams
         lateRT = -60/c1 if c1 != 0 else 0
-    else:
-        interIdx = 0
-        lateRT = 1
-
-    if interIdx == 0:
+    if interIdx == 0 or edcMethod == 'subtractNoise':
         interIdx = -1
-
     truncatedTimeSignal = timeSignal[:interIdx, 0]
     truncatedTimeVector = timeVector[:interIdx]
-
     if lateRT != 0.0:
-        if not bypassLundeby:
+        if edcMethod == 'cutWithCorrection':
             C = samplingRate*BGL*lateRT/(6*np.log(10))
         else:
             C = 0
         sqrInv = truncatedTimeSignal[::-1]**2
+        if edcMethod == 'subtractNoise':
+            sqrInv = sqrInv - BGL
         energyDecayFull = np.cumsum(sqrInv)[::-1] + C
         energyDecay = energyDecayFull/energyDecayFull[0]
     else:
@@ -1444,6 +1522,7 @@ def cumulative_integration(inputSignal,
                            bypassLundeby,
                            plotLundebyResults,
                            suppressWarnings=True,
+                           edcMethod='cutWithCorrection',
                            **kwargs):
     """Cumulative integration with proper corrections."""
 
@@ -1473,6 +1552,7 @@ def plot_lundeby():
                 freqRange=[kwargs['minFreq'],
                            kwargs['maxFreq']])[:, 1]
     listEDC = []
+    listLundeby = []
     for ch in range(hSignal.numChannels):
         signal = hSignal[ch]
         band = bands[ch]
@@ -1485,13 +1565,15 @@ def plot_lundeby():
                                      signal.numChannels,
                                      signal.timeLength,
                                      bypassLundeby,
-                                     suppressWarnings=suppressWarnings)
+                                     suppressWarnings=suppressWarnings,
+                                     edcMethod=edcMethod)
         listEDC.append((energyDecay, energyVector))
-        if plotLundebyResults:  # Placed here because Numba can't handle plots.
+        listLundeby.append(lundebyParams)
+        if plotLundebyResults and lundebyParams is not None:
             # plot_lundeby(band, timeVector, timeSignal,  samplingRate,
             #             lundebyParams)
             plot_lundeby()
-    return listEDC, hSignal
+    return listEDC, hSignal, listLundeby
 
 # @njit
 def reverb_time_regression(energyDecay, energyVector, upperLim, lowerLim):

tests/test_room_edc.py

@@ -0,0 +1,123 @@
+"""
+Tests for the energy decay curve methods and the reverberation time
+uncertainty of pytta.RoomAnalysis.
+
+The impulse response is a synthetic exponential decay of white noise with
+T = 1 s, over a steady background noise 45 dB below its start, so the
+methods of noise compensation of ISO 3382-1 (Guski and Vorlaender, 2013)
+can be told apart. Each test pins one rule.
+"""
+import os
+import unittest
+import warnings
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+import scipy.signal as ss
+
+import pytta
+
+FS = 48000
+T = 1.0
+
+
+def decay_with_noise():
+    rng = np.random.default_rng(0)
+    t = np.arange(int(2.5 * FS)) / FS
+    h = rng.standard_normal(t.size) * 10 ** (-3 * t / T) \
+        + 10 ** (-45 / 20) * rng.standard_normal(t.size)
+    return pytta.SignalObj(h, 'time', FS, signalType='energy')
+
+
+def analyse(ir, **kwargs):
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore')
+        return pytta.RoomAnalysis(ir, nthOct=1, minFreq=500., maxFreq=2000.,
+                                  **kwargs)
+
+
+class TestRoomEDC(unittest.TestCase):
+
+    @classmethod
+    def setUpClass(cls):
+        ir = decay_with_noise()
+        cls.T30 = {m: analyse(ir, edcMethod=m).T30 for m in
+                   ('cutWithCorrection', 'justCut', 'noCut', 'subtractNoise')}
+        cls.bypass = analyse(ir, bypassLundeby=True)
+        for m, v in cls.T30.items():
+            print(f'\n  {m:18s} T30 = {np.round(v, 3)}')
+
+    def test_compensated_methods_find_the_true_decay(self):
+        for m in ('cutWithCorrection', 'subtractNoise'):
+            self.assertTrue(np.allclose(self.T30[m], T, rtol=0.05), m)
+
+    def test_uncompensated_methods_are_biased_as_expected(self):
+        # The background noise flattens the tail of the full curve, and a
+        # truncation without correction makes the curve fall too fast.
+        self.assertTrue(np.all(self.T30['noCut'] > 1.10 * T))
+        self.assertTrue(np.all(self.T30['justCut']
+                               < self.T30['cutWithCorrection']))
+
+    def test_bypass_lundeby_is_the_no_cut_method(self):
+        self.assertTrue(np.array_equal(self.bypass.T30, self.T30['noCut']))
+
+    def test_uncertainty_follows_iso_3382_1(self):
+        # ISO 3382-1:2009, Equations (4) and (5), with the bandwidth of an
+        # octave filter B = 0.71 fc and n = 10 for an impulse response
+        # (clause 7.2).
+        ra = analyse(decay_with_noise())
+        fc = np.array([500., 1000., 2000.])
+        for param, k, c in (('T20', 0.88, 1.90), ('T30', 0.55, 1.52)):
+            rt = getattr(ra, param)
+            for N in (1, 4):
+                iso = k * rt * np.sqrt((1 + c / 10) / (N * 0.71 * fc * rt))
+                got = ra.uncertainty(param, numPositions=N)
+                self.assertTrue(np.allclose(got, iso, rtol=0.005), param)
+        with self.assertRaises(ValueError):
+            ra.uncertainty('EDT')
+
+
+if __name__ == '__main__':
+    unittest.main()
+
+
+class TestRoomExtras(unittest.TestCase):
+
+    def test_psnr_follows_the_noise_floor(self):
+        # The same decay over background noise 45 and 55 dB below its start:
+        # the peak-to-noise ratio of each band grows by 10 dB. Without the
+        # Lundeby estimate (noCut) there is no noise level, so no PSNR.
+        rng = np.random.default_rng(0)
+        t = np.arange(int(2.5 * FS)) / FS
+        decay = rng.standard_normal(t.size) * 10 ** (-3 * t / T)
+        noise = rng.standard_normal(t.size)
+        psnr = []
+        for floor in (45, 55):
+            ir = pytta.SignalObj(decay + 10 ** (-floor / 20) * noise, 'time',
+                                 FS, signalType='energy')
+            psnr.append(analyse(ir).PSNR)
+        print(f'\n  PSNR difference by band: {np.round(psnr[1] - psnr[0], 2)}')
+        self.assertTrue(np.allclose(psnr[1] - psnr[0], 10, atol=0.5))
+        self.assertTrue(np.all(np.isnan(analyse(ir, edcMethod='noCut').PSNR)))
+
+    def test_bass_and_treble_ratios(self):
+        # Octave bands decaying with T of 2.0, 1.8, 1.5, 1.5, 1.2 and 1.0 s
+        # from 125 Hz to 4 kHz: BR = (2.0 + 1.8) / 3.0, TR = (1.2 + 1.0) / 3.0.
+        rng = np.random.default_rng(1)
+        t = np.arange(4 * FS) / FS
+        h = np.zeros(t.size)
+        for fc, Tb in zip((125, 250, 500, 1000, 2000, 4000),
+                          (2.0, 1.8, 1.5, 1.5, 1.2, 1.0)):
+            sos = ss.butter(8, [fc / 2 ** 0.5, fc * 2 ** 0.5], 'bandpass',
+                            fs=FS, output='sos')
+            h += ss.sosfilt(sos, rng.standard_normal(t.size)) \
+                * 10 ** (-3 * t / Tb)
+        ra = pytta.RoomAnalysis(pytta.SignalObj(h, 'time', FS,
+                                                signalType='energy'),
+                                nthOct=1, minFreq=125., maxFreq=4000.)
+        BR, TR = ra.timbre_ratios('T30')
+        print(f'\n  BR {BR:.3f} (expected {3.8 / 3:.3f}), TR {TR:.3f} '
+              f'(expected {2.2 / 3:.3f})')
+        self.assertAlmostEqual(BR / (3.8 / 3), 1, delta=0.05)
+        self.assertAlmostEqual(TR / (2.2 / 3), 1, delta=0.05)