diff --git a/pytta/classes/analysis.py b/pytta/classes/analysis.py index 7f15924..33326fb 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,9 @@ 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, listLundeby = cumulative_integration( + ir, bypassLundeby, plotLundeby, suppressWarnings, + edcMethod=edcMethod, **kwargs) params = {} params['rms'] = fhSignal.rms() params['SPL'] = fhSignal.spl() @@ -1005,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 @@ -1063,6 +1076,53 @@ 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: + """ + 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 +1461,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 +1496,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 +1522,7 @@ def cumulative_integration(inputSignal, bypassLundeby, plotLundebyResults, suppressWarnings=True, + edcMethod='cutWithCorrection', **kwargs): """Cumulative integration with proper corrections.""" @@ -1473,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] @@ -1485,13 +1565,15 @@ 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. + 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 new file mode 100644 index 0000000..970ded3 --- /dev/null +++ b/tests/test_room_edc.py @@ -0,0 +1,123 @@ +""" +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 scipy.signal as ss + +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() + + +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)