"""Room acoustic parameters
"""
import re
import numpy as np
[docs]
def reverberation_time_linear_regression(
energy_decay_curve, T='T20', return_intercept=False):
"""Estimate the reverberation time from a given energy decay curve.
The linear regression is performed using least squares error minimization
according to the ISO standard 3382 [#]_.
Parameters
----------
energy_decay_curve : ndarray, double
Energy decay curve. The time needs to be the arrays last dimension.
times : ndarray, double
Time vector corresponding to each sample of the EDC.
T : 'T20', 'T30', 'T40', 'T50', 'T60', 'EDT', 'LDT'
Decay interval to be used for the reverberation time extrapolation. EDT
corresponds to the early decay time extrapolated from the interval
[0, -10] dB, LDT corresponds to the late decay time extrapolated from
the interval [-25, -35] dB.
normalize : bool, True
Normalize the EDC to the steady state energy level
Returns
-------
reverberation_time : double
The reverberation time
References
----------
.. [#] ISO 3382, Acoustics - Measurement of the reverberation time of
rooms with reference to other acoustical parameters.
Examples
--------
Estimate the reverberation time from an energy decay curve.
>>> import numpy as np
>>> import pyfar as pf
>>> import pyrato as ra
>>> from pyrato.analytic import rectangular_room_rigid_walls
...
>>> L = np.array([8, 5, 3])/10
>>> source_pos = np.array([5, 3, 1.2])/10
>>> receiver_pos = np.array([1, 1, 1.2])/10
>>> rir, _ = rectangular_room_rigid_walls(
... L, source_pos, receiver_pos,
... reverberation_time=1, max_freq=1.5e3, n_samples=2**12,
... speed_of_sound=343.9, samplingrate=3e3)
>>> rir = rir/rir.time.max()
...
>>> awgn = pf.signals.noise(
... rir.n_samples, rms=10**(-50/20),
... sampling_rate=rir.sampling_rate)
>>> rir = rir + awgn
...
>>> edc = ra.energy_decay_curve_chu_lundeby(rir)
>>> t_20 = ra.reverberation_time_linear_regression(edc, 'T20')
>>> t_20
... array([0.99526253])
"""
intervals = [20, 30, 40, 50, 60]
if T == 'EDT':
upper = -0.1
lower = -10.1
elif T == 'LDT':
upper = -25.
lower = -35.
else:
try:
(int(re.findall(r'\d+', T)[0]) in intervals)
except IndexError:
raise ValueError(
"{} is not a valid interval for the regression.".format(T))
upper = -5
lower = -np.double(re.findall(r'\d+', T)) + upper
edcs_db = 10*np.log10(np.abs(energy_decay_curve.time))
times = energy_decay_curve.times
reverberation_times = np.zeros(energy_decay_curve.cshape, dtype=float)
intercepts = np.zeros(energy_decay_curve.cshape, dtype=float)
for ch in np.ndindex(energy_decay_curve.cshape):
edc_db = edcs_db[ch]
idx_upper = np.nanargmin(np.abs(upper - edc_db))
idx_lower = np.nanargmin(np.abs(lower - edc_db))
A = np.vstack(
[times[idx_upper:idx_lower], np.ones(idx_lower - idx_upper)]).T
gradient, const = np.linalg.lstsq(
A, edc_db[..., idx_upper:idx_lower], rcond=None)[0]
reverberation_times[ch] = -60 / gradient
intercepts[ch] = 10**(const/10)
if return_intercept is True:
return reverberation_times, intercepts
else:
return reverberation_times
