[FEATURE] add late lateral sound level L_J to pytta.iso3382

pytta/__init__.py

@@ -108,6 +108,7 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import iso3382
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +131,7 @@
            'generate',
            'utils',
            'iso3741',
+           'iso3382',
            'rooms',
 
            # Functions

pytta/iso3382.py

@@ -0,0 +1,156 @@
+"""
+PyTTa ISO 3382-1 spatial measures:
+---------------------------------
+Early lateral energy fractions, late lateral sound level and interaural
+cross correlation coefficient of ISO 3382-1:2009, Annexes A and B, from
+impulse responses. Reverberation times, clarity, definition and centre time
+are in pytta.RoomAnalysis.
+
+    >>> JLF, JLFC = pytta.iso3382.lateral_fraction(omni, eight)
+    >>> LJ = pytta.iso3382.late_lateral_level(omni, eight, p10)
+    >>> IACC = pytta.iso3382.iacc(binaural, t1=0, t2=0.08)
+
+The time origin is the start of the impulse response: the first sample
+whose squared value is within 20 dB of the maximum, found on the
+broadband responses and kept for every band.
+"""
+import numpy as np
+
+from pytta.classes import OctFilter
+
+
+def _start(*responses):
+    energy = np.max([r ** 2 for r in responses], axis=0)
+    return int(np.argmax(energy >= energy.max() / 100))
+
+
+def _bands(signalObj, nthOct, minFreq, maxFreq):
+    """Channels of signalObj as arrays, (samples, bands) per channel."""
+    if nthOct is None:
+        return [signalObj.timeSignal[:, [ch]].astype(float)
+                for ch in range(signalObj.numChannels)]
+    filt = OctFilter(order=4, nthOct=nthOct,
+                     samplingRate=signalObj.samplingRate, minFreq=minFreq,
+                     maxFreq=maxFreq, refFreq=1000, base=10)
+    return [band.timeSignal.astype(float) for band in filt.filter(signalObj)]
+
+
+def _span(t1, t2, start, fs, n):
+    first = start + int(round(t1 * fs))
+    last = n if t2 is None else min(n, start + int(round(t2 * fs)))
+    return slice(first, last)
+
+
+def lateral_fraction(omni, eight, nthOct=1, minFreq=125., maxFreq=4000.):
+    """
+    Early lateral energy fractions JLF (A.14) and JLFC (A.15), by band:
+
+        JLF  = sum over 5 to 80 ms of pL**2 / sum over 0 to 80 ms of p**2
+        JLFC = sum over 5 to 80 ms of pL p  / sum over 0 to 80 ms of p**2
+
+    Parameters (default), (type):
+    -----------------------------
+        * omni (), (SignalObj): impulse response from an omnidirectional
+          microphone, one channel;
+        * eight (), (SignalObj): impulse response from a figure-of-eight
+          microphone at the same point, its null towards the source;
+        * nthOct (1), (int or None): bands per octave, or None for the
+          broadband value;
+        * minFreq, maxFreq (125, 4000), (float): band range, in Hz.
+
+    Return: (JLF, JLFC), arrays with one value per band, or floats when
+    nthOct is None.
+    """
+    fs = omni.samplingRate
+    start = _start(omni.timeSignal[:, 0])
+    p = _bands(omni, nthOct, minFreq, maxFreq)[0]
+    pL = _bands(eight, nthOct, minFreq, maxFreq)[0]
+    early = _span(0, 0.080, start, fs, len(p))
+    lateral = _span(0.005, 0.080, start, fs, len(p))
+    total = np.sum(p[early] ** 2, axis=0)
+    JLF = np.sum(pL[lateral] ** 2, axis=0) / total
+    JLFC = np.sum(pL[lateral] * p[lateral], axis=0) / total
+    if nthOct is None:
+        return float(JLF[0]), float(JLFC[0])
+    return JLF, JLFC
+
+
+def late_lateral_level(omni, eight, p10, distance=10., nthOct=1,
+                       minFreq=125., maxFreq=4000.):
+    """
+    Late lateral sound level LJ (A.16), by band, in dB:
+
+        LJ = 10 lg( integral from 80 ms to the end of pL**2
+                    / integral over the whole response of p10**2 )
+
+    Parameters (default), (type):
+    -----------------------------
+        * omni (), (SignalObj): impulse response from an omnidirectional
+          microphone at the measurement point, one channel, used only to
+          find the direct sound, the time origin;
+        * eight (), (SignalObj): impulse response from a figure-of-eight
+          microphone at the same point, its null towards the source;
+        * p10 (), (SignalObj): impulse response of the same source in free
+          field, from an omnidirectional microphone at the given distance,
+          with the same calibration as eight;
+        * distance (10), (float): distance of p10 in m; other than 10 m,
+          LpE,10 = LpE,d + 20 lg(d/10) (A.4);
+        * nthOct, minFreq, maxFreq: as in lateral_fraction.
+
+    Return: array with one value per band, or a float when nthOct is None.
+    """
+    start = _start(omni.timeSignal[:, 0])
+    pL = _bands(eight, nthOct, minFreq, maxFreq)[0]
+    ref = _bands(p10, nthOct, minFreq, maxFreq)[0]
+    late = _span(0.080, None, start, eight.samplingRate, len(pL))
+    free = _span(0, None, _start(p10.timeSignal[:, 0]), p10.samplingRate,
+                 len(ref))
+    LJ = 10 * np.log10(np.sum(pL[late] ** 2, axis=0) / eight.samplingRate
+                       / (np.sum(ref[free] ** 2, axis=0) / p10.samplingRate)) \
+        - 20 * np.log10(distance / 10)
+    return float(LJ[0]) if nthOct is None else LJ
+
+
+def late_lateral_average(LJ):
+    """
+    Frequency-averaged late lateral sound level LJ,avg (A.17), the energy
+    mean of LJ in the octave bands of 125 Hz, 250 Hz, 500 Hz and 1 000 Hz,
+    in dB.
+    """
+    LJ = np.asarray(LJ, dtype=float)
+    if LJ.shape != (4,):
+        raise ValueError('LJ must have the 4 octave bands from 125 Hz to '
+                         '1 000 Hz.')
+    return float(10 * np.log10(np.mean(10 ** (LJ / 10))))
+
+
+def iacc(binaural, t1=0., t2=None, nthOct=1, minFreq=125., maxFreq=4000.):
+    """
+    Interaural cross correlation coefficient (B.1, B.2) by band: the largest
+    magnitude of the normalized cross correlation of the left and right ear
+    responses over t1 to t2, for lags shorter than 1 ms.
+
+    Parameters (default), (type):
+    -----------------------------
+        * binaural (), (SignalObj): left ear in channel 0, right in 1;
+        * t1, t2 (0, None), (float): time window in s; None is the end of
+          the response. 0 to 0.08 gives the early IACC, 0.08 to None the
+          late one;
+        * nthOct, minFreq, maxFreq: as in lateral_fraction.
+
+    Return: array with one value per band, or a float when nthOct is None.
+    """
+    fs = binaural.samplingRate
+    start = _start(binaural.timeSignal[:, 0], binaural.timeSignal[:, 1])
+    left, right = _bands(binaural, nthOct, minFreq, maxFreq)
+    window = _span(t1, t2, start, fs, len(left))
+    pl, pr = left[window], right[window]
+    maxLag = int(np.ceil(0.001 * fs)) - 1          # |tau| < 1 ms
+    padded = np.pad(right, ((maxLag, maxLag), (0, 0)))
+    offset = window.start + maxLag
+    iacf = [np.sum(pl * padded[offset + tau:offset + tau + len(pl)], axis=0)
+            for tau in range(-maxLag, maxLag + 1)]
+    norm = np.sqrt(np.sum(pl ** 2, axis=0) * np.sum(pr ** 2, axis=0))
+    with np.errstate(divide='ignore', invalid='ignore'):
+        values = np.max(np.abs(iacf), axis=0) / norm
+    return float(values[0]) if nthOct is None else values

tests/test_iso3382_lj.py

@@ -0,0 +1,92 @@
+"""
+Tests for pytta.iso3382.late_lateral_level and late_lateral_average, the
+late lateral sound level LJ of ISO 3382-1:2009, A.2.5. The impulse responses
+are made of single reflections, so the expected values are exact to the
+single precision of SignalObj. Each test pins one rule.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+
+import pytta
+
+FS = 48000
+N = FS // 2
+
+
+def impulses(*arrivals, delay=0.0):
+    """Impulse response with (time in s, amplitude) arrivals."""
+    h = np.zeros(N)
+    for t, a in arrivals:
+        h[int(round((t + delay) * FS))] += a
+    return h
+
+
+class TestLateLateralLevel(unittest.TestCase):
+
+    # Direct sound at 0 (in the null of the figure of eight), lateral
+    # reflections at 30 ms and 70 ms (before the 80 ms of A.16), 85 ms and
+    # 150 ms. The free-field response at 10 m has amplitude 0.5.
+    omni = [(0, 1.0), (0.030, 0.4), (0.070, 0.5), (0.085, 0.3), (0.150, 0.2)]
+    eight = [(0.030, 0.4), (0.070, 0.5), (0.085, 0.3), (0.150, 0.2)]
+    expected = 10 * np.log10((0.3 ** 2 + 0.2 ** 2) / 0.5 ** 2)  # -2.840 dB
+
+    def level(self, delay=0.0, amplitude=0.5, **kwargs):
+        omni = pytta.SignalObj(impulses(*self.omni, delay=delay), 'time', FS)
+        eight = pytta.SignalObj(impulses(*self.eight, delay=delay), 'time', FS)
+        p10 = pytta.SignalObj(impulses((0.029, amplitude)), 'time', FS)
+        return pytta.iso3382.late_lateral_level(omni, eight, p10,
+                                                nthOct=None, **kwargs)
+
+    def test_lj_follows_a16_from_80_ms(self):
+        LJ = self.level()
+        print(f'\nLJ = {LJ:.4f} dB, difference {LJ - self.expected:.2e} dB')
+        self.assertAlmostEqual(LJ, self.expected, 6)
+
+    def test_80_ms_count_from_the_direct_sound(self):
+        # 10 ms of silence before the direct sound: counted from the first
+        # sample, the window would open at 70 ms and take the 70 ms
+        # reflection in.
+        LJ = self.level(delay=0.010)
+        print(f'\nLJ = {LJ:.4f} dB, difference {LJ - self.expected:.2e} dB')
+        self.assertAlmostEqual(LJ, self.expected, 6)
+
+    def test_reference_at_another_distance_follows_a4(self):
+        # The free-field response at 5 m is twice that at 10 m;
+        # LpE,10 = LpE,5 + 20 lg(5/10) gives the same LJ.
+        LJ = self.level(amplitude=1.0, distance=5.0)
+        print(f'\nLJ = {LJ:.4f} dB, difference {LJ - self.expected:.2e} dB')
+        self.assertAlmostEqual(LJ, self.expected, 6)
+
+    def test_one_value_per_octave_band(self):
+        # The late lateral part is the free-field response at half the
+        # amplitude, so LJ = 20 lg 0.5 in every band.
+        omni = pytta.SignalObj(impulses((0, 1.0), (0.2, 0.5)), 'time', FS)
+        eight = pytta.SignalObj(impulses((0.2, 0.5)), 'time', FS)
+        p10 = pytta.SignalObj(impulses((0.05, 1.0)), 'time', FS)
+        LJ = pytta.iso3382.late_lateral_level(omni, eight, p10, nthOct=1,
+                                              minFreq=125, maxFreq=4000)
+        print(f'\nLJ = {LJ} dB, difference {LJ - 20 * np.log10(0.5)} dB')
+        self.assertEqual(len(LJ), 6)
+        np.testing.assert_allclose(LJ, 20 * np.log10(0.5), atol=1e-3)
+
+
+class TestLateLateralAverage(unittest.TestCase):
+
+    def test_average_is_energetic_over_four_octaves(self):
+        # A.17: 10 lg(0.25 (3 x 10**-1 + 1)) = -4.881 dB, where the
+        # arithmetic mean would be -7.5 dB. Other band counts are refused.
+        LJavg = pytta.iso3382.late_lateral_average([-10, -10, -10, 0])
+        expected = 10 * np.log10(0.325)
+        print(f'\nLJ,avg = {LJavg:.4f} dB, difference '
+              f'{LJavg - expected:.2e} dB')
+        self.assertAlmostEqual(LJavg, expected, 12)
+        with self.assertRaises(ValueError):
+            pytta.iso3382.late_lateral_average([-10, -10, -10, 0, -3, -5])
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_iso3382_spatial.py

@@ -0,0 +1,91 @@
+"""
+Tests for pytta.iso3382.lateral_fraction and pytta.iso3382.iacc, the early
+lateral energy fractions (ISO 3382-1:2009, A.2.4) and the interaural cross
+correlation coefficient (Annex B). The impulse responses are made of single
+reflections, so the expected values are exact. Each test pins one rule.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+
+import pytta
+
+FS = 48000
+N = FS // 2
+
+
+def impulses(*arrivals, delay=0.0):
+    """Impulse response with (time in s, amplitude) arrivals."""
+    h = np.zeros(N)
+    for t, a in arrivals:
+        h[int(round((t + delay) * FS))] += a
+    return h
+
+
+class TestLateralFraction(unittest.TestCase):
+
+    # Direct sound at 0 (in the null of the figure of eight), lateral
+    # reflections at 3 ms (before the 5 ms of A.14), 20 ms and 40 ms (the
+    # last 60 degrees off the lateral axis, so cos = 0.5), and one at 100 ms,
+    # after the 80 ms window.
+    omni = [(0, 1.0), (0.003, 0.3), (0.020, 0.5), (0.040, 0.4), (0.100, 0.6)]
+    eight = [(0.003, 0.3), (0.020, 0.5), (0.040, 0.2), (0.100, 0.6)]
+    energy = 1 + 0.3 ** 2 + 0.5 ** 2 + 0.4 ** 2
+
+    def fractions(self, delay=0.0):
+        omni = pytta.SignalObj(impulses(*self.omni, delay=delay), 'time', FS)
+        eight = pytta.SignalObj(impulses(*self.eight, delay=delay), 'time', FS)
+        return pytta.iso3382.lateral_fraction(omni, eight, nthOct=None)
+
+    def test_jlf_and_jlfc_follow_a14_and_a15(self):
+        JLF, JLFC = self.fractions()
+        self.assertAlmostEqual(JLF, (0.5 ** 2 + 0.2 ** 2) / self.energy, 6)
+        self.assertAlmostEqual(JLFC, (0.5 * 0.5 + 0.2 * 0.4) / self.energy, 6)
+
+    def test_time_starts_at_the_direct_sound(self):
+        # Silence before the direct sound does not change the fractions.
+        self.assertTrue(np.allclose(self.fractions(delay=0.010),
+                                    self.fractions(), atol=1e-7))
+
+
+class TestIACC(unittest.TestCase):
+
+    def iacc(self, left, right, **kwargs):
+        binaural = pytta.SignalObj(np.column_stack((left, right)), 'time', FS)
+        return pytta.iso3382.iacc(binaural, nthOct=None, **kwargs)
+
+    def test_lag_range_is_below_1_ms(self):
+        # The same arrivals 0.5 ms apart correlate fully; 2 ms apart, not
+        # at all. A sign change does not matter (B.2 takes |IACF|).
+        left = impulses((0, 1.0), (0.012, 0.7), (0.031, -0.4))
+        self.assertAlmostEqual(self.iacc(left, left), 1, 6)
+        self.assertAlmostEqual(self.iacc(left, -left), 1, 6)
+        self.assertAlmostEqual(self.iacc(left, np.roll(left, 24)), 1, 6)
+        self.assertAlmostEqual(self.iacc(left, np.roll(left, 96)), 0, 6)
+
+    def test_time_window_of_early_and_late_parts(self):
+        # Identical early parts and late parts 3 ms apart: IACC is 1 for
+        # 0 to 80 ms and 0 for 80 ms to the end.
+        early = [(0, 1.0), (0.015, 0.5)]
+        late = [(0.150, 0.3)]
+        left = impulses(*early, *late)
+        right = impulses(*early, *[(t + 0.003, a) for t, a in late])
+        self.assertAlmostEqual(self.iacc(left, right, t1=0, t2=0.08), 1, 6)
+        self.assertAlmostEqual(self.iacc(left, right, t1=0.08), 0, 6)
+
+    def test_one_value_per_octave_band(self):
+        rng = np.random.default_rng(0)
+        t = np.arange(N) / FS
+        h = rng.standard_normal(N) * np.exp(-t / 0.1)
+        binaural = pytta.SignalObj(np.column_stack((h, h)), 'time', FS)
+        values = pytta.iso3382.iacc(binaural, nthOct=1, minFreq=125,
+                                    maxFreq=4000)
+        self.assertEqual(len(values), 6)
+        self.assertTrue(np.allclose(values, 1, atol=1e-6))
+
+
+if __name__ == '__main__':
+    unittest.main()