Feature/mixing time

pyrato/parameters.py

@@ -5,6 +5,7 @@
 import re
 import numpy as np
 import pyfar as pf
+from math import erfc
 import warnings
 
 def reverberation_time_linear_regression(
@@ -662,6 +663,170 @@ def sound_strength(energy_decay_curve_room,
                                      energy_decay_curve_room))
 
 
+def mixing_time(rir, window_length=None, lindau_regression=False,
+                onset_threshold_db=-30.0, peak_secure_margin=100):
+    r"""Estimate the perceptual mixing time using the echo density.
+
+    The mixing time is estimated from the echo density profile according to
+    criterion I of Abel and Huang [#abel]_, optionally followed by the
+    perceptual regression of Lindau et al. [#lindau]_.
+
+    The direct-sound onset is detected per channel as the first sample
+    exceeding ``onset_threshold_db`` below the peak, and the impulse response
+    is cut from ``onset - peak_secure_margin``, the margin dropping to zero if
+    the onset falls inside it. Multichannel input is cut at the earliest onset
+    across channels, which preserves inter-channel time differences: a channel
+    is then measured from the first arrival across all channels, so its
+    transition time includes its propagation delay relative to the earliest
+    channel. Pass channels individually to measure each from its own onset.
+
+    The echo density counts the samples with :math:`|h| > \sigma` inside a
+    rectangular window (growing at the start, of length ``N`` in the middle,
+    shrinking at the end), divided by ``N`` in every case and by the Gaussian
+    expectation :math:`\mathrm{erfc}(1/\sqrt{2}) \approx 0.317`. Here
+    :math:`\sigma` is the sample standard deviation of the window. The
+    transition time is
+
+    .. math::
+
+        t_\mathrm{Abel} = (d - \mathrm{margin}) / f_s,
+
+    where :math:`d` is the first sample at which the echo density exceeds
+    unity, i.e. the field is statistically diffuse. If ``lindau_regression``
+    is ``True``, the perceptual mixing time follows from the regression
+
+    .. math::
+
+        t_\mathrm{mp50} = 0.8 \cdot t_\mathrm{Abel} - 8\,\mathrm{ms},
+
+    with negative results clipped to 1 ms.
+
+    Parameters
+    ----------
+    rir : pyfar.Signal
+        Room impulse response.
+    window_length : float, optional
+        Sliding window length in seconds, rounded down to an even number of
+        samples and at least 4. The default ``None`` uses 1024 samples (about
+        21 ms at 48 kHz), matching the window recommended by Abel and Huang
+        [#abel]_.
+    lindau_regression : bool, optional
+        If ``True``, apply the Lindau et al. [#lindau]_ ``tmp50`` regression
+        to obtain the perceptual mixing time. If ``False``, return the raw
+        transition time :math:`t_\mathrm{Abel}`. The default is ``False``.
+    onset_threshold_db : float, optional
+        Peak criterion for onset detection, in dB below the peak. The default
+        is ``-30.0``.
+    peak_secure_margin : int, optional
+        Safety margin in samples kept before the detected onset, and
+        subtracted from the transition time so it is onset-relative. The
+        default is ``100``.
+
+    Returns
+    -------
+    mixing_time : numpy.ndarray
+        Estimated mixing time in seconds, of shape ``rir.cshape`` with one
+        value per channel.
+
+    Raises
+    ------
+    TypeError
+        If ``rir`` is not a :py:class:`pyfar.Signal`.
+    ValueError
+        If a channel is silent, if the impulse response is shorter than the
+        analysis window, or if the echo density never exceeds unity.
+
+    References
+    ----------
+    .. [#abel] J. S. Abel and P. Huang, "A simple, robust measure of
+        reverberation echo density," in Proc. of the 121st AES Convention,
+        San Francisco, 2006.
+    .. [#lindau] A. Lindau, L. Kosanke, and S. Weinzierl, "Perceptual
+        evaluation of model- and signal-based predictors of the mixing time
+        in binaural room impulse responses," J. Audio Eng. Soc., vol. 60,
+        no. 11, pp. 887-898, 2012.
+
+    Examples
+    --------
+    Estimate the perceptual mixing time of a measured room impulse response:
+
+    >>> import pyfar as pf
+    >>> import pyrato as ra
+    ...
+    >>> rir = pf.signals.files.room_impulse_response(sampling_rate=48000)
+    >>> t_mix = ra.parameters.mixing_time(rir, lindau_regression=True)
+    """
+    if not isinstance(rir, pf.Signal):
+        raise TypeError("Input data must be a pyfar.Signal.")
+
+    fs = rir.sampling_rate
+    cshape = rir.cshape
+    data = rir.time.reshape(-1, rir.n_samples)
+    n_channels = data.shape[0]
+
+    if window_length is None:
+        N = 1024
+    else:
+        # Must be even: the window holds 2 * (N // 2) samples, normalised by N.
+        N = round(window_length * float(fs))
+        N = max(4, N - N % 2)
+    half = N // 2
+
+    peaks = np.max(np.abs(data), axis=-1)
+    if np.any(peaks <= 0.0):
+        raise ValueError("cannot detect an onset: a channel is silent")
+    threshold = 10.0 ** (onset_threshold_db / 20.0)
+    # No onset guard needed: the peak sample always clears the criterion.
+    above = np.abs(data) > peaks[:, np.newaxis] * threshold
+    min_onset = int(np.argmax(above, axis=-1).min())
+    # Margin is all-or-nothing, and decoupled from the one subtracted from
+    # t_abel below, which is always peak_secure_margin.
+    cut_margin = (
+        0 if (min_onset + 1) <= peak_secure_margin else peak_secure_margin)
+    data = data[:, min_onset - cut_margin:]
+    L = data.shape[-1]
+    if L < N:
+        raise ValueError(
+            f"IR shorter than analysis window length ({N} samples). "
+            "Provide at least an IR of some 100 msec.")
+
+    # Normalises the outlier fraction so Gaussian noise gives echo density 1.
+    p_gauss_norm = 1.0 / erfc(1.0 / np.sqrt(2.0))
+
+    out = np.empty(n_channels)
+    for c in range(n_channels):
+        x = data[c]
+        d = None
+        for k in range(1, L + 1):                  # 1-based window centre
+            if k <= half + 1:                      # growing window
+                window = x[0 : k + half - 1]
+            elif k <= L - half + 1:                # constant window (length N)
+                window = x[k - half - 1 : k + half - 1]
+            else:                                  # shrinking window
+                window = x[k - half - 1 : L]
+            sigma = np.std(window, ddof=1)
+            n_out = np.count_nonzero(np.abs(window) > sigma)
+            echo_dens = n_out / N * p_gauss_norm
+            if echo_dens > 1.0:
+                d = k
+                break
+
+        if d is None:
+            raise ValueError(
+                "Mixing time not found within given temporal limits. "
+                "Try again with extended stopping criterion.")
+
+        t_abel_I = (d - peak_secure_margin) / float(fs)   # onset-relative
+
+        if lindau_regression:
+            t_abel_I = 0.8 * t_abel_I - 0.008
+            if t_abel_I < 0:
+                t_abel_I = 0.001
+        out[c] = t_abel_I
+
+    return out.reshape(cshape)
+
+
 def speech_transmission_index_indirect(
         rir, rir_type="acoustical", level=None, snr=np.inf,
         ambient_noise_correction=True):

tests/test_parameters_mixing_time.py

@@ -0,0 +1,185 @@
+import numpy as np
+import pytest
+import pyfar as pf
+import numpy.testing as npt
+
+from pyrato.parameters import mixing_time
+
+
+def _noise_rir(n_samples=20000, onset=500, seed=0, rt=0.25):
+    """Exponentially decaying Gaussian noise with a direct sound at ``onset``.
+
+    The samples before ``onset`` are zero, so the onset detection returns
+    ``onset`` rather than an early noise sample. The echo density of Gaussian
+    noise crosses unity shortly after the onset, so the transition time is
+    well defined.
+    """
+    rng = np.random.default_rng(seed)
+    n_tail = n_samples - onset
+    decay = np.exp(-np.arange(n_tail) / n_tail / rt)
+    rir = np.zeros(n_samples)
+    rir[onset:] = rng.standard_normal(n_tail) * decay
+    rir[onset] += 10.0
+    return rir
+
+
+@pytest.mark.parametrize("cshape", [(1,), (2,), (2, 3)])
+def test_mixing_time_returns_correct_shape(cshape):
+    """Return one value per channel, matching the input cshape."""
+    n_channels = int(np.prod(cshape))
+    data = np.stack([_noise_rir(seed=i) for i in range(n_channels)])
+    rir = pf.Signal(data.reshape(*cshape, -1), 48000)
+
+    result = mixing_time(rir)
+
+    assert isinstance(result, np.ndarray)
+    assert result.shape == rir.cshape
+
+
+def test_mixing_time_rejects_non_signal_input():
+    """TypeError is raised when input data is not a pyfar.Signal."""
+    match = "Input data must be a pyfar.Signal."
+    with pytest.raises(TypeError, match=match):
+        mixing_time(_noise_rir())
+
+
+def test_mixing_time_rejects_silent_channel():
+    """ValueError is raised when a channel carries no signal."""
+    rir = pf.Signal(np.zeros(20000), 48000)
+    match = "cannot detect an onset: a channel is silent"
+    with pytest.raises(ValueError, match=match):
+        mixing_time(rir)
+
+
+def test_mixing_time_rejects_ir_shorter_than_window():
+    """ValueError is raised when the IR is shorter than the analysis window."""
+    rir = pf.Signal(_noise_rir(n_samples=512, onset=10), 48000)
+    match = "IR shorter than analysis window length"
+    with pytest.raises(ValueError, match=match):
+        mixing_time(rir)
+
+
+def test_mixing_time_rejects_ir_without_diffuse_field():
+    """ValueError is raised when the echo density never exceeds unity.
+
+    A Dirac in silence is maximally sparse: inside every window the sample
+    standard deviation is set by the single non-zero sample, so the outlier
+    fraction never reaches the Gaussian expectation.
+    """
+    rir = pf.signals.impulse(5000)
+    match = "Mixing time not found within given temporal limits"
+    with pytest.raises(ValueError, match=match):
+        mixing_time(rir)
+
+
+def test_mixing_time_matches_reference_values():
+    """Transition time and tmp50 against the reference implementation.
+
+    Cross-checked against AKcutIRmixingTime + AKmixingTimeAbel +
+    AKdataBasedMixingTime (AKtools, TU Berlin) on the same impulse response.
+    Both implementations agree, so the tolerance is tight.
+    """
+    time = np.loadtxt('./tests/test_data/room_impulse_response_with_noise.csv')
+    rir = pf.Signal(time, 48000)
+
+    npt.assert_allclose(mixing_time(rir), 72.8333333333e-3, atol=1e-9)
+    npt.assert_allclose(
+        mixing_time(rir, lindau_regression=True), 50.2666666667e-3, atol=1e-9)
+
+
+def test_mixing_time_lindau_regression_maps_transition_time():
+    """tmp50 is an affine map of the transition time: 0.8 * t - 8 ms."""
+    rir = pf.Signal(_noise_rir(), 48000)
+
+    t_abel = mixing_time(rir)
+    t_mp50 = mixing_time(rir, lindau_regression=True)
+
+    assert np.all(0.8 * t_abel - 0.008 > 0)  # clipping branch not exercised
+    npt.assert_allclose(t_mp50, 0.8 * t_abel - 0.008, atol=1e-12)
+
+
+def test_mixing_time_lindau_regression_clips_negative_to_one_ms():
+    """A negative tmp50 is clipped to 1 ms."""
+    # The transition is found at the same sample regardless of the sampling
+    # rate, so a high rate drives the transition time low enough for the
+    # regression to go negative.
+    rir = pf.Signal(_noise_rir(onset=20, seed=3), 192000)
+
+    assert 0.8 * mixing_time(rir)[0] - 0.008 < 0
+    npt.assert_allclose(mixing_time(rir, lindau_regression=True), 0.001)
+
+
+def test_mixing_time_rounds_window_down_to_even_samples():
+    """An odd window length is rounded down to an even sample count.
+
+    The analysis window spans ``2 * (N // 2)`` samples while the echo density
+    normalises by ``N``, so an odd ``N`` would bias the density low.
+
+    The measured impulse response is used rather than the synthetic one: the
+    bias is small, and only a signal whose echo density crosses unity close to
+    a sample boundary resolves 1024 from an unrounded 1025.
+    """
+    time = np.loadtxt('./tests/test_data/room_impulse_response_with_noise.csv')
+    rir = pf.Signal(time, 48000)
+
+    npt.assert_allclose(
+        mixing_time(rir, window_length=1025 / 48000),
+        mixing_time(rir, window_length=1024 / 48000))
+
+
+def test_mixing_time_window_length_changes_result():
+    """The window length is not silently ignored."""
+    rir = pf.Signal(_noise_rir(), 48000)
+
+    assert not np.isclose(
+        mixing_time(rir, window_length=512 / 48000), mixing_time(rir))
+
+
+def test_mixing_time_window_has_a_lower_bound_of_four_samples():
+    """Window lengths below 4 samples are clipped to 4."""
+    rir = pf.Signal(_noise_rir(), 48000)
+
+    npt.assert_allclose(
+        mixing_time(rir, window_length=1 / 48000),
+        mixing_time(rir, window_length=4 / 48000))
+
+
+def test_mixing_time_cuts_all_channels_at_earliest_onset():
+    """Multichannel input is cut at the earliest onset across channels.
+
+    The channel holding the earliest onset is unaffected by the presence of
+    the other, while a later-onset channel is not aligned to its own direct
+    sound. The delayed channel comes first so that the earliest onset is not
+    simply the onset of channel zero.
+
+    The 4.83 ms of the delayed channel is not a meaningful mixing time: the
+    shared cut places its analysis window in the pre-onset noise floor of the
+    measured impulse response, whose Gaussian statistics cross unity at once.
+    It pins the cut behaviour, and matches AKtools on the same input.
+    """
+    time = np.loadtxt('./tests/test_data/room_impulse_response_with_noise.csv')
+    delayed = np.concatenate([np.zeros(1500), time[:-1500]])
+
+    both = mixing_time(pf.Signal(np.stack([delayed, time]), 48000))
+    delayed_alone = mixing_time(pf.Signal(delayed, 48000))
+    time_alone = mixing_time(pf.Signal(time, 48000))
+
+    # the earliest onset belongs to the second channel, which is unaffected
+    npt.assert_allclose(both[1], time_alone[0])
+    assert not np.isclose(both[0], delayed_alone[0])
+
+    npt.assert_allclose(both, [4.8333333333e-3, 72.8333333333e-3], atol=1e-9)
+
+
+def test_mixing_time_keeps_full_margin_when_onset_precedes_it():
+    """An onset inside the first ``peak_secure_margin`` samples drops the cut
+    margin to zero, while the full margin is still subtracted from the
+    transition time.
+
+    This decoupling is inherited from the reference implementation, where the
+    cut function zeroes only its local copy of the margin. Cross-checked
+    against AKtools, which returns the same value.
+    """
+    rir = pf.Signal(_noise_rir(onset=20, seed=3), 48000)
+
+    npt.assert_allclose(mixing_time(rir), 11.4583333333e-3, atol=1e-9)