From 2d3627ce2c68da7f7a5bad6d8ff95773510402de Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 20:11:10 +0100 Subject: [PATCH 1/2] [FEATURE] EDC noise compensation methods and T20/T30 uncertainty RoomAnalysis takes edcMethod, one of the noise compensations of ISO 3382-1 compared by Guski and Vorlaender (2013): 'cutWithCorrection', the method used so far and the default, 'justCut', 'noCut' and 'subtractNoise' (Chu, 1978). bypassLundeby=True is the 'noCut' method; it raised UnboundLocalError because the Lundeby parameters were returned without being set. RoomAnalysis.uncertainty gives the standard deviation of T20 or T30 by band from ISO 3382-1:2009, Equations (4) and (5), for N positions and n decays per position, with n = 10 for an impulse response (clause 7.2). Tests use a synthetic decay with T = 1 s over background noise 45 dB below its start: the compensated methods find T, 'noCut' overestimates it, 'justCut' falls below 'cutWithCorrection', bypassLundeby equals 'noCut', and the uncertainty follows the equations of the standard. Co-Authored-By: Claude Fable 5.1 --- pytta/classes/analysis.py | 83 +++++++++++++++++++++++++++++++-------- tests/test_room_edc.py | 81 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 148 insertions(+), 16 deletions(-) create mode 100644 tests/test_room_edc.py diff --git a/pytta/classes/analysis.py b/pytta/classes/analysis.py index 7f15924..3e8f24b 100644 --- a/pytta/classes/analysis.py +++ b/pytta/classes/analysis.py @@ -957,7 +957,8 @@ def __init__(self, ir: SignalObj, nthOct: int = 1, plotLundeby: bool = False, bypassLundeby: bool = False, suppressWarnings: bool = True, - ircut: float = None, **kwargs): + ircut: float = None, + edcMethod: str = 'cutWithCorrection', **kwargs): _ir = ir.IR if type(ir) == ImpulsiveResponse else ir minBand = freq_to_band(minFreq, nthOct, 1000, 10) maxBand = freq_to_band(maxFreq, nthOct, 1000, 10) @@ -966,6 +967,7 @@ def __init__(self, ir: SignalObj, nthOct: int = 1, self.ir = crop_IR(_ir, ircut) self._params, self.listEDC, self.fhSignal = self.estimate_energy_parameters(self.ir, self.bands, plotLundeby, bypassLundeby, suppressWarnings, + edcMethod=edcMethod, nthOct=nthOct, minFreq=minFreq, maxFreq=maxFreq) return @@ -974,7 +976,9 @@ def __init__(self, ir: SignalObj, nthOct: int = 1, def estimate_energy_parameters(ir: SignalObj, bands: np.ndarray, plotLundeby: bool = False, bypassLundeby: bool = False, - suppressWarnings: bool = False, **kwargs): + suppressWarnings: bool = False, + edcMethod: str = 'cutWithCorrection', + **kwargs): """ Estimate the Impulse Response energy parameters. @@ -993,7 +997,8 @@ def estimate_energy_parameters(ir: SignalObj, bands: np.ndarray, A dict with parameters by name. """ - listEDC, fhSignal = cumulative_integration(ir, bypassLundeby, plotLundeby, suppressWarnings, **kwargs) + listEDC, fhSignal = cumulative_integration(ir, bypassLundeby, plotLundeby, suppressWarnings, + edcMethod=edcMethod, **kwargs) params = {} params['rms'] = fhSignal.rms() params['SPL'] = fhSignal.spl() @@ -1063,6 +1068,32 @@ def T30(self): """Reverberation time with 30 dB decay, by frequency `band`.""" return self._params['T30'] + def uncertainty(self, param: str = 'T30', numPositions: int = 1, + numDecays: int = 10) -> np.ndarray: + """ + Standard deviation of T20 or T30 by band, from ISO 3382-1:2009, + Equations (4) and (5): sigma(T) = k T sqrt((1 + c / n) / (N B T)), + with k, c = 0.88, 1.90 for T20 and 0.55, 1.52 for T30. + + * param ('T30'), (str): 'T20' or 'T30'; + * numPositions (1), (int): N, independent source and receiver + positions averaged; + * numDecays (10), (int): n, decays measured in each position; an + impulse response counts as n = 10 (clause 7.2). + + B is the bandwidth of the band filter, fc (2**(1/(2 nthOct)) - + 2**(-1/(2 nthOct))): 0.71 fc for octaves and 0.23 fc for thirds. + """ + if param not in ('T20', 'T30'): + raise ValueError("param must be 'T20' or 'T30'.") + k, c = (0.88, 1.90) if param == 'T20' else (0.55, 1.52) + rt = self._params[param] + half = 2 ** (1 / (2 * self.nthOct)) + B = np.asarray(self.bands) * (half - 1 / half) + with np.errstate(divide='ignore', invalid='ignore'): + return k * rt * np.sqrt((1 + c / numDecays) + / (numPositions * B * rt)) + # @property # def BR(self): # """Reverberation time with 30 dB decay, by frequency `band`.""" @@ -1401,9 +1432,31 @@ def _Lundeby_correction(band, timeSignal, samplingRate, numSamples, # @njit def energy_decay_calculation(band, timeSignal, timeVector, samplingRate, numSamples, numChannels, timeLength, bypassLundeby, - suppressWarnings=True): - """Calculate the Energy Decay Curve.""" - if bypassLundeby is False: + suppressWarnings=True, + edcMethod='cutWithCorrection'): + """ + Calculate the Energy Decay Curve by backward integration, with one of + the noise compensations of ISO 3382-1 (Guski and Vorlaender, 2013): + + * 'cutWithCorrection': truncate at the intersection of the decay + with the background noise (Lundeby et al., 1995) and add the + energy of the decay beyond it; + * 'justCut': truncate at the intersection, with no correction; + * 'noCut': integrate the whole response, noise included; + * 'subtractNoise': subtract the background noise level from the + squared response and integrate the whole of it (Chu, 1978). + + bypassLundeby=True is the 'noCut' method. + """ + methods = ('cutWithCorrection', 'justCut', 'noCut', 'subtractNoise') + if edcMethod not in methods: + raise ValueError(f'edcMethod must be one of {methods}.') + if bypassLundeby: + edcMethod = 'noCut' + if edcMethod == 'noCut': + lundebyParams = None + interIdx, lateRT, BGL = 0, 1, 0 + else: lundebyParams = \ _Lundeby_correction(band, timeSignal, @@ -1414,22 +1467,18 @@ def energy_decay_calculation(band, timeSignal, timeVector, samplingRate, suppressWarnings=suppressWarnings) _, c1, interIdx, BGL = lundebyParams lateRT = -60/c1 if c1 != 0 else 0 - else: - interIdx = 0 - lateRT = 1 - - if interIdx == 0: + if interIdx == 0 or edcMethod == 'subtractNoise': interIdx = -1 - truncatedTimeSignal = timeSignal[:interIdx, 0] truncatedTimeVector = timeVector[:interIdx] - if lateRT != 0.0: - if not bypassLundeby: + if edcMethod == 'cutWithCorrection': C = samplingRate*BGL*lateRT/(6*np.log(10)) else: C = 0 sqrInv = truncatedTimeSignal[::-1]**2 + if edcMethod == 'subtractNoise': + sqrInv = sqrInv - BGL energyDecayFull = np.cumsum(sqrInv)[::-1] + C energyDecay = energyDecayFull/energyDecayFull[0] else: @@ -1444,6 +1493,7 @@ def cumulative_integration(inputSignal, bypassLundeby, plotLundebyResults, suppressWarnings=True, + edcMethod='cutWithCorrection', **kwargs): """Cumulative integration with proper corrections.""" @@ -1485,9 +1535,10 @@ def plot_lundeby(): signal.numChannels, signal.timeLength, bypassLundeby, - suppressWarnings=suppressWarnings) + suppressWarnings=suppressWarnings, + edcMethod=edcMethod) listEDC.append((energyDecay, energyVector)) - if plotLundebyResults: # Placed here because Numba can't handle plots. + if plotLundebyResults and lundebyParams is not None: # plot_lundeby(band, timeVector, timeSignal, samplingRate, # lundebyParams) plot_lundeby() diff --git a/tests/test_room_edc.py b/tests/test_room_edc.py new file mode 100644 index 0000000..560565f --- /dev/null +++ b/tests/test_room_edc.py @@ -0,0 +1,81 @@ +""" +Tests for the energy decay curve methods and the reverberation time +uncertainty of pytta.RoomAnalysis. + +The impulse response is a synthetic exponential decay of white noise with +T = 1 s, over a steady background noise 45 dB below its start, so the +methods of noise compensation of ISO 3382-1 (Guski and Vorlaender, 2013) +can be told apart. Each test pins one rule. +""" +import os +import unittest +import warnings + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import numpy as np + +import pytta + +FS = 48000 +T = 1.0 + + +def decay_with_noise(): + rng = np.random.default_rng(0) + t = np.arange(int(2.5 * FS)) / FS + h = rng.standard_normal(t.size) * 10 ** (-3 * t / T) \ + + 10 ** (-45 / 20) * rng.standard_normal(t.size) + return pytta.SignalObj(h, 'time', FS, signalType='energy') + + +def analyse(ir, **kwargs): + with warnings.catch_warnings(): + warnings.simplefilter('ignore') + return pytta.RoomAnalysis(ir, nthOct=1, minFreq=500., maxFreq=2000., + **kwargs) + + +class TestRoomEDC(unittest.TestCase): + + @classmethod + def setUpClass(cls): + ir = decay_with_noise() + cls.T30 = {m: analyse(ir, edcMethod=m).T30 for m in + ('cutWithCorrection', 'justCut', 'noCut', 'subtractNoise')} + cls.bypass = analyse(ir, bypassLundeby=True) + for m, v in cls.T30.items(): + print(f'\n {m:18s} T30 = {np.round(v, 3)}') + + def test_compensated_methods_find_the_true_decay(self): + for m in ('cutWithCorrection', 'subtractNoise'): + self.assertTrue(np.allclose(self.T30[m], T, rtol=0.05), m) + + def test_uncompensated_methods_are_biased_as_expected(self): + # The background noise flattens the tail of the full curve, and a + # truncation without correction makes the curve fall too fast. + self.assertTrue(np.all(self.T30['noCut'] > 1.10 * T)) + self.assertTrue(np.all(self.T30['justCut'] + < self.T30['cutWithCorrection'])) + + def test_bypass_lundeby_is_the_no_cut_method(self): + self.assertTrue(np.array_equal(self.bypass.T30, self.T30['noCut'])) + + def test_uncertainty_follows_iso_3382_1(self): + # ISO 3382-1:2009, Equations (4) and (5), with the bandwidth of an + # octave filter B = 0.71 fc and n = 10 for an impulse response + # (clause 7.2). + ra = analyse(decay_with_noise()) + fc = np.array([500., 1000., 2000.]) + for param, k, c in (('T20', 0.88, 1.90), ('T30', 0.55, 1.52)): + rt = getattr(ra, param) + for N in (1, 4): + iso = k * rt * np.sqrt((1 + c / 10) / (N * 0.71 * fc * rt)) + got = ra.uncertainty(param, numPositions=N) + self.assertTrue(np.allclose(got, iso, rtol=0.005), param) + with self.assertRaises(ValueError): + ra.uncertainty('EDT') + + +if __name__ == '__main__': + unittest.main() From 4d037ac2b2e409cf5853541c2851a4aae19bb28d Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 21:29:36 +0100 Subject: [PATCH 2/2] [FEATURE] peak-to-noise ratio and timbre ratios in RoomAnalysis RoomAnalysis.PSNR gives, by band, the peak of the squared band response over the background noise level found by Lundeby (nan without it), and RoomAnalysis.timbre_ratios the bass and treble ratios of Beranek from the reverberation times of the octave bands. cumulative_integration also returns the Lundeby parameters of each band. Tests check that the PSNR grows by 10 dB when the noise floor drops by 10 dB, and BR and TR of a synthetic response with known decay per band. Co-Authored-By: Claude Fable 5.1 --- pytta/classes/analysis.py | 37 +++++++++++++++++++++++++++++++--- tests/test_room_edc.py | 42 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+), 3 deletions(-) diff --git a/pytta/classes/analysis.py b/pytta/classes/analysis.py index 3e8f24b..33326fb 100644 --- a/pytta/classes/analysis.py +++ b/pytta/classes/analysis.py @@ -997,8 +997,9 @@ def estimate_energy_parameters(ir: SignalObj, bands: np.ndarray, A dict with parameters by name. """ - listEDC, fhSignal = cumulative_integration(ir, bypassLundeby, plotLundeby, suppressWarnings, - edcMethod=edcMethod, **kwargs) + listEDC, fhSignal, listLundeby = cumulative_integration( + ir, bypassLundeby, plotLundeby, suppressWarnings, + edcMethod=edcMethod, **kwargs) params = {} params['rms'] = fhSignal.rms() params['SPL'] = fhSignal.spl() @@ -1010,6 +1011,13 @@ def estimate_energy_parameters(ir: SignalObj, bands: np.ndarray, params['EDT'] = reverberation_time('EDT', listEDC) params['T20'] = reverberation_time(20, listEDC) params['T30'] = reverberation_time(30, listEDC) + # peak-to-noise ratio: peak of the squared band response over the + # background noise level found by Lundeby (nan without it) + noise = np.array([p[3] if p is not None and p[3] > 0 else np.nan + for p in listLundeby], dtype=float) + with np.errstate(divide='ignore', invalid='ignore'): + params['PSNR'] = 10 * np.log10( + np.max(fhSignal.timeSignal ** 2, axis=0) / noise) # self._params['BR'], self._params['TR'] = timbre_ratios(self.T20) return params, listEDC, fhSignal @@ -1068,6 +1076,27 @@ def T30(self): """Reverberation time with 30 dB decay, by frequency `band`.""" return self._params['T30'] + @property + def PSNR(self): + """Peak-to-noise ratio by band, in dB: the peak of the squared band + response over the background noise level found by Lundeby; nan + with edcMethod='noCut'.""" + return self._params['PSNR'] + + def timbre_ratios(self, param: str = 'T30'): + """ + Bass ratio and treble ratio (Beranek) from the reverberation time + `param` of the octave bands: BR = (T125 + T250) / (T500 + T1000), + TR = (T2000 + T4000) / (T500 + T1000). nan when a band is missing. + """ + bands, rt = np.asarray(self.bands), self._params[param] + + def get(f): # band within 6 % of the nominal f + near = np.flatnonzero(np.isclose(bands, f, rtol=0.06)) + return rt[near[0]] if near.size else np.nan + mid = get(500) + get(1000) + return (get(125) + get(250)) / mid, (get(2000) + get(4000)) / mid + def uncertainty(self, param: str = 'T30', numPositions: int = 1, numDecays: int = 10) -> np.ndarray: """ @@ -1523,6 +1552,7 @@ def plot_lundeby(): freqRange=[kwargs['minFreq'], kwargs['maxFreq']])[:, 1] listEDC = [] + listLundeby = [] for ch in range(hSignal.numChannels): signal = hSignal[ch] band = bands[ch] @@ -1538,11 +1568,12 @@ def plot_lundeby(): suppressWarnings=suppressWarnings, edcMethod=edcMethod) listEDC.append((energyDecay, energyVector)) + listLundeby.append(lundebyParams) if plotLundebyResults and lundebyParams is not None: # plot_lundeby(band, timeVector, timeSignal, samplingRate, # lundebyParams) plot_lundeby() - return listEDC, hSignal + return listEDC, hSignal, listLundeby # @njit def reverb_time_regression(energyDecay, energyVector, upperLim, lowerLim): diff --git a/tests/test_room_edc.py b/tests/test_room_edc.py index 560565f..970ded3 100644 --- a/tests/test_room_edc.py +++ b/tests/test_room_edc.py @@ -14,6 +14,7 @@ os.environ.setdefault('MPLBACKEND', 'Agg') import numpy as np +import scipy.signal as ss import pytta @@ -79,3 +80,44 @@ def test_uncertainty_follows_iso_3382_1(self): if __name__ == '__main__': unittest.main() + + +class TestRoomExtras(unittest.TestCase): + + def test_psnr_follows_the_noise_floor(self): + # The same decay over background noise 45 and 55 dB below its start: + # the peak-to-noise ratio of each band grows by 10 dB. Without the + # Lundeby estimate (noCut) there is no noise level, so no PSNR. + rng = np.random.default_rng(0) + t = np.arange(int(2.5 * FS)) / FS + decay = rng.standard_normal(t.size) * 10 ** (-3 * t / T) + noise = rng.standard_normal(t.size) + psnr = [] + for floor in (45, 55): + ir = pytta.SignalObj(decay + 10 ** (-floor / 20) * noise, 'time', + FS, signalType='energy') + psnr.append(analyse(ir).PSNR) + print(f'\n PSNR difference by band: {np.round(psnr[1] - psnr[0], 2)}') + self.assertTrue(np.allclose(psnr[1] - psnr[0], 10, atol=0.5)) + self.assertTrue(np.all(np.isnan(analyse(ir, edcMethod='noCut').PSNR))) + + def test_bass_and_treble_ratios(self): + # Octave bands decaying with T of 2.0, 1.8, 1.5, 1.5, 1.2 and 1.0 s + # from 125 Hz to 4 kHz: BR = (2.0 + 1.8) / 3.0, TR = (1.2 + 1.0) / 3.0. + rng = np.random.default_rng(1) + t = np.arange(4 * FS) / FS + h = np.zeros(t.size) + for fc, Tb in zip((125, 250, 500, 1000, 2000, 4000), + (2.0, 1.8, 1.5, 1.5, 1.2, 1.0)): + sos = ss.butter(8, [fc / 2 ** 0.5, fc * 2 ** 0.5], 'bandpass', + fs=FS, output='sos') + h += ss.sosfilt(sos, rng.standard_normal(t.size)) \ + * 10 ** (-3 * t / Tb) + ra = pytta.RoomAnalysis(pytta.SignalObj(h, 'time', FS, + signalType='energy'), + nthOct=1, minFreq=125., maxFreq=4000.) + BR, TR = ra.timbre_ratios('T30') + print(f'\n BR {BR:.3f} (expected {3.8 / 3:.3f}), TR {TR:.3f} ' + f'(expected {2.2 / 3:.3f})') + self.assertAlmostEqual(BR / (3.8 / 3), 1, delta=0.05) + self.assertAlmostEqual(TR / (2.2 / 3), 1, delta=0.05)