RAP Lateral Level based on energy ratio

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"""

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)