From 0777cf54e8b5fd4be4ae9810833148f58d24dee9 Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sat, 24 Jan 2026 00:09:42 +0900 Subject: [PATCH 1/9] Minor refactoring in libroom. --- pyroomacoustics/libroom_src/common.hpp | 5 +++++ pyroomacoustics/libroom_src/libroom.cpp | 1 + pyroomacoustics/libroom_src/room.cpp | 10 +++++++--- 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/pyroomacoustics/libroom_src/common.hpp b/pyroomacoustics/libroom_src/common.hpp index 87c82d69..8d1ef382 100644 --- a/pyroomacoustics/libroom_src/common.hpp +++ b/pyroomacoustics/libroom_src/common.hpp @@ -169,6 +169,11 @@ class Histogram2D { return array; } + + Eigen::ArrayXXi get_counts() const + { + return counts; + } }; #endif // __COMMON_HPP__ diff --git a/pyroomacoustics/libroom_src/libroom.cpp b/pyroomacoustics/libroom_src/libroom.cpp index 444c6be2..a01f8d1b 100644 --- a/pyroomacoustics/libroom_src/libroom.cpp +++ b/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 diff --git a/pyroomacoustics/libroom_src/room.cpp b/pyroomacoustics/libroom_src/room.cpp index 1a5521c6..498506b0 100644 --- a/pyroomacoustics/libroom_src/room.cpp +++ b/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::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::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::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) From fa46a5b00b1c31d5dea277aaa5b9e0649a7f4946 Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sat, 24 Jan 2026 00:22:52 +0900 Subject: [PATCH 2/9] Fixes the DC shelf filter in the octave bands and the list of bandwidths. --- pyroomacoustics/acoustics.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/pyroomacoustics/acoustics.py b/pyroomacoustics/acoustics.py index 8a5bb5a5..1f165030 100644 --- a/pyroomacoustics/acoustics.py +++ b/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 ---------- From d8cb3f57be22bbadad34262b6d5ac01438d3f588 Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sat, 24 Jan 2026 00:25:23 +0900 Subject: [PATCH 3/9] Rewrites Sabine's formula in room.py to be the harmonic mean over walls and frequency bands. --- pyroomacoustics/room.py | 60 +++++++++++++++++++---------------------- 1 file changed, 27 insertions(+), 33 deletions(-) diff --git a/pyroomacoustics/room.py b/pyroomacoustics/room.py index 528efd58..b56f1f88 100644 --- a/pyroomacoustics/room.py +++ b/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 not self.air_absorption 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, From c88e6d0162f224be930ff35db55ac80f054e5e92 Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sat, 24 Jan 2026 00:27:31 +0900 Subject: [PATCH 4/9] Use fit in RT60 computation function. --- pyroomacoustics/experimental/rt60.py | 90 ++++++++++++++++++++++++---- 1 file changed, 77 insertions(+), 13 deletions(-) diff --git a/pyroomacoustics/experimental/rt60.py b/pyroomacoustics/experimental/rt60.py index d9df7ee2..6458256a 100644 --- a/pyroomacoustics/experimental/rt60.py +++ b/pyroomacoustics/experimental/rt60.py @@ -30,10 +30,58 @@ .. [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 measure_rt60(h, fs=1, decay_db=60, energy_thres=1.0, plot=False, rt60_tgt=None): +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 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 +105,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 +137,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 +172,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() From ad57c5e5f96627098fe58276d51a2e7450cf0eab Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sat, 24 Jan 2026 00:30:20 +0900 Subject: [PATCH 5/9] Fixes the RT histogram to rir function. --- pyroomacoustics/simulation/rt.py | 42 +-- .../tests/test_rt_multiband_energy.py | 286 ++++++++++++++++++ 2 files changed, 309 insertions(+), 19 deletions(-) create mode 100644 pyroomacoustics/tests/test_rt_multiband_energy.py diff --git a/pyroomacoustics/simulation/rt.py b/pyroomacoustics/simulation/rt.py index 8ee6f3c2..4abb3d53 100644 --- a/pyroomacoustics/simulation/rt.py +++ b/pyroomacoustics/simulation/rt.py @@ -30,6 +30,7 @@ import numpy as np from scipy.interpolate import interp1d +from scipy.signal import fftconvolve def sequence_generation(volume, duration, c, fs, max_rate=10000): @@ -78,11 +79,10 @@ def interp_hist(hist, N): return out -def seq_bin_power(seq, hbss): - seq_rot = seq.reshape((-1, hbss)) # shape 72,64 - power = np.sum(seq_rot**2, axis=1) - power[power <= 0.0] = 1.0 - return power +def seq_rolling_power(seq, hbss, filter_length_mult=1): + """Smooth local energy of the sequence.""" + filter = np.ones(filter_length_mult * hbss) / filter_length_mult + return fftconvolve(seq**2, filter, mode="same") def compute_rt_rir( @@ -96,6 +96,7 @@ def compute_rt_rir( octave_bands, air_abs_coeffs=None, ): + # get the maximum length from the histograms # Sum vertically across octave band for each value in # histogram (7,2500) -> (2500) -> np .nonzero( @@ -125,37 +126,40 @@ def compute_rt_rir( # this is the random sequence for the tail generation seq = sequence_generation(volume_room, N / fs, c, fs) seq = seq[:N] # take values according to N as seq is larger + seq_original_power = np.mean(seq**2) n_bands = histograms[0].shape[0] bws = octave_bands.get_bw() if n_bands > 1 else [fs / 2] + hist_bands = histograms[0][:, :n_bins] + rir_bands = np.zeros((n_bands, N)) - for b, bw in enumerate(bws): # Loop through every band + for b, bw in enumerate(np.sort(bws)): # Loop through every band if n_bands > 1: seq_bp = octave_bands.analysis(seq, band=b) else: seq_bp = seq.copy() - # Take only those bins which have some non-zero values for that specific octave bands. - hist = histograms[0][b, :n_bins] + # This accounts for the relative energy weight of this band in the full sequence. + band_weight = np.mean(seq_bp**2) / seq_original_power - # we normalize the histogram by the sequence power in that bin - seq_power = seq_bin_power(seq_bp, hbss) + # We normalize the histogram by the sequence power in that bin. + seq_power = seq_rolling_power(seq_bp, hbss, filter_length_mult=2) + nonzero = seq_power > 0.0 + seq_power_norm = np.where(nonzero, 1.0 / np.where(nonzero, seq_power, 1.0), 0.0) - # we linarly interpolate the histogram to smoothly cover all the samples - hist_lin_interp = interp_hist(hist / seq_power, N) + # We linarly interpolate the histogram to smoothly cover all the samples. + hist = interp_hist(hist_bands[b], N) - # Normalize the band power - # the (bw / fs * 2.0) is necessary to normalize the band power - # this is the contribution of the octave band to total energy - seq_bp *= np.sqrt(bw / fs * 2.0 * hist_lin_interp) + # Compute the gain to be applied to the sequence. + gain = np.sqrt(hist * seq_power_norm * band_weight) - # Impulse response for every octave band for each microphone - rir_bands[b] = seq_bp + # Impulse response for every octave band for each microphone. + rir_bands[b] = seq_bp * gain if air_abs_coeffs is not None: if rir_bands.shape[0] == 1: - # do the octave band analysis if it was not done in the first step + # Do the octave band analysis if it was not done in the first step. rir_bands = np.array( [ octave_bands.analysis(rir_bands[0], band=b) diff --git a/pyroomacoustics/tests/test_rt_multiband_energy.py b/pyroomacoustics/tests/test_rt_multiband_energy.py new file mode 100644 index 00000000..9589bde2 --- /dev/null +++ b/pyroomacoustics/tests/test_rt_multiband_energy.py @@ -0,0 +1,286 @@ +import numpy as np +import matplotlib.pyplot as plt +import pyroomacoustics as pra +import pytest +from scipy.optimize import curve_fit + +pra.constants.set("octave_bands_keep_dc", True) + + +@pytest.fixture +def samplerate(): + return 16_000 + + +def get_scene_geometry(): + room_dim = np.array([3.0, 4.0, 5.0]) + source_loc = room_dim / 2.0 + np.array([0.001, -0.002, 0.003]) + co, az = np.pi / 2.5, 0.1 * np.pi + unit_vec = np.array([np.cos(az) * np.sin(co), np.sin(az) * np.sin(co), np.cos(co)]) + mic_loc = source_loc + 1.5 * unit_vec + return room_dim, source_loc, mic_loc + + +@pytest.fixture +def scene_geometry(): + return get_scene_geometry() + + +def get_wall_material(multiband=True): + # Material with maximum energy absorption + if multiband: + full_eabs = { + "description": "Multi-band fully absorbant", + "coeffs": [0.11, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16], + "center_freqs": 125 * 2.0 ** np.arange(7), + } + return pra.Material(energy_absorption=full_eabs) + else: + return pra.Material(energy_absorption=0.15) + + +@pytest.fixture(params=[(True,), (False,)], ids=["Multiband", "Single-band"]) +def wall_material(request): + return get_wall_material(multiband=request.param) + + +def energy_hist_2_rt60(hist, hist_bin_size, decay_db): + # Integrate the energy from the tail. + energy = np.cumsum(hist[::-1])[::-1] + schroeder = 10.0 * np.log10(energy) + # Remove the first 5 dB. + schroeder += 5.0 - schroeder[0] + + t0 = np.where(schroeder < 0.0)[0][0] + t1 = np.where(schroeder < -decay_db)[0][0] + N = t1 - t0 + + data = schroeder[t0:t1] + data -= data[0] + + t = np.arange(N) * hist_bin_size + X = np.column_stack((t, np.ones(N))) + p, *_ = np.linalg.lstsq(X, data) + + rt60_1 = -60.0 / p[0] + + # 1. Define your model function + def model(x, a, b): + return b * np.exp(a * x) + + def jac(x, a, b): + u = np.exp(a * x) + return np.column_stack((x * b * u, u)) + + # 2. Provide an initial guess [a, b, c] + # This is crucial for convergence! + initial_guess = [p[0] / (10.0 * np.log10(np.e)), 10.0 ** (p[1] / 10.0)] + + # 3. Perform the fit + popt, pcov, info_dict, *_ = curve_fit( + model, t, 10.0 ** (data / 10.0), p0=initial_guess, jac=jac, full_output=True + ) + + rt60_2 = -6.0 * np.log(10.0) / popt[0] + + return rt60_1, rt60_2 + + +def simulate(room_dim, source_loc, mic_loc, material, use_rt, fs): + hist_bin_size = 0.004 + + room = pra.ShoeBox( + room_dim, + fs=fs, + materials=material, + max_order=-1 if use_rt else 100, + ray_tracing=use_rt, + ) + if use_rt: + room.set_ray_tracing( + n_rays=250_000, receiver_radius=0.5, hist_bin_size=hist_bin_size + ) + + room.add_source(source_loc, signal=np.zeros(10)) + room.add_microphone(mic_loc) + + room.simulate() + rir = room.rir[0][0] + + bin_size = int(hist_bin_size * fs) + n_bins = rir.shape[0] // bin_size + n_max = n_bins * bin_size + rir_hist = rir[:n_max] + rir_hist = rir_hist.reshape((n_bins, bin_size)) + hist = np.sum(rir_hist**2, axis=1) + + bands_energy = np.mean(room.octave_bands.analysis(rir) ** 2, axis=0) + bands_energy = bands_energy / bands_energy.sum() + norm_bw = room.octave_bands.get_bw() / (room.fs / 2.0) + + rt60 = room.measure_rt60(decay_db=20.0)[0][0] + N60 = int(rt60 * room.fs) + + rir_energy = np.sum(rir[:N60] ** 2) + + return room, rir[:N60], hist, bands_energy, norm_bw, rt60, rir_energy + + +def compute_tail_decay_curve(t0, rir, fs, hist_bin_size, n_max): + rir = rir[t0:n_max] + bin_size = int(hist_bin_size * fs) + n_bins = rir.shape[0] // bin_size + n_max = n_bins * bin_size + rir_hist = rir[:n_max] + rir_hist = rir_hist.reshape((n_bins, bin_size)) + return np.sum(rir_hist**2, axis=1) + + +def ism_to_rt_energy_ratio_db(t0, rir_rt, rir_ism, fs, hist_bin_size=0.004): + n_max = min(rir_rt.shape[0], rir_ism.shape[0]) + hist_rt = compute_tail_decay_curve(t0, rir_rt, fs, hist_bin_size, n_max) + hist_ism = compute_tail_decay_curve(t0, rir_ism, fs, hist_bin_size, n_max) + return 10.0 * np.log10(np.mean(hist_ism**2) / np.mean((hist_rt - hist_ism) ** 2)) + + +def test_energy_decay_ism_vs_rt(samplerate, scene_geometry, wall_material): + np.random.seed(42) + + room_dim, source_loc, mic_loc = scene_geometry + fs = samplerate + + room, rir_ism, *_ = simulate( + room_dim, source_loc, mic_loc, wall_material, use_rt=False, fs=fs + ) + room, rir_rt, *_ = simulate( + room_dim, source_loc, mic_loc, wall_material, use_rt=True, fs=fs + ) + + c = pra.constants.get("c") + t0 = int( + fs * np.linalg.norm(source_loc - mic_loc) / c + + pra.constants.get("frac_delay_length") + ) + + # Compare the energy decay only for the full band, because low + # frequencies are not reliable with a shoebox room. + irer_db = ism_to_rt_energy_ratio_db(t0, rir_ism, rir_rt, fs, hist_bin_size=0.04) + assert irer_db >= 10.0, f"Failed IRER {irer_db:.2f} dB < 10.0" + + bws = room.octave_bands.get_bw() + tols = [0.25, 0.1, 0.08] + for b, bw in enumerate(np.sort(bws)): # Loop through every band + rir_band_ism = room.octave_bands.analysis(rir_ism, band=b) + rir_band_rt = room.octave_bands.analysis(rir_rt, band=b) + + # Check the RT60 are within 15% of each other. + rt60_band_ism = pra.experimental.measure_rt60(rir_band_ism, fs=fs) + rt60_band_rt = pra.experimental.measure_rt60(rir_band_rt, fs=fs) + rt60_delta = np.abs(rt60_band_ism - rt60_band_rt) + rt60_rel_error = rt60_delta / rt60_band_ism + tol = tols[min(b, 2)] + assert rt60_rel_error < tol, f"Failed {rt60_rel_error=:.3} > {tol}" + + +if __name__ == "__main__": + # Room dimensions + fs = 16_000 + hist_bin_size = 0.004 + room_dim, source_loc, mic_loc = get_scene_geometry() + material = get_wall_material(multiband=True) + + c = pra.constants.get("c") + t0 = int( + fs * np.linalg.norm(source_loc - mic_loc) / c + + pra.constants.get("frac_delay_length") + ) + + fig, axes = plt.subplots(2, 2, figsize=(12.0, 8.0)) + + rirs = {} + + for label, use_rt in {"RT": True, "ISM": False}.items(): + + room, rir, hist, bands_energy, norm_bw, rt60, rir_energy = simulate( + room_dim, source_loc, mic_loc, material, use_rt=use_rt, fs=fs + ) + rirs[label] = rir + + d = int(t0 / fs / room.rt_args["hist_bin_size"]) + d = 0 + et60_1, et60_2 = energy_hist_2_rt60( + hist[d + 1 :], room.rt_args["hist_bin_size"], 20.0 + ) + + print(f"{label}:") + print(f" T60={rt60:.3f}") + print(f" ET60 1={et60_1:.3f}") + print(f" ET60 2={et60_2:.3f}") + for formula in ["sabine", "eyring"]: + rt60_th = room.rt60_theory(formula) + print(f" T60={rt60_th:.3} ({formula})") + + print(f" Energy(RIR)={rir_energy:.5f}") + + axes[0, 0].plot( + np.arange(hist.shape[0]) * hist_bin_size, hist.T**0.5, label=label + ) + + axes[0, 1].plot(np.arange(rir.shape[0]) / room.fs, rir, label=label) + + axes[1, 0].magnitude_spectrum(rir, Fs=room.fs, scale="dB") + + axes[1, 1].plot(bands_energy, label=f"{label} measured") + axes[1, 1].plot(norm_bw, label=f"{label} expected") + + irer_fb = ism_to_rt_energy_ratio_db( + t0, rirs["ISM"], rirs["RT"], room.fs, hist_bin_size=0.04 + ) + print(f"Full band IRER: {irer_fb:.2f} dB") + bws = room.octave_bands.get_bw() + fig2, axes2 = plt.subplots(2, len(bws)) + for b, bw in enumerate(np.sort(bws)): # Loop through every band + # Compute the energy difference between the RIR at different bands. + rir_band_ism = room.octave_bands.analysis(rirs["ISM"][t0:], band=b) + rir_band_rt = room.octave_bands.analysis(rirs["RT"][t0:], band=b) + n_max = min(rir_band_ism.shape[0], rir_band_rt.shape[0]) + + rt60_band_ism = pra.experimental.measure_rt60( + rir_band_ism, fs=fs, decay_db=60.0 + ) + rt60_band_rt = pra.experimental.measure_rt60(rir_band_rt, fs=fs, decay_db=60.0) + rt60_delta = np.abs(rt60_band_ism - rt60_band_rt) + rt60_rel_error = rt60_delta / rt60_band_ism + print( + f"Band {b}: RT60 error abs={rt60_delta:.3f} rel={rt60_rel_error:.3f}. ISM: {rt60_band_ism:.3f} RT: {rt60_band_rt:.3f}" + ) + + for lbl, ir in zip(["ism", "rt"], [rir_band_ism, rir_band_rt]): + axes2[0, b].plot( + np.arange(ir.shape[0]) / fs, ir, label=lbl if b == 0 else None + ) + axes2[0, b].set_title(f"band {bw:.0f}") + + hist_band = compute_tail_decay_curve(t0, ir, fs, 0.04, n_max) + axes2[1, b].plot( + np.arange(hist_band.shape[0]) * 0.04, + hist_band, + label=lbl if b == 0 else None, + ) + axes2[1, b].set_title(f"band energy {bw:.0f}") + + irer_db = ism_to_rt_energy_ratio_db( + 0, rir_band_ism, rir_band_rt, fs, hist_bin_size=0.04 + ) + print(f"Band {b}: IRER={irer_db:.2f} dB") + + fig2.legend() + + axes[0, 0].legend() + axes[0, 1].set_title("RIR") + axes[0, 0].set_title("Magnitude histogram") + axes[1, 0].set_title("Magnitude spectrum of the RIR") + axes[1, 1].set_title("Band energies") + axes[1, 1].legend() + + plt.show() From 4b4eeba2026b18f351c69d184f5f49fccb5fab57 Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sat, 24 Jan 2026 00:37:41 +0900 Subject: [PATCH 6/9] Updates the RT60 measurement in an example. --- CHANGELOG.rst | 10 ++++++++++ examples/room_shoebox_3d_rt.py | 2 +- pyroomacoustics/room.py | 2 +- 3 files changed, 12 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 7ac40b28..898701ca 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -13,6 +13,16 @@ adheres to `Semantic Versioning `_. Bugfix: ~~~~~~ +- Fixes the computation of the RIR based on ray tracing. In particular the sequence + energy shaping and the band weights. + +- Fixes the lowest band shelf filter in the octave bands. + +- Fixes the computation of the octave band widths. + +- Fixes `pra.experimental.measure_rt60`: Compute the T60 using a fit. Default is + log-domain. Adds option to fit in linear domain. + - In `doa.py`, the `ax.xaxis.grid` and `ax.yaxis.grid` parameters were changed from `b` to `visible`. - Fixes MicrophoneArray.append(MicrophoneArray) diff --git a/examples/room_shoebox_3d_rt.py b/examples/room_shoebox_3d_rt.py index 665d4da8..793c80c2 100644 --- a/examples/room_shoebox_3d_rt.py +++ b/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") diff --git a/pyroomacoustics/room.py b/pyroomacoustics/room.py index b56f1f88..72d0da39 100644 --- a/pyroomacoustics/room.py +++ b/pyroomacoustics/room.py @@ -2951,7 +2951,7 @@ def rt60_theory(self, formula="sabine"): # and walls. inv_rt60 = 0.0 for i, bw in enumerate(bandwidths): - m = 0.0 if not self.air_absorption else self.air_absorption[i] + m = 0.0 if self.air_absorption is None else self.air_absorption[i] for w in self.walls: inv_rt60 += bw / rt60_fn(w.area(), V, w.absorption[i], m, c) From a0d2c17184f3683783556900aecccf517d8ee329 Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sat, 24 Jan 2026 00:38:33 +0900 Subject: [PATCH 7/9] Lint --- pyroomacoustics/experimental/rt60.py | 1 + pyroomacoustics/tests/test_rt_multiband_energy.py | 5 +++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/pyroomacoustics/experimental/rt60.py b/pyroomacoustics/experimental/rt60.py index 6458256a..d3e18288 100644 --- a/pyroomacoustics/experimental/rt60.py +++ b/pyroomacoustics/experimental/rt60.py @@ -31,6 +31,7 @@ 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 diff --git a/pyroomacoustics/tests/test_rt_multiband_energy.py b/pyroomacoustics/tests/test_rt_multiband_energy.py index 9589bde2..0aa93f63 100644 --- a/pyroomacoustics/tests/test_rt_multiband_energy.py +++ b/pyroomacoustics/tests/test_rt_multiband_energy.py @@ -1,9 +1,10 @@ -import numpy as np import matplotlib.pyplot as plt -import pyroomacoustics as pra +import numpy as np import pytest from scipy.optimize import curve_fit +import pyroomacoustics as pra + pra.constants.set("octave_bands_keep_dc", True) From adfdf5924463e0fb771ea819111a601423360211 Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sun, 25 Jan 2026 14:14:02 +0900 Subject: [PATCH 8/9] Revert to using theoretical weights for the bands. --- pyroomacoustics/simulation/rt.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pyroomacoustics/simulation/rt.py b/pyroomacoustics/simulation/rt.py index 4abb3d53..52628450 100644 --- a/pyroomacoustics/simulation/rt.py +++ b/pyroomacoustics/simulation/rt.py @@ -126,7 +126,6 @@ def compute_rt_rir( # this is the random sequence for the tail generation seq = sequence_generation(volume_room, N / fs, c, fs) seq = seq[:N] # take values according to N as seq is larger - seq_original_power = np.mean(seq**2) n_bands = histograms[0].shape[0] bws = octave_bands.get_bw() if n_bands > 1 else [fs / 2] @@ -140,9 +139,6 @@ def compute_rt_rir( else: seq_bp = seq.copy() - # This accounts for the relative energy weight of this band in the full sequence. - band_weight = np.mean(seq_bp**2) / seq_original_power - # We normalize the histogram by the sequence power in that bin. seq_power = seq_rolling_power(seq_bp, hbss, filter_length_mult=2) nonzero = seq_power > 0.0 @@ -151,6 +147,9 @@ def compute_rt_rir( # We linarly interpolate the histogram to smoothly cover all the samples. hist = interp_hist(hist_bands[b], N) + # This accounts for the relative energy weight of this band in the full sequence. + band_weight = bw / fs * 2.0 + # Compute the gain to be applied to the sequence. gain = np.sqrt(hist * seq_power_norm * band_weight) From c972d739ddc7b615ca8d174c80652817c76582bf Mon Sep 17 00:00:00 2001 From: Robin Scheibler Date: Sun, 25 Jan 2026 21:47:21 +0900 Subject: [PATCH 9/9] Changelog --- CHANGELOG.rst | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index fe66ab31..cd824b1b 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -10,23 +10,29 @@ adheres to `Semantic Versioning `_. `Unreleased`_ ------------- -Bugfix: -~~~~~~ - -- Fixes the computation of the RIR based on ray tracing. In particular the sequence - energy shaping and the band weights. -- Fixes the lowest band shelf filter in the octave bands. +Changed +~~~~~~~ -- Fixes the computation of the octave band widths. +- Refactor the way the RIR is weighted with the histogram in simulation/rt.py. -- Fixes `pra.experimental.measure_rt60`: Compute the T60 using a fit. Default is +- Improves ``pra.experimental.measure_rt60``: Compute the T60 using a fit. Default is log-domain. Adds option to fit in linear domain. -- In `doa.py`, the `ax.xaxis.grid` and `ax.yaxis.grid` parameters were changed from `b` to `visible`. +Bugfix +~~~~~~ + +- 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.