Refactor the parametric energy decay curve function

pyrato/parametric.py

@@ -4,6 +4,9 @@
 such as Sabine's theory of sound in rooms.
 """
 import numpy as np
+from numpy.typing import NDArray
+import pyfar as pf
+
 
 def energy_decay_curve_analytic(
         surfaces, alphas, volume, times, source=None,
@@ -81,6 +84,64 @@ def energy_decay_curve_analytic(
     return energy_decay_curve
 
 
+def energy_decay_curve(
+        times : np.ndarray[float],
+        reverberation_time : float | np.ndarray[float],
+        energy : float | np.ndarray[float] = 1,
+    ) -> pf.TimeData:
+    r"""Calculate the energy decay curve from the reverberation time and energy.
+
+    The energy decay curve is calculated as
+
+    .. math::
+        E(t) = E_0 e^{-\frac{6 \ln(10)}{T_{60}} t}
+
+    where :math:`E_0` is the initial energy, :math:`T_{60}` the reverberation
+    time, and :math:`t` the time [#]_.
+
+    Parameters
+    ----------
+    times : numpy.ndarray[float]
+        The times at which the energy decay curve is evaluated.
+    reverberation_time : float | numpy.ndarray[float]
+        The reverberation time in seconds.
+    energy : float | numpy.ndarray[float], optional
+        The initial energy of the sound field, by default 1
+
+    Returns
+    -------
+    pyfar.TimeData
+        The energy decay curve.
+
+    Example
+    -------
+    Calculate and plot an energy decay curve with a reverberation time of
+    2 seconds.
+
+    .. plot::
+
+        >>> import numpy as np
+        >>> import pyrato
+        >>> import pyfar as pf
+        >>>
+        >>> times = np.linspace(0, 3, 50)
+        >>> T_60 = 2
+        >>> edc = pyrato.parametric.energy_decay_curve(times, T_60)
+        >>> pf.plot.time(edc, log_prefix=10, dB=True)
+
+
+    References
+    ----------
+    .. [#] H. Kuttruff, Room acoustics, 4th Ed. Taylor & Francis, 2009.
+
+    """
+
+    damping_term = 3*np.log(10) / reverberation_time
+    edc = energy * np.exp(-2*damping_term*times)
+
+    return pf.TimeData(edc, times)
+
+
 def air_attenuation_coefficient(
         frequency,
         temperature=20,

tests/test_edc.py

@@ -62,3 +62,36 @@ def test_edc_sabine():
         6.37107964e-04, 5.46555336e-04,
         ])
     npt.assert_almost_equal(edc, truth)
+
+
+def test_edc_function():
+    times = np.linspace(0, 0.25, 50)
+
+    c = 343.4
+
+    volume = 2*2*2
+    alphas = np.asarray([0.9, 0.1])
+    surfaces = np.asarray([2, 5*2])
+    surface_room = np.sum(surfaces)
+    alpha_mean = np.sum(surfaces*alphas) / surface_room
+
+    rt = 24*np.log(10)/c * volume / (surface_room*alpha_mean)
+    energy = 1
+    edc = ra.parametric.energy_decay_curve(times, rt, energy)
+
+    truth = array([
+        1.00000000e+00, 8.57869258e-01, 7.35939663e-01, 6.31340013e-01,
+        5.41607188e-01, 4.64628156e-01, 3.98590212e-01, 3.41938289e-01,
+        2.93338346e-01, 2.51645949e-01, 2.15879324e-01, 1.85196235e-01,
+        1.58874157e-01, 1.36293255e-01, 1.16921793e-01, 1.00303612e-01,
+        8.60473853e-02, 7.38174065e-02, 6.33256837e-02, 5.43251573e-02,
+        4.66038824e-02, 3.99800380e-02, 3.42976455e-02, 2.94228957e-02,
+        2.52409977e-02, 2.16534759e-02, 1.85758513e-02, 1.59356518e-02,
+        1.36707058e-02, 1.17276782e-02, 1.00608146e-02, 8.63086356e-03,
+        7.40415251e-03, 6.35179482e-03, 5.44900951e-03, 4.67453774e-03,
+        4.01014222e-03, 3.44017773e-03, 2.95122272e-03, 2.53176324e-03,
+        2.17192185e-03, 1.86322499e-03, 1.59840344e-03, 1.37122117e-03,
+        1.17632849e-03, 1.00913605e-03, 8.65706791e-04, 7.42663242e-04,
+        6.37107964e-04, 5.46555336e-04,
+        ])
+    npt.assert_almost_equal(edc.time, np.atleast_2d(truth))