RAP Definition based on energy ratio

pyrato/parameters.py

@@ -269,6 +269,95 @@ def early_lateral_energy_fraction(energy_decay_curve_omni,
                          energy_decay_curve_omni,
                          energy_decay_curve_lateral)
 
+
+
+def definition(energy_decay_curve, early_time_limit=50):
+    r"""
+    Calculate the definition from the energy decay curve (EDC).
+
+    The definition parameter (D50) is defined as the ratio of early-to-total
+    arriving energy in an impulse response and is a measure for how defined
+    speech or music can be perceived in a room. The early-to-total boundary is
+    typically set at 50 ms (D50) [#isoDefinition]_.
+
+    Definition is calculated as:
+
+    .. math::
+
+        D_{t_\mathrm{e}} = \frac{
+            \displaystyle \int_0^{t_\mathrm{e}} p^2(t) \, dt
+        }{
+            \displaystyle \int_{0}^{\infty} p^2(t) \, dt
+        }
+
+    where :math:`t_e` is the early time limit and :math:`p(t)` is the pressure
+    of a room impulse response. Here, the definition is efficiently computed
+    from the EDC :math:`e(t)` directly by:
+
+    .. math::
+
+        D_{t_\mathrm{e}} = \frac{e(0) - e(t_\mathrm{e})}{e(0) - e(\infty)}
+                = 1 - \left( \frac{e(t_\mathrm{e})}{e(0)} \right),
+
+    where :math:`e(\infty) = 0` by definition of the EDC.
+
+    Parameters
+    ----------
+    energy_decay_curve : pyfar.TimeData
+        Energy decay curve (EDC) of the room impulse response
+        (time-domain signal). The EDC must start at time zero.
+    early_time_limit : float, optional
+        Early time limit (:math:`t_\mathrm{e}`) in milliseconds. Defaults to
+        typical value 50 (D50) [#isoDefinition]_.
+
+    Returns
+    -------
+    definition : numpy.ndarray[float]
+        Definition index (early-to-total energy ratio),
+        shaped according to the channel shape of the input EDC.
+
+    References
+    ----------
+    .. [#isoDefinition] ISO 3382, Acoustics — Measurement of the reverberation
+        time of rooms with reference to other acoustical parameters.
+
+    Examples
+    --------
+    Estimate the defintion from a real room impulse response filtered in
+    octave bands:
+
+    >>> import pyfar as pf
+    >>> import pyrato
+    ...
+    >>> rir = pf.signals.files.room_impulse_response(sampling_rate=44100)
+    >>> rir = pf.dsp.filter.fractional_octave_bands(
+    >>>     rir, num_fractions=1, frequency_range=(125, 20e3))
+    >>> edc = pyrato.edc.energy_decay_curve_lundeby(rir)
+    ...
+    >>> D50 = pyrato.parameters.definition(edc, early_time_limit=50)
+    >>> D50
+    ...     [[0.25984852]
+    ...     [0.50208742]
+    ...     [0.66722359]
+    ...     [0.73528532]
+    ...     [0.87801455]
+    ...     [0.82757594]
+    ...     [0.86536142]
+    ...     [0.87374988]]
+    """
+
+    if not isinstance(early_time_limit, (int, float)):
+        raise TypeError('early_time_limit must be a number.')
+
+    # Convert milliseconds to seconds
+    early_time_limit_sec = early_time_limit / 1000
+
+    limits = np.array([0.0, np.inf, 0.0, early_time_limit_sec])
+
+    return _energy_ratio(limits,
+                         energy_decay_curve,
+                         energy_decay_curve)
+
 def _energy_ratio(limits, energy_decay_curve1, energy_decay_curve2):
     r"""
     Calculate the energy ratio for the time limits from two energy

tests/test_parameters.py

@@ -6,6 +6,7 @@
 import re
 
 from pyrato.parameters import clarity
+from pyrato.parameters import definition
 from pyrato.parameters import early_lateral_energy_fraction
 from pyrato.parameters import _energy_ratio
 
@@ -93,6 +94,99 @@ def test_clarity_for_exponential_decay(make_edc):
     expected_dB = 10 * np.log10(expected_ratio)
     np.testing.assert_allclose(result, expected_dB, atol=1e-6)
 
+
+# parameters definition tests
+@pytest.mark.parametrize(
+    ("energy", "expected_shape"),
+    [
+        # 1D single channel
+        (np.linspace(1, 0, 1000), (1,)),
+        # 2D two channels
+        (np.linspace((1, 0.5), (0, 0), 1000).T, (2,)),
+        # 3D multichannel (2x3 channels)
+        (np.arange(2 * 3 * 1000).reshape(2, 3, 1000), (2, 3)),
+    ],
+)
+def test_definition_accepts_timedata_and_returns_correct_shape(
+    energy, expected_shape, make_edc,
+):
+    """Test return shape and type of pyfar.TimeData input."""
+    edc = make_edc(energy=energy, sampling_rate=1000)
+    result = definition(edc)
+    assert isinstance(result, (float, np.ndarray))
+    assert result.shape == expected_shape
+    assert result.shape == edc.cshape
+
+def test_definition_rejects_non_numeric_early_time_limit(make_edc):
+    """Rejects non-number type early_time_limit."""
+    edc = make_edc()
+    invalid_time_limit = "not_a_number"
+    expected_error_message = "early_time_limit must be a number."
+
+    with pytest.raises(TypeError, match=re.escape(expected_error_message)):
+        definition(edc, invalid_time_limit)
+
+def test_definition_returns_zero_for_zero_signal(make_edc):
+    """
+    Definition must return 0.0 for a zero-energy signal.
+
+    The fixture clips all values to a small noise floor (min_energy),
+    so the EDC is flat. The denominator becomes min_energy - 0 (since
+    np.inf maps to zero in _energy_ratio), and the numerator becomes
+    min_energy - min_energy = 0. The result is therefore 0.0, not NaN.
+    """
+    edc = make_edc(energy=np.zeros(128), sampling_rate=1000)
+    result = definition(edc)
+    assert np.isclose(result, 0.0)
+
+def test_definition_calculates_known_reference_value(make_edc):
+    """
+    Linear decay → early_time_limit at index 2 -> ratio = 0.5
+    Monotonic decay, 1 sample = 1ms.
+    """
+    edc_vals = np.array([1.0, 0.75, 0.5, 0.25, 0.0])  # monotonic decay
+    edc = make_edc(energy=edc_vals[np.newaxis, :], sampling_rate=1000)
+
+    result = definition(edc, early_time_limit=2)
+    expected = 0.5
+    np.testing.assert_allclose(result, expected, atol=1e-5)
+
+def test_definition_for_exponential_decay(make_edc):
+    """Definition validation for analytical solution from exponential decay."""
+    rt60 = 2.0  # seconds
+    sampling_rate = 1000
+    total_samples = 2000
+    early_cutoff = 80  # ms
+
+    # Generate EDC
+    edc = make_edc(rt=rt60,
+                   sampling_rate=sampling_rate,
+                   total_samples=total_samples)
+    result = definition(edc, early_time_limit=early_cutoff)
+
+    # Analytical expected value
+    te = early_cutoff / 1000  # convert ms to seconds
+    a = 13.8155 / rt60
+    expected_ratio = 1- np.exp(-a * te)
+    np.testing.assert_allclose(result, expected_ratio, atol=1e-5)
+
+
+def test_definition_values_for_given_ratio(make_edc):
+    """Definition validation for a given ratio from analytical baseline."""
+    energy_early = 1
+    energy_late = .5
+    energy = np.zeros((3, 1000))
+    edc = make_edc(energy=energy,
+                   sampling_rate=1000,
+                   dynamic_range = 120.0)
+    edc.time[..., 10] = energy_early
+    edc.time[..., 100] = energy_late
+    edc = ra.edc.schroeder_integration(edc, is_energy=True)
+    edc = pf.dsp.normalize(edc, reference_method='max')
+    result = definition(edc, early_time_limit=80)
+    definition_value = energy_early/(energy_late+energy_early)
+    np.testing.assert_allclose(result, definition_value, atol=1e-5)
+
 # _energy_ratio tests
 @pytest.mark.parametrize(
     "limits",