RAP Clarity based on energy ratio

pyrato/parameters.py

@@ -104,15 +104,14 @@ def reverberation_time_linear_regression(
     else:
         return reverberation_times
 
-
 def clarity(energy_decay_curve, early_time_limit=80):
     r"""
     Calculate the clarity from the energy decay curve (EDC).
 
     The clarity parameter (C50 or C80) is defined as the ratio of early-to-late
     arriving energy in an impulse response and is a measure for how clearly
     speech or music can be perceived in a room. The early-to-late boundary is
-    typically set at 50 ms (C50) or 80 ms (C80) [#iso]_.
+    typically set at 50 ms (C50) or 80 ms (C80) [#isoClarity]_.
 
     Clarity is calculated as:
 
@@ -126,11 +125,14 @@ def clarity(energy_decay_curve, early_time_limit=80):
 
     where :math:`t_e` is the early time limit and :math:`p(t)` is the pressure
     of a room impulse response. Here, the clarity is efficiently computed
-    from the EDC :math:`e(t)` directly by:
+    from the EDC :math:`e(t)`
 
     .. math::
 
-        C_{t_e} = 10 \log_{10} \left( \frac{e(0)}{e(t_e)} - 1 \right).
+        C_{t_e} = 10 \log_{10} \frac{e(0) - e(t_e)}{e(t_e) - e(\infty)}
+                = 10 \log_{10} \left( \frac{e(0)}{e(t_e)} - 1 \right),
+
+    where :math:`e(\infty) = 0` by definition of the EDC.
 
     Parameters
     ----------
@@ -139,7 +141,7 @@ def clarity(energy_decay_curve, early_time_limit=80):
         (time-domain signal). The EDC must start at time zero.
     early_time_limit : float, optional
         Early time limit (:math:`t_e`) in milliseconds. Defaults to 80 (C80).
-        Typical values are 50 ms (C50) or 80 ms (C80) [#iso]_.
+        Typical values are 50 ms (C50) or 80 ms (C80) [#isoClarity]_.
 
     Returns
     -------
@@ -149,59 +151,50 @@ def clarity(energy_decay_curve, early_time_limit=80):
 
     References
     ----------
-    .. [#iso] ISO 3382, Acoustics — Measurement of the reverberation time of
-        rooms with reference to other acoustical parameters.
+    .. [#isoClarity] ISO 3382, Acoustics — Measurement of the reverberation
+        time of rooms with reference to other acoustical parameters.
 
     Examples
     --------
     Estimate the clarity from a real room impulse response filtered in
     octave bands:
 
-    >>> import numpy as np
     >>> import pyfar as pf
     >>> import pyrato as ra
+    ...
     >>> rir = pf.signals.files.room_impulse_response(sampling_rate=44100)
-    >>> rir = pf.dsp.filter.fractional_octave_bands(rir)
+    >>> rir = pf.dsp.filter.fractional_octave_bands(rir, num_fractions=1)
     >>> edc = ra.edc.energy_decay_curve_lundeby(rir)
-    >>> C80 = clarity(edc, early_time_limit=80)
+    ...
+    >>> C80 = ra.parameters.clarity(edc, early_time_limit=80)
+    >>> C80
+    ...     [[-55.57140506]
+    ...     [-11.75657677]
+    ...     [ -3.21150787]
+    ...     [  2.76276817]
+    ...     [  4.70786211]
+    ...     [  5.98148157]
+    ...     [  9.66764094]
+    ...     [  9.08687417]
+    ...     [ 14.14550646]
+    ...     [ 21.60048332]]
     """
 
-    # Check input type
-    if not isinstance(energy_decay_curve, pf.TimeData):
-        raise TypeError("Input must be a pyfar.TimeData object.")
     if not isinstance(early_time_limit, (int, float)):
         raise TypeError('early_time_limit must be a number.')
 
-    # Validate time range
-    if (early_time_limit > energy_decay_curve.signal_length * 1000) or (
-            early_time_limit <= 0):
-        raise ValueError(
-            "early_time_limit must be in the range of 0"
-            f"and {energy_decay_curve.signal_length * 1000}.",
-            )
-
-    # Raise error if TimeData is complex
-    if energy_decay_curve.complex:
-        raise ValueError(
-            "Complex-valued input detected. Clarity is"
-            "only defined for real TimeData.",
-        )
-
     # Convert milliseconds to seconds
-    early_time_limit_sec = early_time_limit / 1000.0
-
-    start_vals_energy_decay_curve = energy_decay_curve.time[..., 0]
+    early_time_limit_sec = early_time_limit / 1000
 
-    idx_early_time_limit = int(
-        energy_decay_curve.find_nearest_time(early_time_limit_sec))
-    vals_energy_decay_curve = \
-        energy_decay_curve.time[..., idx_early_time_limit]
-    vals_energy_decay_curve[vals_energy_decay_curve == 0] = np.nan
+    limits = np.array([early_time_limit_sec,
+                       np.inf,
+                       0.0,
+                       early_time_limit_sec])
 
-    clarity = start_vals_energy_decay_curve / vals_energy_decay_curve - 1
-    clarity_db = 10 * np.log10(clarity)
+    return 10*np.log10(_energy_ratio(limits,
+                                     energy_decay_curve,
+                                     energy_decay_curve))
 
-    return clarity_db
 
 def _energy_ratio(limits, energy_decay_curve1, energy_decay_curve2):
     r"""

tests/test_parameters.py

@@ -9,91 +9,71 @@
 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)
-
-    result = clarity(edc, early_time_limit=4)  # 4 ms
+@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_clarity_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 = clarity(edc)
     assert isinstance(result, (float, np.ndarray))
+    assert result.shape == expected_shape
     assert result.shape == edc.cshape
 
-def test_clarity_rejects_non_timedata_input():
-    invalid_input = np.array([1, 2, 3])
-    expected_error_message = "Input must be a pyfar.TimeData object."
-
-    with pytest.raises(TypeError, match=re.escape(expected_error_message)):
-        clarity(invalid_input)
-
 def test_clarity_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)):
         clarity(edc, invalid_time_limit)
 
-def test_clarity_rejects_complex_timedata():
-    complex_data = pf.TimeData(np.array([1+1j, 2+2j, 3+3j]),
-                               np.arange(3) / 1000, is_complex=True)
-    expected_error_message = "Complex-valued input detected. Clarity is"
-    "only defined for real TimeData."
-
-    with pytest.raises(ValueError, match=re.escape(expected_error_message)):
-        clarity(complex_data, early_time_limit=2)
-
-def test_clarity_rejects_invalid_time_range(make_edc):
-    energy = np.zeros(128)
-    edc = make_edc(energy=energy)
-    actual_signal_length_ms = edc.signal_length * 1000
-
-    # Test negative time limit
-    expected_error_message = "early_time_limit must be in the range of 0"
-    f"and {actual_signal_length_ms}."
-    with pytest.raises(ValueError, match=re.escape(expected_error_message)):
-        clarity(edc, early_time_limit=-1)
-
-    # Test time limit beyond signal length
-    with pytest.raises(ValueError, match=re.escape(expected_error_message)):
-        clarity(edc, early_time_limit=200000)
-
-def test_clarity_preserves_multichannel_shape(make_edc):
-    energy = np.ones((2,2,10)) / (1+np.arange(10))
-    edc = make_edc(energy=energy, sampling_rate=10)
-    output = clarity(edc, early_time_limit=80)
-    assert edc.cshape == output.shape
-
-
 def test_clarity_returns_nan_for_zero_signal():
+    """Correct return of NaN for zero signal."""
     edc = pf.TimeData(np.zeros((1, 128)), np.arange(128) / 1000)
     result = clarity(edc)
     assert np.isnan(result)
 
-
 def test_clarity_calculates_known_reference_value(make_edc):
-    # Linear decay → early_time_limit at 1/2 energy -> ratio = 1 -> 0 dB
-    edc_vals = np.array([1.0, 0.75, 0.5, 0.0])  # monotonic decay
+    """
+    Linear decay → early_time_limit at 1/2 energy -> ratio = 1 -> 0 dB.
+    Monotonic decay, 1 sample = 1ms.
+    """
+    edc_vals = np.array([1.0, 0.75, 0.5, 0.0])
     edc = make_edc(energy=edc_vals, sampling_rate=1000)
 
     result = clarity(edc, early_time_limit=2)
-    np.testing.assert_allclose(result, 0.0, atol=1e-6)
-
+    np.testing.assert_allclose(result, 0.0, atol=1e-5)
 
 def test_clarity_values_for_given_ratio(make_edc):
+    """Clarity 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)
+                   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 = clarity(edc, early_time_limit=80)
     clarity_value_db = 10 * np.log10(energy_early/energy_late)
-    npt.assert_allclose(result, clarity_value_db, atol=1e-6)
+    npt.assert_allclose(result, clarity_value_db, atol=1e-5)
 
 def test_clarity_for_exponential_decay(make_edc):
+    """Clarity validation for analytical solution from exponential decay."""
     rt60 = 2.0  # seconds
     sampling_rate = 1000
     total_samples = 2000