Energy Ratio base function for room-acoustic parameters

pyrato/parameters.py

@@ -202,3 +202,111 @@ def clarity(energy_decay_curve, early_time_limit=80):
     clarity_db = 10 * np.log10(clarity)
 
     return clarity_db
+
+def _energy_ratio(limits, energy_decay_curve1, energy_decay_curve2):
+    r"""
+    Calculate the energy ratio for the time limits from two energy
+    decay curves (EDC).
+
+    A variety of room-acoustic parameters are defined by energy ratios derived
+    from one or two time-domain Energy Decay Curves (EDCs). These parameters
+    distinguish between time regions using the four provided limits, and some,
+    such as Strength (:math:`G`), Early lateral sound (:math:`J_\mathrm{LF}`),
+    and Late lateral sound (:math:`L_J`), require EDCs obtained from different
+    impulse-response measurements [#iso]_.
+
+    Energy-Ratio is calculated as:
+    .. math::
+        ER = \frac{
+            \displaystyle \int_{lim3}^{lim4} p_2^2(t) \, dt
+        }{
+            \displaystyle \int_{lim1}^{lim2} p_1^2(t) \, dt
+        }
+    where :math:`[lim1, ..., lim4]` are the time limits and :math:`p(t)` is the
+    pressure of a room impulse response. Here, the energy ratio is
+    efficiently computed from the EDC :math:`e(t)` directly by:
+    .. math::
+        ER = \frac{
+            \displaystyle e_2(lim3) - e_2(lim4)
+        }{
+            \displaystyle e_1(lim1) - e_1(lim2)
+        }.
+
+    Parameters
+    ----------
+    limits : np.ndarray, list or tuple
+        Four time limits (:math:`t_e`) in seconds, shape (4,)
+        in ascending order.
+    energy_decay_curve1 : pyfar.TimeData
+        Energy decay curve 1 (EDC1) of the room impulse response
+        (time-domain signal). The EDC must start at time zero.
+    energy_decay_curve2 : pyfar.TimeData
+        Energy decay curve 2 (EDC2) of the room impulse response
+        (time-domain signal). The EDC must start at time zero.
+
+    Returns
+    -------
+    energy ratio : numpy.ndarray[float]
+        energy-ratio index (early-to-late energy ratio),
+        shaped according to the channel shape of the input EDC.
+
+    References
+    ----------
+    .. [#iso] ISO 3382, Acoustics — Measurement of the reverberation time of
+        rooms with reference to other acoustical parameters.
+    """
+
+    # Check input type
+    if not isinstance(energy_decay_curve1, pf.TimeData):
+        raise TypeError(
+            "energy_decay_curve1 must be a pyfar.TimeData or derived object.")
+    if not isinstance(energy_decay_curve2, pf.TimeData):
+        raise TypeError(
+            "energy_decay_curve2 must be a pyfar.TimeData or derived object.")
+
+    if isinstance(limits, (list, tuple)):
+        limits = np.asarray(limits)
+
+    if not isinstance(limits, np.ndarray):
+        raise TypeError("limits must be a numpy ndarray, list, or tuple.")
+
+    # Check shape
+    if limits.shape != (4,):
+        raise ValueError(
+            "limits must have shape (4,), " \
+            "containing [lim1, lim2, lim3, lim4].",
+            )
+
+    # Check if limits are within valid time range
+    if (
+        np.any(limits[0:2] < 0)
+        or np.any(limits[0:2] > energy_decay_curve1.signal_length)
+    ):
+        raise ValueError(
+            f"limits[0:2] must be between 0 and "
+            f"{energy_decay_curve1.signal_length} seconds.",
+        )
+    if (
+        np.any(limits[2:4] < 0)
+        or np.any(limits[2:4] > energy_decay_curve2.signal_length)
+    ):
+        raise ValueError(
+            f"limits[2:4] must be between 0 and "
+            f"{energy_decay_curve2.signal_length} seconds.",
+        )
+
+    limits_energy_decay_curve1_idx = energy_decay_curve1.find_nearest_time(
+        limits[0:2])
+    limits_energy_decay_curve2_idx = energy_decay_curve2.find_nearest_time(
+        limits[2:4])
+
+    numerator = np.diff(
+        energy_decay_curve2.time[..., limits_energy_decay_curve2_idx],
+        axis=-1)[..., 0]
+    denominator = np.diff(
+        energy_decay_curve1.time[..., limits_energy_decay_curve1_idx],
+        axis=-1)[..., 0]
+
+    energy_ratio = numerator / denominator
+
+    return energy_ratio

tests/test_parameters.py

@@ -2,10 +2,13 @@
 import pytest
 import pyfar as pf
 import pyrato as ra
-from pyrato.parameters import clarity
 import numpy.testing as npt
 import re
 
+from pyrato.parameters import clarity
+from pyrato.parameters import _energy_ratio
+
+# parameter clarity tests
 def test_clarity_accepts_timedata_returns_correct_type(make_edc):
     energy = np.concatenate(([1, 1, 1, 1], np.zeros(124)))
     edc = make_edc(energy=energy)
@@ -108,3 +111,179 @@ def test_clarity_for_exponential_decay(make_edc):
     expected_ratio = np.exp(a * te) - 1
     expected_dB = 10 * np.log10(expected_ratio)
     np.testing.assert_allclose(result, expected_dB, atol=1e-6)
+
+# _energy_ratio tests
+@pytest.mark.parametrize(
+    "limits",
+    [[0.0, 0.001, 0.0, 0.005],
+    (0.0, 0.001, 0.0, 0.005),
+    np.array([0.0, 0.001, 0.0, 0.005])],
+)
+def test_energy_ratio_accepts_timedata_and_limits_returns_correct_shape(limits,
+                                                                        make_edc):
+    """Test return shape of pyfar.TimeData and accepted limits input types."""
+    energy = np.linspace((1,0.5),(0,0),1000).T
+    edc = make_edc(energy=energy, sampling_rate=1000)
+    result = _energy_ratio(limits, edc, edc)
+    assert isinstance(result, np.ndarray)
+    assert result.shape == edc.cshape
+
+def test_energy_ratio_rejects_non_timedata_input():
+    """Reject wrong input type of EDC."""
+    invalid_input = np.arange(10)
+    limits = np.array([0.0, 0.001, 0.0, 0.005])
+    expected_message = "energy_decay_curve1 must be a pyfar.TimeData " \
+    "or derived object."
+    with pytest.raises(TypeError, match=expected_message):
+        _energy_ratio(limits, invalid_input, invalid_input)
+
+def test_energy_ratio_rejects_if_second_edc_is_not_timedata(make_edc):
+    """Reject if second EDC is of wrong type."""
+    edc = make_edc(energy=np.linspace(1, 0, 10), sampling_rate=1000)
+    limits = np.array([0.0, 0.001, 0.0, 0.005])
+    with pytest.raises(
+        TypeError,
+        match="energy_decay_curve2 must be a pyfar.TimeData",
+    ):
+        _energy_ratio(limits, edc, "invalid_type")
+
+def test_energy_ratio_rejects_wrong_shape_limits(make_edc):
+    """Limits array wrong shape."""
+    edc = make_edc(energy=np.linspace(1, 0, 10), sampling_rate=1000)
+    wrong_shape_limits = np.array([0.0, 0.001, 0.005])  # Only 3 elements
+    with pytest.raises(ValueError, match="limits must have shape"):
+        _energy_ratio(wrong_shape_limits, edc, edc)
+
+def test_energy_ratio_rejects_wrong_type_limits(make_edc):
+    """Rejects wrong limits type correctly."""
+    edc = make_edc(energy=np.linspace(1, 0, 10), sampling_rate=1000)
+    wrong_type_limits = "3, 2, 0.5, 1"  # string
+    with pytest.raises(TypeError,
+                       match="limits must be a numpy ndarray"):
+        _energy_ratio(wrong_type_limits, edc, edc)
+
+def test_energy_ratio_computes_known_ratio_correctly(make_edc):
+    """
+    If EDC is linear, energy ratio should be 1.
+
+    numerator = e(lim3)-e(lim4) = (1.0 - 0.75) = 0.25
+    denominator = e(lim1)-e(lim2) = (0.75 - 0.5) = 0.25
+    ratio = 1
+    """
+    edc_vals = np.array([1.0, 0.75, 0.5, 0.25])
+    edc = make_edc(energy=edc_vals, sampling_rate=1000)
+
+    # For linear EDC:
+    limits = np.array([0.0, 0.001, 0.001, 0.002])
+    result = _energy_ratio(limits, edc, edc)
+    npt.assert_allclose(result, 1.0, atol=1e-12)
+
+@pytest.mark.parametrize(
+    "energy",
+    [
+        # 1D, einzelner Kanal
+        np.linspace(1, 0, 10),
+        # 2D, zwei Kanäle
+        np.stack([
+            np.linspace(1, 0, 10),
+            np.linspace(0.5, 0, 10),
+        ]),
+        # 3D – deterministic 2×3×4 “multichannel” structure
+        np.arange(2 * 3 * 4).reshape(2, 3, 4),
+    ],
+)
+def test_energy_ratio_preserves_multichannel_shape_correctly(energy, make_edc):
+    """Preserves any multichannel shape (1,), (2,), (2,3,)."""
+    edc = make_edc(energy=energy, sampling_rate=1000)
+    limits = np.array([0.0, 0.001, 0.0, 0.003])
+
+    result = _energy_ratio(limits, edc, edc)
+
+    # Ergebnis muss exakt dieselbe Kanalform haben wie das EDC
+    assert result.shape == edc.cshape
+
+def test_energy_ratio_returns_nan_for_zero_denominator(make_edc):
+    """If denominator e(lim1)-e(lim2)=0, expect NaN (invalid ratio)."""
+    energy = np.ones(10)
+    edc = make_edc(energy=energy, sampling_rate=1000)
+    limits = np.array([0.0, 0.001, 0.002, 0.003])
+    result = _energy_ratio(limits, edc, edc)
+    assert np.isnan(result)
+
+def test_energy_ratio_matches_reference_case(make_edc):
+    r"""
+    Analytical reference:
+    EDC = exp(-a*t). For exponential decay, ratio known analytically from
+    .. math::
+        ER = \frac{
+            \displaystyle e(lim3) - e(lim4)
+        }{
+            \displaystyle e(lim1) - e(lim2)
+        }.
+    where :math:`[lim1, ..., lim4]` are the time limits and here
+    the energy ratio is efficiently computed from the EDC :math:`e(t)'.
+    """
+    sampling_rate = 1000
+    a = 13.8155  # decay constant
+    times = np.arange(1000) / sampling_rate
+    edc_vals = np.exp(-a * times)
+    edc = make_edc(energy=edc_vals, sampling_rate=sampling_rate)
+
+    limits = np.array([0.0, 0.02, 0.0, 0.05])
+    lim1, lim2, lim3, lim4 = limits
+
+    analytical_ratio = (
+        (np.exp(-a*lim3) - np.exp(-a*lim4)) /
+        (np.exp(-a*lim1) - np.exp(-a*lim2))
+    )
+
+    result = _energy_ratio(limits, edc, edc)
+    npt.assert_allclose(result, analytical_ratio, atol=1e-8)
+
+def test_energy_ratio_works_with_two_different_edcs(make_edc):
+    """
+    Energy ratio between two different EDCs should compute distinct ratio.
+    """
+    edc1 = make_edc(energy=np.linspace(1, 0, 10), sampling_rate=1000)
+    edc2 = make_edc(energy=np.linspace(1, 0, 10) ** 2, sampling_rate=1000)
+
+    limits = np.array([0.0, 0.002, 0.0, 0.004])
+    # Expect a ratio != 1 because edc2 decays faster
+    ratio = _energy_ratio(limits, edc1, edc2)
+    assert not np.allclose(ratio, 1.0)
+
+def test_energy_ratio_rejects_limits_outside_time_range(make_edc):
+    """Limits outside valid time range are rejected."""
+    edc1 = make_edc(energy=np.linspace(1, 0, 100), sampling_rate=1000)
+    edc2 = make_edc(energy=np.linspace(1, 0, 100), sampling_rate=1000)
+    max_time = edc1.times[-1]
+
+    # Test negative limit
+    limits_negative = np.array([-0.01, 0.02, 0.02, 0.05])
+    with pytest.raises(
+        ValueError,
+        match=r"limits\[0:2\] must be between 0 and",
+    ):
+        _energy_ratio(limits_negative, edc1, edc2)
+
+    # Test limit beyond signal length
+    limits_too_large = np.array([0.0, 0.02, 0.02, max_time + 0.01])
+    with pytest.raises(
+        ValueError,
+        match=r"limits\[2:4\] must be between 0 and",
+    ):
+        _energy_ratio(limits_too_large, edc1, edc2)
+
+def test_energy_ratio_handles_different_edc_lengths(make_edc):
+    """Validation uses the shorter EDC's time range."""
+    edc1 = make_edc(energy=np.linspace(1, 0, 100), sampling_rate=1000)
+    edc2 = make_edc(energy=np.linspace(1, 0, 50), sampling_rate=1000)
+
+    # Limit valid for edc1 but not edc2
+    limits = np.array([0.0, 0.02, 0.02, 0.06])  # 0.06s > edc2.times[-1]
+
+    with pytest.raises(
+        ValueError,
+        match=r"limits\[2:4\] must be between 0 and",
+    ):
+        _energy_ratio(limits, edc1, edc2)