diff --git a/pyrato/parameters.py b/pyrato/parameters.py index 3b6c96ce..1e07725b 100644 --- a/pyrato/parameters.py +++ b/pyrato/parameters.py @@ -509,6 +509,85 @@ def _energy_ratio(limits, energy_decay_curve1, energy_decay_curve2): return energy_ratio +def late_lateral_sound_level(energy_decay_curve_free_field, + energy_decay_curve_room_lateral): + r""" + Calculate the late lateral sound level. + + The late lateral sound level :math:`L_\mathrm{J}` quantifies the strength + of late-arriving lateral sound energy. According to ISO 3382-1 [#isoLLJ]_, + it is defined as the level ratio between the late lateral sound energy + captured with a figure-of-eight microphone and the total sound energy of a + reference impulse response measured with an omnidirectional microphone at + a distance of 10 m in the free field. It is a measure of listener + envelopment. + + The parameter is defined as + + .. math:: + + L_\mathrm{J} = + 10 \log_{10} + \frac{ + \displaystyle \int_{0.08}^{\infty} p_\mathrm{L}^2(t)\,\mathrm{d}t + }{ + \displaystyle \int_{0}^{\infty} p_{10}^2(t)\,\mathrm{d}t + } + + where :math:`p_\mathrm{L}(t)` is the lateral sound pressure measured with a + figure-eight microphone whose zero axis is oriented towards the source, + and :math:`p_{10}(t)` is the instantaneous sound pressure of the + impulse response measured with an omnidirectional microphone + at 10 m distance in the free field. + + Using the energy decay curves of the reference response + :math:`e_{10}(t)` and the lateral response :math:`e_\mathrm{L}(t)`, + the parameter can be computed efficiently as + + .. math:: + + L_\mathrm{J} = + 10 \log_{10} + \frac{ + e_\mathrm{L}(0.08) + }{ + e_{10}(0) + }. + + Parameters + ---------- + energy_decay_curve_free_field : pyfar.TimeData + Energy decay curve of the reference free field impulse response + measured vwith an omnidirectional microphone at 10 m distance in the + free field. The EDC must start at time zero. + + energy_decay_curve_room_lateral : pyfar.TimeData + Energy decay curve of the room impulse response measured with a + figure-eight microphone oriented according to [#isoLLJ]_ + (zero axis pointing towards the source). The EDC must start at + time zero. + + Both EDCs must have identical ``signal.cshape``. + + Returns + ------- + Late Lateral Sound Level : numpy.ndarray + Late lateral sound level (:math:`L_\mathrm{J}`) in decibels, + shaped according to the channel shape of the input EDCs. + + References + ---------- + .. [#isoLLJ] ISO 3382-1, ISO 3382, Acoustics — Measurement of the + reverberation time of rooms with reference to other acoustical + parameters. + """ + + limits = np.array([0.0, np.inf, 0.08, np.inf]) + + return 10 * np.log10(_energy_ratio(limits, + energy_decay_curve_free_field, + energy_decay_curve_room_lateral)) + def sound_strength(energy_decay_curve_room, energy_decay_curve_free_field): r""" diff --git a/tests/test_parameters.py b/tests/test_parameters.py index a67ab3d4..b1856488 100644 --- a/tests/test_parameters.py +++ b/tests/test_parameters.py @@ -10,6 +10,7 @@ from pyrato.parameters import definition from pyrato.parameters import early_lateral_energy_fraction from pyrato.parameters import _energy_ratio +from pyrato.parameters import late_lateral_sound_level # parameter clarity tests @pytest.mark.parametrize( @@ -724,3 +725,88 @@ def test_JLF_for_exponential_decay_analytical(make_edc): expected = expected_lateral / expected_omni np.testing.assert_allclose(result, expected, atol=1e-5) + + +# late_lateral_level tests +@pytest.mark.parametrize( + ("energy", "expected_shape"), + [ + (np.linspace(1, 0, 1000), (1,)), + (np.linspace((1, 0.5), (0, 0), 1000).T, (2,)), + (np.arange(2 * 3 * 1000).reshape(2, 3, 1000), (2, 3)), + ], +) +def test_LJ_accepts_timedata_and_returns_correct_shape( + energy, expected_shape, make_edc, +): + """Return type and shape of pyfar.TimeData input for identical edcs.""" + edc = make_edc(energy=energy, sampling_rate=1000) + result = late_lateral_sound_level(edc, edc) + + assert isinstance(result, (float, np.ndarray)) + assert result.shape == expected_shape + assert result.shape == edc.cshape + +def test_LJ_returns_nan_for_zero_denominator_signal(): + """Correct return of NaN for division by zero signal.""" + edc_ref = pf.TimeData(np.zeros((1, 128)), np.arange(128) / 1000) + edc_lat = pf.TimeData(np.ones((1, 128)), np.arange(128) / 1000) + + result = late_lateral_sound_level(edc_ref, edc_lat) + assert np.isnan(result) + +def test_LJ_calculates_known_reference_value(make_edc): + """ + Construct simple deterministic EDCs: + e_10(0) = 1 + e_L(0.08) = 0.5 + Expected: + LJ = 10*log10(0.5) = -3.0103 dB. + """ + + pad = np.zeros(200) + + edc_ref = np.concatenate(([1.0], pad)) + + # 80 ms at 1 kHz sampling rate -> index 80 + edc_lateral = np.concatenate(( + np.zeros(80), + np.array([0.5]), + pad, + )) + + edc_ref = make_edc(energy=edc_ref, sampling_rate=1000) + edc_lateral = make_edc(energy=edc_lateral, + sampling_rate=1000, normalize=False) + + result = late_lateral_sound_level(edc_ref, edc_lateral) + + expected = 10 * np.log10(0.5) + np.testing.assert_allclose(result, expected, atol=1e-5) + +def test_LJ_for_exponential_decay_analytical(make_edc): + """ + LJ validation for analytical exponential decay: + e(t) = exp(-a*t) + LJ = 10*log10(e_L(0.08) / e_10(0)). + """ + + sampling_rate = 1000 + total_samples = 2000 + + # emulating free field impulse response + edc_ref = make_edc(rt=0.01, + sampling_rate=sampling_rate, + total_samples=total_samples) + + edc_lat = make_edc(rt=2.2, + sampling_rate=sampling_rate, + total_samples=total_samples) + + result = late_lateral_sound_level(edc_ref, edc_lat) + + a_lat = 13.8155 / 2.2 + + expected = 10 * np.log10(np.exp(-a_lat * 0.08)) + + np.testing.assert_allclose(result, expected, atol=1e-5)