RAP center time

pyrato/parameters.py

@@ -6,6 +6,7 @@
 import numpy as np
 import pyfar as pf
 import warnings
+from scipy.integrate import trapezoid
 
 def reverberation_time_linear_regression(
         energy_decay_curve, T='T20', return_intercept=False):
@@ -1082,3 +1083,81 @@ def _sti_calc(mtf):
 
     # limit STI to 1 (IEC 60268-16:2020, Section A.2.1)
     return min(sti, 1.0)
+
+
+def center_time(energy_decay_curve):
+    r"""
+    Calculate the room-acoustic center time (:math:`T_s`).
+
+    The center time :math:`T_s` is the time of the centroid of the squared
+    impulse response. It quantifies the balance between early and late
+    sound energy [#isoTs]_.
+
+    The parameter is defined as
+
+    .. math::
+
+        T_s =
+        \frac{
+            \displaystyle \int_{0}^{\infty} t \cdot p^2(t)\,\mathrm{d}t
+        }{
+            \displaystyle \int_{0}^{\infty} p^2(t)\,\mathrm{d}t
+        }
+
+    where :math:`p(t)` is the room impulse response sound pressure.
+
+    Using the energy decay curve :math:`e(t)`, the parameter can be
+    computed efficiently via the Schroeder backward integral identity
+    :math:`e(t) = \int_t^\infty p^2(\tau)\,\mathrm{d}\tau` [#schroeder]_
+    and Fubini's theorem as
+
+    .. math::
+
+        T_s =
+        \frac{
+            \displaystyle \int_{0}^{\infty} e(t)\,\mathrm{d}t
+        }{
+            e(0)
+        }.
+
+    Parameters
+    ----------
+    energy_decay_curve : pyfar.TimeData
+        Energy decay curve of the room impulse response.
+
+    Returns
+    -------
+    center_time : numpy.ndarray
+        Center time (:math:`T_s`) in seconds,
+        shaped according to the channel shape of the input EDC.
+
+    References
+    ----------
+    .. [#isoTs] ISO 3382, Acoustics — Measurement of the reverberation
+        time of rooms with reference to other acoustical parameters.
+    .. [#schroeder] M. R. Schroeder, "New Method of Measuring Reverberation
+        Time," J. Acoust. Soc. Am., vol. 37, no. 3, pp. 409-412, 1965.
+        doi:10.1121/1.1909343
+    """
+
+    if not isinstance(energy_decay_curve, pf.TimeData):
+        raise TypeError(
+            "energy_decay_curve must be a pyfar.TimeData or derived object.")
+
+    if not np.isclose(energy_decay_curve.times[0], 0.0):
+        raise ValueError("energy_decay_curve must start at time zero.")
+
+    if np.any(energy_decay_curve.time[..., 0] == 0):
+        raise ValueError(
+            "Initial energy of energy_decay_curve must not be zero.")
+
+    initial_energy = energy_decay_curve.time[..., 0]
+    # NaN values in the tail (e.g. from Lundeby/Chu truncation) are treated as
+    # zero, since the EDC is assumed to converge to zero.
+    edc = np.nan_to_num(energy_decay_curve.time, nan=0.0)
+    center_time = (
+        trapezoid(edc, energy_decay_curve.times, axis=-1)
+        / initial_energy
+    )
+
+    return center_time

tests/test_parameters.py

@@ -16,6 +16,7 @@
 from pyrato.parameters import modulation_transfer_function
 from pyrato.parameters import _sti_calc
 from pyrato.parameters import _ambient_noise_correction
+from pyrato.parameters import center_time
 # parameter clarity tests
 @pytest.mark.parametrize(
     ("energy", "expected_shape"),
@@ -1391,3 +1392,115 @@ def test_sti_ir_level_snr():
         sti_test = speech_transmission_index_indirect(
             ir, rir_type="acoustical", level=level, snr=snr)
     np.testing.assert_allclose(sti_test, sti_expected, atol=0.07)
+
+@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_center_time_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 = center_time(edc)
+    assert isinstance(result, np.ndarray)
+    assert result.shape == expected_shape
+    assert result.shape == edc.cshape
+
+def test_center_time_rejects_non_timedata():
+    """Reject wrong input type."""
+    with pytest.raises(TypeError,
+                       match="energy_decay_curve must be a pyfar.TimeData"):
+        center_time(np.ones(100))
+
+def test_center_time_rejects_edc_not_starting_at_zero():
+    """Reject EDC whose time axis does not start at zero."""
+    edc = pf.TimeData(np.ones((1, 100)), np.arange(1, 101) / 1000)
+    with pytest.raises(ValueError, match="must start at time zero"):
+        center_time(edc)
+
+def test_center_time_rejects_zero_initial_energy():
+    """Reject EDC with zero initial energy (would cause division by zero)."""
+    edc = pf.TimeData(np.zeros((1, 100)), np.arange(100) / 1000)
+    with pytest.raises(ValueError, match="Initial energy"):
+        center_time(edc)
+
+def test_center_time_accepts_non_uniform_time_spacing():
+    """Accept EDC with non-uniform time spacing (trapezoid integration)."""
+    # monotonically increasing but not uniform
+    times = np.concatenate([[0, 0.001, 0.003], np.arange(4, 101) / 1000])
+    edc = pf.TimeData(np.ones((1, len(times))), times)
+    result = center_time(edc)
+    assert np.all(np.isfinite(result))
+
+def test_center_time_non_uniform_spacing_analytical():
+    r"""center_time() is correct for non-uniform spacing with known solution.
+
+    For a linear EDC e(t) = 1 - t/T over [0, T]:
+
+        T_s = integral(e(t), 0, T) / e(0) = (T/2) / 1 = T/2
+    """
+    T = 0.1  # total duration in seconds
+    # non-uniform time grid: dense at start, coarse at end
+    times = np.concatenate([
+        np.linspace(0, 0.02, 20, endpoint=False),
+        np.linspace(0.02, T, 10),
+    ])
+    edc_values = 1 - times / T
+    edc = pf.TimeData(edc_values[np.newaxis, :], times)
+    result = center_time(edc)
+    npt.assert_allclose(result, T / 2, rtol=1e-6)
+
+def test_center_time_exponential_decay_analytical(make_edc):
+    r"""Center time for exponential EDC matches analytical solution.
+
+    For e(t) = exp(-alpha * t) with alpha = 13.8155 / RT60:
+
+        T_s = integral(e(t), 0, inf) / e(0) = 1 / alpha
+    """
+    rt60 = 2.0
+    sampling_rate = 1000
+    total_samples = 5000
+    edc = make_edc(rt=rt60, sampling_rate=sampling_rate,
+                   total_samples=total_samples)
+    result = center_time(edc)
+
+    # Analytical expected value
+    a = 13.8155 / rt60
+    expected = 1 / a
+    npt.assert_allclose(result, expected, atol=1e-3)
+
+def test_center_time_multichannel(make_edc):
+    """Each channel is computed independently and results differ."""
+    energy = np.stack([
+        np.exp(-13.8155 / 1.0 * np.arange(2000) / 1000),
+        np.exp(-13.8155 / 2.0 * np.arange(2000) / 1000),
+    ])
+    edc = make_edc(energy=energy, sampling_rate=1000)
+    result = center_time(edc)
+    assert result[1] > result[0]
+
+def test_center_time_nan_tail_returns_finite_result(make_edc):
+    """center_time() returns a finite result when the EDC tail is NaN.
+
+    Lundeby/Chu methods set the noise tail to NaN. The finite head should
+    still produce a valid Ts.
+    """
+    rt60 = 1.0
+    sampling_rate = 1000
+    total_samples = 2000
+    edc = make_edc(rt=rt60, sampling_rate=sampling_rate,
+                   total_samples=total_samples)
+    # Simulate Lundeby/Chu truncation: set the last 500 samples to NaN
+    edc.time[..., 1500:] = np.nan
+
+    result = center_time(edc)
+
+    assert np.all(np.isfinite(result))