From 7adcc93e495c33c85ba28f1c389ba238d9f63f2f Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 20:39:15 +0100 Subject: [PATCH 1/2] [FEATURE] add pytta.iso3382 with lateral energy fractions and IACC pytta.iso3382.lateral_fraction gives the early lateral energy fractions JLF and JLFC of ISO 3382-1:2009 (A.14, A.15) from the impulse responses of an omnidirectional and a figure-of-eight microphone, and pytta.iso3382.iacc the interaural cross correlation coefficient (B.1, B.2) of a binaural response over any time window. Both work by octave or fractional octave band, or broadband. The time origin is the start of the response, the first sample within 20 dB of the maximum. Tests use responses made of single reflections, so JLF, JLFC and IACC are exact, and pin the 5 ms and 80 ms limits, the lags below 1 ms, the early and late windows and the band output. Co-Authored-By: Claude Fable 5.1 --- pytta/__init__.py | 2 + pytta/iso3382.py | 105 ++++++++++++++++++++++++++++++++++ tests/test_iso3382_spatial.py | 91 +++++++++++++++++++++++++++++ 3 files changed, 198 insertions(+) create mode 100644 pytta/iso3382.py create mode 100644 tests/test_iso3382_spatial.py diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..b768a9d 100644 --- a/pytta/__init__.py +++ b/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 diff --git a/pytta/iso3382.py b/pytta/iso3382.py new file mode 100644 index 0000000..5be9467 --- /dev/null +++ b/pytta/iso3382.py @@ -0,0 +1,105 @@ +""" +PyTTa ISO 3382-1 spatial measures: +--------------------------------- +Early lateral energy fractions 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) + >>> 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 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 diff --git a/tests/test_iso3382_spatial.py b/tests/test_iso3382_spatial.py new file mode 100644 index 0000000..bd28cd9 --- /dev/null +++ b/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() From f3ff054078725bcef8ac0e12c6f44c0f1d843b75 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 21:56:35 +0100 Subject: [PATCH 2/2] [FEATURE] add late lateral sound level L_J to pytta.iso3382 pytta.iso3382.late_lateral_level(omni, eight, p10) gives the late lateral sound level LJ of ISO 3382-1:2009, A.2.5, Equation (A.16), by band: the energy of the figure-of-eight response from 80 ms to the end over the energy of the free-field response of the same source at 10 m. The time origin is the direct sound found on the omnidirectional response, as in lateral_fraction. A free-field response at another distance d can be used with distance=d, which applies LpE,10 = LpE,d + 20 lg(d/10) of A.4. pytta.iso3382.late_lateral_average(LJ) gives LJ,avg of Equation (A.17), the energy mean over the octave bands from 125 Hz to 1000 Hz, which Table A.1 asks for in place of the arithmetic mean used for the other measures. Co-Authored-By: Claude Fable 5.1 --- pytta/iso3382.py | 57 +++++++++++++++++++++++-- tests/test_iso3382_lj.py | 92 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+), 3 deletions(-) create mode 100644 tests/test_iso3382_lj.py diff --git a/pytta/iso3382.py b/pytta/iso3382.py index 5be9467..961fa3b 100644 --- a/pytta/iso3382.py +++ b/pytta/iso3382.py @@ -1,11 +1,13 @@ """ PyTTa ISO 3382-1 spatial measures: --------------------------------- -Early lateral energy fractions 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. +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 @@ -73,6 +75,55 @@ def lateral_fraction(omni, eight, nthOct=1, minFreq=125., maxFreq=4000.): 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 diff --git a/tests/test_iso3382_lj.py b/tests/test_iso3382_lj.py new file mode 100644 index 0000000..9af4c96 --- /dev/null +++ b/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()