robustifies the RT60 computation function (#347)

Author: fakufaku

CHANGELOG.rst

@@ -11,7 +11,12 @@ adheres to `Semantic Versioning <http://semver.org/spec/v2.0.0.html>`_.
 `Unreleased`_
 -------------
 
-Nothing yet
+Changed
+~~~~~~~
+
+- Extra parameter ``energy_thresh`` added to ``pyroomacoustics.experimental.measure_rt60``.
+  The energy tail beyond this threshold is discarded which is useful for noisy RIR
+  measurements. The default value is 0.95.
 
 `0.7.4`_ - 2024-04-25
 ---------------------

pyroomacoustics/bss/fastmnmf.py

@@ -161,9 +161,7 @@ def separate():
         Y_FTM = np.einsum("nft, nfm -> ftm", lambda_NFT, G_NFM)
 
         # update G_NFM (diagonal element of spatial covariance matrices)
-        numerator = np.einsum(
-            "nft, ftm -> nfm", lambda_NFT, Qx_power_FTM / (Y_FTM**2)
-        )
+        numerator = np.einsum("nft, ftm -> nfm", lambda_NFT, Qx_power_FTM / (Y_FTM**2))
         denominator = np.einsum("nft, ftm -> nfm", lambda_NFT, 1 / Y_FTM) + eps
         G_NFM *= np.sqrt(numerator / denominator)
         Y_FTM = np.einsum("nft, nfm -> ftm", lambda_NFT, G_NFM)

pyroomacoustics/experimental/rt60.py

@@ -33,7 +33,7 @@
 import numpy as np
 
 
-def measure_rt60(h, fs=1, decay_db=60, plot=False, rt60_tgt=None):
+def measure_rt60(h, fs=1, decay_db=60, energy_thres=0.95, plot=False, rt60_tgt=None):
     """
     Analyze the RT60 of an impulse response. Optionaly plots some useful information.
 
@@ -47,6 +47,10 @@ def measure_rt60(h, fs=1, decay_db=60, plot=False, rt60_tgt=None):
         The decay in decibels for which we actually estimate the time. Although
         we would like to estimate the RT60, it might not be practical. Instead,
         we measure the RT20 or RT30 and extrapolate to RT60.
+    energy_thres: float
+        This should be a value between 0.0 and 1.0.
+        If provided, the fit will be done using a fraction energy_thres of the
+        whole energy. This is useful when there is a long noisy tail for example.
     plot: bool, optional
         If set to ``True``, the power decay and different estimated values will
         be plotted (default False).
@@ -60,21 +64,36 @@ def measure_rt60(h, fs=1, decay_db=60, plot=False, rt60_tgt=None):
 
     # The power of the impulse response in dB
     power = h**2
+    # Backward energy integration according to Schroeder
     energy = np.cumsum(power[::-1])[::-1]  # Integration according to Schroeder
 
+    if energy_thres < 1.0:
+        assert 0.0 < energy_thres < 1.0
+        energy -= energy[0] * (1.0 - energy_thres)
+        energy = np.maximum(energy, 0.0)
+
     # remove the possibly all zero tail
     i_nz = np.max(np.where(energy > 0)[0])
     energy = energy[:i_nz]
     energy_db = 10 * np.log10(energy)
     energy_db -= energy_db[0]
 
+    min_energy_db = -np.min(energy_db)
+    if min_energy_db - 5 < decay_db:
+        decay_db = min_energy_db
+
     # -5 dB headroom
-    i_5db = np.min(np.where(-5 - energy_db > 0)[0])
+    try:
+        i_5db = np.min(np.where(energy_db < -5)[0])
+    except ValueError:
+        return 0.0
     e_5db = energy_db[i_5db]
     t_5db = i_5db / fs
-
     # after decay
-    i_decay = np.min(np.where(-5 - decay_db - energy_db > 0)[0])
+    try:
+        i_decay = np.min(np.where(energy_db < -5 - decay_db)[0])
+    except ValueError:
+        i_decay = len(energy_db)
     t_decay = i_decay / fs
 
     # compute the decay time