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..961fa3b --- /dev/null +++ b/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 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() 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()