diff --git a/pyrato/parameters.py b/pyrato/parameters.py index 1430a800..1c1f992e 100644 --- a/pyrato/parameters.py +++ b/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 diff --git a/tests/test_parameters.py b/tests/test_parameters.py index b9c38bfb..6c71acd6 100644 --- a/tests/test_parameters.py +++ b/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)