Fix octave band filters and their bandwidths

CHANGELOG.rst

@@ -10,13 +10,29 @@ adheres to `Semantic Versioning <http://semver.org/spec/v2.0.0.html>`_.
 
 `Unreleased`_
 -------------
-Bugfix:
+
+Changed
+~~~~~~~
+
+- Refactor the way the RIR is weighted with the histogram in simulation/rt.py.
+
+- Improves ``pra.experimental.measure_rt60``: Compute the T60 using a fit. Default is
+  log-domain. Adds option to fit in linear domain.
+
+Bugfix
 ~~~~~~
 
-- In `doa.py`, the `ax.xaxis.grid` and `ax.yaxis.grid` parameters were changed from `b` to `visible`.
+- Fixes the lowest octave band filter that was malformed when using
+  ``octave_bands_keep_dc=True``.
+
+- Fixes the computation of the octave band widths that were not correct for the lowest
+  and high bands.
+
+- In ``doa.py``, the ``ax.xaxis.grid`` and ``ax.yaxis.grid`` parameters were changed from ``b`` to ``visible``.
+
+- Fixes ``MicrophoneArray.append(MicrophoneArray) AttributeError: 'MicrophoneArray'
+  object has no attribute 'shape'.``
 
-- Fixes MicrophoneArray.append(MicrophoneArray)
-  AttributeError: 'MicrophoneArray' object has no attribute 'shape'.
 - Fixes issue #421: When generating a highpass filtered room impulse response make
   sure that the output is a memory contiguous NumPy array.
 

examples/room_shoebox_3d_rt.py

@@ -104,7 +104,7 @@ def chrono(f, *args, **kwargs):
         plt.legend()
 
         plt.figure(2)
-        rt60 = shoebox.measure_rt60(plot=True, decay_db=60)
+        rt60 = shoebox.measure_rt60(plot=True, decay_db=30.0, linear_domain_fit=True)
 
         print(f"  RT60:")
         print(f"    - Eyring   {rt60_eyring:.3f} s")

pyroomacoustics/acoustics.py

@@ -180,7 +180,12 @@ def __init__(self, base_frequency=125.0, fs=16000, n_fft=512, keep_dc=False):
 
     def get_bw(self):
         """Returns the bandwidth of the bands"""
-        return np.array([b2 - b1 for b1, b2 in self.bands])
+        bandwidths = []
+        for b, (b1, b2) in enumerate(self.bands):
+            lo = 0.0 if b == 0 and self.keep_dc else b1
+            hi = min(b2, self.fs / 2.0)
+            bandwidths.append(hi - lo)
+        return np.array(bandwidths)
 
     def analysis(self, x, band=None):
         """
@@ -290,14 +295,17 @@ def _make_filters(self):
                 make_one = freq < center
                 freq_resp[make_one, b] = 1.0
 
-            lo = np.logical_and(band[0] <= freq, freq < center)
-
-            freq_resp[lo, b] = 0.5 * (1 + np.cos(2 * np.pi * freq[lo] / center))
+            else:
+                # Rising side.
+                lo = np.logical_and(band[0] <= freq, freq < center)
+                freq_resp[lo, b] = 0.5 * (1 + np.cos(2 * np.pi * freq[lo] / center))
 
             if b != n - 1:
+                # Falling side.
                 hi = np.logical_and(center <= freq, freq < band[1])
                 freq_resp[hi, b] = 0.5 * (1 - np.cos(2 * np.pi * freq[hi] / band[1]))
             else:
+                # Shelve for the last octave band.
                 hi = center <= freq
                 freq_resp[hi, b] = 1.0
 
@@ -606,7 +614,7 @@ def rt60_eyring(S, V, a, m, c):
 
 def rt60_sabine(S, V, a, m, c):
     """
-    This is the Eyring formula for estimation of the reverberation time.
+    This is the Sabine formula for estimation of the reverberation time.
 
     Parameters
     ----------

pyroomacoustics/experimental/rt60.py

@@ -30,10 +30,59 @@
 .. [1] M. R. Schroeder, "New Method of Measuring Reverberation Time,"
     J. Acoust. Soc. Am., vol. 37, no. 3, pp. 409-412, Mar. 1968.
 """
+import math
+
 import numpy as np
+from scipy.optimize import curve_fit
+
+
+def _fit_exp_and_extrapolate(
+    data_db, fs, extrapolate_value_db=-60.0, linear_domain_fit=False
+):
+    """Non-linear fit to an exponential f(t) = b * np.exp(a * t)."""
+    # Fix the origin.
+    data = data_db - data_db[0]
+
+    N = data.shape[0]
+    t = np.arange(N) / fs
+
+    # We use a least-square fit in log-domain as initialization.
+    X = np.column_stack((t, np.ones(N)))
+    p, *_ = np.linalg.lstsq(X, data)
+
+    if not linear_domain_fit:
+        return extrapolate_value_db / p[0]
+
+    # Non-linear fit using scipy.optimize.curve_fit.
 
+    def model(x, a, b):
+        return b * np.exp(a * x)
 
-def measure_rt60(h, fs=1, decay_db=60, energy_thres=1.0, plot=False, rt60_tgt=None):
+    def jac(x, a, b):
+        u = np.exp(a * x)
+        return np.column_stack((x * b * u, u))
+
+    # Provide an initial guess [a, b]
+    initial_guess = [p[0] / (10.0 * np.log10(np.e)), 10.0 ** (p[1] / 10.0)]
+
+    # Perform the fit in the linear domain.
+    # Empirically, this takes less than 10 iterations.
+    popt, pcov = curve_fit(model, t, 10.0 ** (data / 10.0), p0=initial_guess, jac=jac)
+
+    # This is the estimate of the T60 (or other value) based on the fit.
+    return (extrapolate_value_db / 10.0) * np.log(10.0) / popt[0]
+
+
+def measure_rt60(
+    h,
+    fs=1,
+    decay_db=60,
+    energy_thres=1.0,
+    plot=False,
+    rt60_tgt=None,
+    label=None,
+    linear_domain_fit=False,
+):
     """
     Analyze the RT60 of an impulse response. Optionaly plots some useful information.
 
@@ -57,6 +106,11 @@ def measure_rt60(h, fs=1, decay_db=60, energy_thres=1.0, plot=False, rt60_tgt=No
     rt60_tgt: float
         This parameter can be used to indicate a target RT60 to which we want
         to compare the estimated value.
+    label: str, optional
+        A label to use in a plot.
+    linear_domain_fit: bool, optional
+        When True, applies a direct fit of an exponential to the data in the linear domain.
+        When False, a linear fit is done in the logarithm domain (default).
     """
 
     h = np.array(h)
@@ -84,21 +138,24 @@ def measure_rt60(h, fs=1, decay_db=60, energy_thres=1.0, plot=False, rt60_tgt=No
 
     # -5 dB headroom
     try:
-        i_5db = np.min(np.where(energy_db < -5)[0])
+        i_5db = np.min(np.where(energy_db < -5.0)[0])
     except ValueError:
         return 0.0
     e_5db = energy_db[i_5db]
-    t_5db = i_5db / fs
+
     # after decay
     try:
-        i_decay = np.min(np.where(energy_db < -5 - decay_db)[0])
+        i_decay = np.min(np.where(energy_db < e_5db - decay_db)[0])
     except ValueError:
         i_decay = len(energy_db)
-    t_decay = i_decay / fs
 
-    # compute the decay time
-    decay_time = t_decay - t_5db
-    est_rt60 = (60 / decay_db) * decay_time
+    # Compute the RT60 estimate using a fitting method.
+    est_rt60 = _fit_exp_and_extrapolate(
+        energy_db[i_5db:i_decay],
+        fs,
+        extrapolate_value_db=-60.0,
+        linear_domain_fit=linear_domain_fit,
+    )
 
     if plot:
         import matplotlib.pyplot as plt
@@ -116,18 +173,26 @@ def get_time(x, fs):
 
         T = get_time(power_db, fs)
 
+        # Optional label for the legend.
+        label = f" {label}" if label else ""
+
         # plot power and energy
-        plt.plot(get_time(energy_db, fs), energy_db, label="Energy")
+        plt.plot(get_time(energy_db, fs), energy_db, label=f"Energy{label}")
 
         # now the linear fit
-        plt.plot([0, est_rt60], [e_5db, -65], "--", label="Linear Fit")
-        plt.plot(T, np.ones_like(T) * -60, "--", label="-60 dB")
+        plt.plot([0, est_rt60], [e_5db, -65], "--", label=f"Linear Fit{label}")
+        plt.plot(T, np.ones_like(T) * -60, "--", label=f"-60 dB{label}")
+        plt.plot(T, np.ones_like(T) * -5.0, "--", label=f"-5 dB{label}")
         plt.vlines(
-            est_rt60, energy_db_min, 0, linestyles="dashed", label="Estimated RT60"
+            est_rt60,
+            energy_db_min,
+            0,
+            linestyles="dashed",
+            label=f"Estimated RT60{label}",
         )
 
         if rt60_tgt is not None:
-            plt.vlines(rt60_tgt, energy_db_min, 0, label="Target RT60")
+            plt.vlines(rt60_tgt, energy_db_min, 0, label=f"Target RT60{label}")
 
         plt.legend()
 

pyroomacoustics/libroom_src/common.hpp

@@ -169,6 +169,11 @@ class Histogram2D
     {
       return array;
     }
+
+    Eigen::ArrayXXi get_counts() const
+    {
+      return counts;
+    }
 };
 
 #endif // __COMMON_HPP__

pyroomacoustics/libroom_src/libroom.cpp

@@ -242,6 +242,7 @@ PYBIND11_MODULE(libroom, m) {
       .def("log", &Histogram2D::log)
       .def("bin", &Histogram2D::bin)
       .def("get_hist", &Histogram2D::get_hist)
+      .def("get_counts", &Histogram2D::get_counts)
       .def("reset", &Histogram2D::reset);
 
   // Structure to hold detector hit information

pyroomacoustics/libroom_src/room.cpp

@@ -32,6 +32,10 @@
 const double pi = 3.14159265358979323846;
 const double pi_2 = 1.57079632679489661923;
 
+// Initial energy of a particule.
+// The value 2.0 is necessary to match the scale of the ISM.
+const double energy_0_numerator = 2.0f;
+
 size_t number_image_sources_2(size_t max_order) {
   /*
   ¦* The number of image sources for a given maximum order in 2D
@@ -995,7 +999,7 @@ void Room<D>::ray_tracing(
   )
 {
   // float energy_0 = 2.f / (mic_radius * mic_radius * angles.cols());
-  float energy_0 = 2.f / angles.cols();
+  float energy_0 = energy_0_numerator / angles.cols();
 
   for (int k(0) ; k < angles.cols() ; k++)
   {
@@ -1035,7 +1039,7 @@ void Room<D>::ray_tracing(
 
   // ------------------ INIT --------------------
   // initial energy of one ray
-  float energy_0 = 2.f / (nb_phis * nb_thetas);
+  float energy_0 = energy_0_numerator / (nb_phis * nb_thetas);
 
   // ------------------ RAY TRACING --------------------
 
@@ -1085,7 +1089,7 @@ void Room<D>::ray_tracing(
 
   // ------------------ INIT --------------------
   // initial energy of one ray
-  float energy_0 = 2.f / n_rays;
+  float energy_0 = energy_0_numerator / n_rays;
 
   // ------------------ RAY TRACING --------------------
   if (D == 3)

pyroomacoustics/room.py

@@ -2272,7 +2272,8 @@ def ray_tracing(self):
             # reset all the receivers' histograms
             self.room_engine.reset_mics()
 
-        # Basically, histograms for 2 mics corresponding to each source , the histograms are in each octave bands hence (7,2500) 2500 histogram length
+        # Basically, histograms for 2 mics corresponding to each source , the
+        # histograms are in each octave bands hence (7,2500) 2500 histogram length
         # update the state
         self.simulator_state["rt_done"] = True
 
@@ -2931,46 +2932,34 @@ def rt60_theory(self, formula="sabine"):
         formula: str
             The formula to use for the calculation, 'sabine' (default) or 'eyring'
         """
-
-        rt60 = 0.0
-
         if self.is_multi_band:
             bandwidths = self.octave_bands.get_bw()
         else:
             bandwidths = [1.0]
 
         V = self.volume
-        S = np.sum([w.area() for w in self.walls])
         c = self.c
 
+        if formula == "eyring":
+            rt60_fn = rt60_eyring
+        elif formula == "sabine":
+            rt60_fn = rt60_sabine
+        else:
+            raise ValueError("Only Eyring and Sabine's formulas are supported")
+
+        # The harmonic mean will naturally average the absorption over the Octave bands
+        # and walls.
+        inv_rt60 = 0.0
         for i, bw in enumerate(bandwidths):
-            # average absorption coefficients
-            a = 0.0
+            m = 0.0 if self.air_absorption is None else self.air_absorption[i]
             for w in self.walls:
-                if len(w.absorption) == 1:
-                    a += w.area() * w.absorption[0]
-                else:
-                    a += w.area() * w.absorption[i]
-            a /= S
-
-            try:
-                m = self.air_absorption[i]
-            except:
-                m = 0.0
-
-            if formula == "eyring":
-                rt60_loc = rt60_eyring(S, V, a, m, c)
-            elif formula == "sabine":
-                rt60_loc = rt60_sabine(S, V, a, m, c)
-            else:
-                raise ValueError("Only Eyring and Sabine's formulas are supported")
-
-            rt60 += rt60_loc * bw
+                inv_rt60 += bw / rt60_fn(w.area(), V, w.absorption[i], m, c)
 
-        rt60 /= np.sum(bandwidths)
-        return rt60
+        return np.sum(bandwidths) / inv_rt60
 
-    def measure_rt60(self, decay_db=60, plot=False):
+    def measure_rt60(
+        self, decay_db=60, plot=False, label=None, linear_domain_fit=False
+    ):
         """
         Measures the reverberation time (RT60) of the simulated RIR.
 
@@ -2984,6 +2973,8 @@ def measure_rt60(self, decay_db=60, plot=False):
             dB, the RT60 is twice the time measured.
         plot: bool
             Displays a graph of the Schroeder curve and the estimated RT60.
+        label: str, optional
+            Optional label to use in a plot.
 
         Returns
         -------
@@ -2995,8 +2986,14 @@ def measure_rt60(self, decay_db=60, plot=False):
 
         for m in range(self.n_mics):
             for s in range(self.n_sources):
+                label = f"{label}-{m}-{s}" if label else ""
                 rt60[m, s] = measure_rt60(
-                    self.rir[m][s], fs=self.fs, plot=plot, decay_db=decay_db
+                    self.rir[m][s],
+                    fs=self.fs,
+                    plot=plot,
+                    decay_db=decay_db,
+                    label=label,
+                    linear_domain_fit=linear_domain_fit,
                 )
 
         return rt60
@@ -3317,9 +3314,6 @@ def __init__(
         # Anyways, we use the energy_absorption and scattering corresponding to an anechoic room.
         materials = Material(energy_absorption=1.0, scattering=0.0)
 
-        # Set deprecated parameter
-        absorption = None
-
         ShoeBox.__init__(
             self,
             p=p,