pyrato

Functions:

air_attenuation_coefficient(frequency[, ...])	Calculate the attenuation coefficient m for the absorption caused
energy_decay_curve_analytic(surfaces, ...[, ...])	Calculate the energy decay curve analytically by using Eyring's or Sabine's equation [1].
energy_decay_curve_chu(data[, noise_level, ...])	Implementation of the "subtraction of noise"-method after Chu [#] The noise level is estimated and subtracted from the impulse response before backward integration.
energy_decay_curve_chu_lundeby(data[, freq, ...])	This function combines Chu's and Lundeby's methods: The estimated noise level is subtracted before backward integration, the impulse response is truncated at the intersection time, and the correction for the truncation is applied [4], [5], [6]
energy_decay_curve_lundeby(data[, freq, ...])	Lundeby et al. [2] proposed a correction term to prevent the truncation error.
energy_decay_curve_truncation(data[, freq, ...])	This function truncates a given room impulse response by the intersection time after Lundeby and calculates the energy decay curve.
estimate_noise_energy(data[, interval, ...])	This function estimates the noise energy level of a given room impulse response.
find_impulse_response_maximum(impulse_response)	Find the maximum of an impulse response as argmax(h(t)).
find_impulse_response_start(impulse_response)	Find the start sample of an impulse response.
intersection_time_lundeby(data[, freq, ...])	Calculate the intersection time between impulse response and noise.
preprocess_rir(data[, is_energy, shift, ...])	Preprocess the room impulse response for further processing:
reverberation_time_energy_decay_curve(...[, T])	Estimate the reverberation time from a given energy decay curve.
reverberation_time_linear_regression(...[, ...])	Estimate the reverberation time from a given energy decay curve.
schroeder_integration(room_impulse_response)	Calculate the Schroeder integral of a room impulse response [3].
time_shift(signal, shift[, circular_shift, unit])	Apply a time-shift to a signal.

pyrato.air_attenuation_coefficient(frequency, temperature=20, humidity=50, atmospheric_pressure=101325)

Calculate the attenuation coefficient m for the absorption caused

by friction with the surrounding air.

Parameters:

• frequency (double) – The frequency for which the attenuation coefficient is calculated. When processing in fractional octave bands use the center frequency.

• temperature (double) – Temperature in degrees Celsius.

• humidity (double) – Humidity in percent.

• atmospheric_pressure (double) – Atmospheric pressure.

Returns:

attenuation_coefficient – The resulting attenuation coefficient.

Return type:

double

pyrato.energy_decay_curve_analytic(surfaces, alphas, volume, times, source=None, receiver=None, method='eyring', c=343.4, frequency=None, air_absorption=True)

Calculate the energy decay curve analytically by using Eyring’s or Sabine’s equation [7].

Parameters:

• surfaces (ndarray, double) – Surface areas of all surfaces in the room

• alphas (ndarray, double) – Absorption coefficients corresponding to each surface

• volume (double) – Room volume

• times (ndarray, double) – Time vector for which the decay curve is calculated

• source (Coordinates) – Coordinate object with the source coordinates

• receiver (Coordinates) – Coordinate object with the receiver coordinates

• method ('eyring', 'sabine') – Use either Eyring’s or Sabine’s equation

• c (double) – Speed of sound

• frequency (double, optional) – Center frequency of the respective octave band. This is only used for the air absorption calculation.

Returns:

energy_decay_curve – The energy decay curve

Return type:

ndarray, double

References

[1]

8. Kuttruff, Room acoustics, 4th Ed. Taylor & Francis, 2009.

pyrato.energy_decay_curve_chu(data, noise_level='auto', is_energy=False, time_shift=True, channel_independent=False, normalize=True, threshold=10, plot=False)

Implementation of the “subtraction of noise”-method after Chu [#] The noise level is estimated and subtracted from the impulse response before backward integration.

Parameters:

• data (ndarray, double) – The room impulse response with dimension […, n_samples]

• noise_level (ndarray, double OR string) – If not specified, the noise level is calculated based on the last 10 percent of the RIR. Otherwise specify manually for each channel as array.

• is_energy (boolean) – Defines, if the data is already squared.

• time_shift (boolean) – Defines, if the silence at beginning of the RIR should be removed.

• channel_independent (boolean) – Defines, if the time shift and normalizsation is done channel-independently or not.

• normalize (boolean) – Defines, if the energy decay curve should be normalized in the end or not.

• threshold (float, None) – Defines a peak-signal-to-noise ratio based threshold in dB for final truncation of the EDC. Values below the sum of the threshold level and the peak-signal-to-noise ratio in dB are discarded. The default is 10 dB. If None, the decay curve will not be truncated further.

• plot (Boolean) – Specifies, whether the results should be visualized or not.

Returns:

energy_decay_curve – Returns the noise handeled edc.

Return type:

ndarray, double

References

[2]

W. T. Chu. “Comparison of reverberation measurements using Schroeder’s impulse method and decay-curve averaging method”. In: Journal of the Acoustical Society of America 63.5 (1978), pp. 1444–1450.

Examples

>>> 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=1e3, n_samples=2**16,
...     speed_of_sound=343.9)
...
>>> awgn = pf.signals.noise(
...     rir.n_samples, rms=rir.time.max()*10**(-40/20),
...     sampling_rate=rir.sampling_rate)
>>> rir = rir + awgn
>>> edc = ra.energy_decay_curve_chu(rir)
...
>>> pf.plot.time(rir/np.abs(rir.time).max(), dB=True, label='RIR')
>>> ax = pf.plot.time(
...     edc/edc.time[..., 0], dB=True, log_prefix=10, label='EDC')
>>> ax.set_ylim(-65, 5)
>>> ax.legend()

(Source code, png, hires.png, pdf)

pyrato.energy_decay_curve_chu_lundeby(data, freq='broadband', noise_level='auto', is_energy=False, time_shift=True, channel_independent=False, normalize=True, plot=False)

This function combines Chu’s and Lundeby’s methods: The estimated noise level is subtracted before backward integration, the impulse response is truncated at the intersection time, and the correction for the truncation is applied [4], [5], [6]

Parameters:

• data (pyfar.Signal) – The room impulse response.

• freq (integer OR string) – The frequency band. If set to ‘broadband’, the time window of the Lundeby-algorithm will not be set in dependence of frequency.

• noise_level (ndarray, double OR string) – If not specified, the noise level is calculated based on the last 10 percent of the RIR. Otherwise specify manually for each channel as array.

• is_energy (boolean) – Defines, if the data is already squared.

• time_shift (boolean) – Defines, if the silence at beginning of the RIR should be removed.

• channel_independent (boolean) – Defines, if the time shift and normalizsation is done channel-independently or not.

• normalize (boolean) – Defines, if the energy decay curve should be normalized in the end or not.

• plot (Boolean) – Specifies, whether the results should be visualized or not.

Returns:

Returns the noise handeled edc.

Return type:

pyfar.TimeData

References

[4] (1,2)

Lundeby, Virgran, Bietz and Vorlaender - Uncertainties of Measurements in Room Acoustics - ACUSTICA Vol. 81 (1995)

[5] (1,2)

W. T. Chu. “Comparison of reverberation measurements using Schroeder’s impulse method and decay-curve averaging method”. In: Journal of the Acoustical Society of America 63.5 (1978), pp. 1444–1450.

[6] (1,2)

M. Guski, “Influences of external error sources on measurements of room acoustic parameters,” 2015.

Examples

Calculate and plot the EDC using a combination of Chu’s and Lundeby’s methods.

>>> 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=1e3, n_samples=2**16,
...     speed_of_sound=343.9)
>>> rir = rir/rir.time.max()
...
>>> awgn = pf.signals.noise(
...     rir.n_samples, rms=rir.time.max()*10**(-50/20),
...     sampling_rate=rir.sampling_rate)
>>> rir = rir + awgn
>>> edc = ra.energy_decay_curve_chu_lundeby(rir)
...
>>> ax = pf.plot.time(rir, dB=True, label='RIR')
>>> pf.plot.time(edc, dB=True, log_prefix=10, label='EDC')
>>> ax.set_ylim(-65, 5)
>>> ax.legend()

(Source code, png, hires.png, pdf)

pyrato.energy_decay_curve_lundeby(data, freq='broadband', noise_level='auto', is_energy=False, time_shift=True, channel_independent=False, normalize=True, plot=False)

Lundeby et al. [8] proposed a correction term to prevent the truncation error. The missing signal energy from truncation time to infinity is estimated and added to the truncated integral.

Parameters:

• data (pyfar.Signal) – The room impulse response.

• freq (integer OR string) – The frequency band. If set to ‘broadband’, the time window of the Lundeby-algorithm will not be set in dependence of frequency.

• noise_level (ndarray, double OR string) – If not specified, the noise level is calculated based on the last 10 percent of the RIR. Otherwise specify manually for each channel as array.

• is_energy (boolean) – Defines, if the data is already squared.

• time_shift (boolean) – Defines, if the silence at beginning of the RIR should be removed.

• channel_independent (boolean) – Defines, if the time shift and normalizsation is done channel-independently or not.

• normalize (boolean) – Defines, if the energy decay curve should be normalized in the end or not.

• plot (Boolean) – Specifies, whether the results should be visualized or not.

Returns:

Returns the noise handeled edc.

Return type:

pyfar.TimeData

References

[3]

Lundeby, Virgran, Bietz and Vorlaender - Uncertainties of Measurements in Room Acoustics - ACUSTICA Vol. 81 (1995)

Examples

Plot the RIR and the EDC calculated after Lundeby.

>>> 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=1e3, n_samples=2**16,
...     speed_of_sound=343.9)
>>> rir = rir/rir.time.max()
...
>>> awgn = pf.signals.noise(
...     rir.n_samples, rms=rir.time.max()*10**(-50/20),
...     sampling_rate=rir.sampling_rate)
>>> rir = rir + awgn
>>> edc = ra.energy_decay_curve_lundeby(rir)
...
>>> ax = pf.plot.time(rir, dB=True, label='RIR')
>>> pf.plot.time(edc, dB=True, log_prefix=10, label='EDC')
>>> ax.set_ylim(-65, 5)
>>> ax.legend()

(Source code, png, hires.png, pdf)

pyrato.energy_decay_curve_truncation(data, freq='broadband', noise_level='auto', is_energy=False, time_shift=True, channel_independent=False, normalize=True, threshold=15, plot=False)

This function truncates a given room impulse response by the intersection time after Lundeby and calculates the energy decay curve.

Parameters:

• data (pyfar.Signal) – The room impulse response.

• freq (integer OR string) – The frequency band. If set to ‘broadband’, the time window of the Lundeby-algorithm will not be set in dependence of frequency.

• noise_level (ndarray, double OR string) – If not specified, the noise level is calculated based on the last 10 percent of the RIR. Otherwise specify manually for each channel as array.

• is_energy (boolean) – Defines, if the data is already squared.

• time_shift (boolean) – Defines, if the silence at beginning of the RIR should be removed.

• channel_independent (boolean) – Defines, if the time shift and normalization is done channel-independently or not.

• normalize (boolean) – Defines, if the energy decay curve should be normalized in the end or not.

• threshold (float) – Defines a peak-signal-to-noise ratio based threshold in dB for final truncation of the EDC. Values below the sum of the threshold level and the peak-signal-to-noise ratio in dB are discarded. The default is 15 dB, which is in correspondence with ISO 3382-1:2009 [9].

• plot (Boolean) – Specifies, whether the results should be visualized or not.

Returns:

Returns the noise compensated EDC.

Return type:

pyfar.TimeData

Examples

Plot the RIR and the EDC calculated truncating the integration at the intersection time.

>>> 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=1e3, n_samples=2**16,
...     speed_of_sound=343.9)
>>> rir = rir/rir.time.max()
...
>>> awgn = pf.signals.noise(
...     rir.n_samples, rms=rir.time.max()*10**(-50/20),
...     sampling_rate=rir.sampling_rate)
>>> rir = rir + awgn
>>> edc = ra.energy_decay_curve_truncation(rir)
...
>>> ax = pf.plot.time(rir, dB=True, label='RIR')
>>> pf.plot.time(edc, dB=True, log_prefix=10, label='EDC')
>>> ax.set_ylim(-65, 5)
>>> ax.legend()

(Source code, png, hires.png, pdf)

References

[7]

International Organization for Standardization, “EN ISO 3382-1:2009 Acoustics - Measurement of room acoustic parameters,” 2009.

pyrato.estimate_noise_energy(data, interval=[0.9, 1.0], is_energy=False)

This function estimates the noise energy level of a given room impulse response. The noise is assumed to be Gaussian.

Parameters:

• data (np.array) – The room impulse response with dimension […, n_samples]

• interval (tuple, float, [0.9, 1.]) – Defines the interval of the RIR to be evaluated for the estimation. The interval is relative to the length of the RIR [0 = 0%, 1=100%)]

• is_energy (Boolean) – Defines if the data is already squared.

Returns:

noise_energy – The energy of the background noise.

Return type:

float

pyrato.find_impulse_response_maximum(impulse_response, threshold=20, noise_energy='auto')

Find the maximum of an impulse response as argmax(h(t)). Performs an initial SNR check according to a defined threshold level in dB.

Parameters:

• impulse_response (ndarray, double) – The impulse response

• threshold (double, optional) – Threshold SNR value in dB

Returns:

max_sample – Sample at which the impulse response starts

Return type:

int

Note

The function tries to estimate the SNR in the IR based on the signal energy in the last 10 percent of the IR.

pyrato.find_impulse_response_start(impulse_response, threshold=20)

Find the start sample of an impulse response. The start sample is identified as the first sample which is below the threshold level relative to the maximum level of the impulse response. For room impulse responses, ISO 3382 [10] specifies a threshold of 20 dB. This function is primary intended to be used when processing room impulse responses.

Parameters:

• impulse_response (pyfar.Signal) – The impulse response

• threshold (float, optional) – The threshold level in dB, by default 20, which complies with ISO 3382.

Returns:

start_sample – Sample at which the impulse response starts

Return type:

numpy.ndarray, int

Notes

The function tries to estimate the PSNR in the IR based on the signal power in the last 10 percent of the IR. The automatic estimation may fail if the noise spectrum is not white or the impulse response contains non-linear distortions. If the PSNR is lower than the specified threshold, the function will issue a warning.

References

[8]

ISO 3382-1:2009-10, Acoustics - Measurement of the reverberation time of rooms with reference to other acoustical parameters. pp. 22

Examples

Create a band-limited impulse shifted by 0.5 samples and estimate the starting sample of the impulse and plot. .. plot:

>>> import pyfar as pf
>>> import numpy as np
>>> n_samples = 256
>>> delay_samples = n_samples // 2 + 1/2
>>> ir = pf.signals.impulse(n_samples)
>>> ir = pf.dsp.linear_phase(ir, delay_samples, unit='samples')
>>> start_samples = pf.dsp.find_impulse_response_start(ir)
>>> ax = pf.plot.time(ir, unit='ms', label='impulse response', dB=True)
>>> ax.axvline(
...     start_samples/ir.sampling_rate*1e3,
...     color='k', linestyle='-.', label='start sample')
>>> ax.axhline(
...     20*np.log10(np.max(np.abs(ir.time)))-20,
...     color='k', linestyle=':', label='threshold')
>>> ax.legend()

Create a train of weighted impulses with levels below and above the threshold, serving as a very abstract room impulse response. The starting sample is identified as the last sample below the threshold relative to the maximum of the impulse response. .. plot:

>>> import pyfar as pf
>>> import numpy as np
>>> n_samples = 64
>>> delays = np.array([14, 22, 26, 30, 33])
>>> amplitudes = np.array([-35, -22, -6, 0, -9], dtype=float)
>>> ir = pf.signals.impulse(n_samples, delays, 10**(amplitudes/20))
>>> ir.time = np.sum(ir.time, axis=0)
>>> start_sample_est = pf.dsp.find_impulse_response_start(
...     ir, threshold=20)
>>> ax = pf.plot.time(
...     ir, dB=True, unit='samples',
...     label=f'peak samples: {delays}')
>>> ax.axvline(
...     start_sample_est, linestyle='-.', color='k',
...     label=f'ir start sample: {start_sample_est}')
>>> ax.axhline(
...     20*np.log10(np.max(np.abs(ir.time)))-20,
...     color='k', linestyle=':', label='threshold')
>>> ax.legend()

pyrato.intersection_time_lundeby(data, freq='broadband', initial_noise_power='auto', is_energy=False, time_shift=False, channel_independent=False, plot=False)

Calculate the intersection time between impulse response and noise.

This function uses the algorithm after Lundeby et al. [11] to calculate the intersection time, lundeby reverberation time, and noise level estimation.

Parameters:

• data (pyfar.Signal) – The room impulse response

• freq (integer OR string) – The frequency band. If set to ‘broadband’, the time window of the Lundeby-algorithm will not be set in dependence of frequency.

• noise_level (ndarray, double OR string) – If not specified, the noise level is calculated based on the last 10 percent of the RIR. Otherwise specify manually for each channel as array.

• is_energy (boolean) – Defines, if the data is already squared.

• time_shift (boolean) – Defines, if the silence at beginning of the RIR should be removed.

• channel_independent (boolean) – Defines, if the time shift and normalizsation is done channel-independently or not.

• plot (Boolean) – Specifies, whether the results should be visualized or not.

Returns:

• intersection_time (ndarray, float) – Returns the Lundeby intersection time for each channel.

• reverberation_time (ndarray, float) – Returns the Lundeby reverberation time for each channel.

• noise_level (ndarray, float) – Returns the noise level estimation for each channel.

References

[9]

Lundeby, Virgran, Bietz and Vorlaender - Uncertainties of Measurements in Room Acoustics - ACUSTICA Vol. 81 (1995)

Examples

Estimate the intersection time and plot the RIR and the estimated noise power.

>>> 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=1e3, n_samples=2**16,
...     speed_of_sound=343.9)
>>> rir = rir/np.abs(rir.time).max()
...
>>> awgn = pf.signals.noise(
...     rir.n_samples, rms=rir.time.max()*10**(-40/20),
...     sampling_rate=rir.sampling_rate)
>>> rir = rir + awgn
>>> inter_time, _, noise_power = ra.intersection_time_lundeby(rir)
...
>>> ax = pf.plot.time(rir, dB=True, label='RIR')
>>> ax.axvline(inter_time, c='k', linestyle='--', label='$T_i$')
>>> ax.axhline(
...     10*np.log10(noise_power), c='k', linestyle=':', label='Noise')
>>> ax.set_ylim(-65, 5)
>>> ax.legend()

(Source code, png, hires.png, pdf)

pyrato.preprocess_rir(data, is_energy=False, shift=False, channel_independent=False)

Preprocess the room impulse response for further processing:

• Square data

• Shift the RIR to the first sample of the array, compensating for the delay of the time of arrival of the direct sound. The time shift is performed as a non-cyclic shift, adding numpy.nan values in the end of the RIR corresponding to the number of samples the data is shifted by.

• The time shift can be done channel-independent or not.

Parameters:

• data (ndarray, double) – The room impulse response with dimension […, n_samples]

• is_energy (boolean) – Defines, if the data is already squared.

• shift (boolean) – Defines, if the silence at beginning of the RIR should be removed.

• channel_independent (boolean) – Defines, if the time shift is done channel-independent or not.

Returns:

energy_data – The preprocessed RIR

Return type:

ndarray, double

pyrato.reverberation_time_energy_decay_curve(energy_decay_curve, T='T20')

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 [12].

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 – The reverberation time

Return type:

double

References

[10]

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_energy_decay_curve(edc, 'T20')
>>> t_20
...     array([0.99526253])

pyrato.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 – The reverberation time

Return type:

double

References

[11]

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])

pyrato.schroeder_integration(room_impulse_response, is_energy=False)

Calculate the Schroeder integral of a room impulse response [#]_. The result is the energy decay curve for the given room impulse response.

Parameters:

• room_impulse_response (pyfar.Signal) – Room impulse response as array

• is_energy (boolean, optional) – Whether the input represents energy data or sound pressure values.

Returns:

energy_decay_curve – The energy decay curve

Return type:

pyfar.TimeData

Note

This function does not apply any compensation of measurement noise and integrates the full length of the input signal. It should only be used if no measurement noise or artifacts are present in the data.

References

[12]

M. R. Schroeder, “New Method of Measuring Reverberation Time,” The Journal of the Acoustical Society of America, vol. 37, no. 6, pp. 1187-1187, 1965.

Example

Calculate the Schroeder integral of a simulated RIR and plot.

>>> 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=1e3, n_samples=2**16,
...     speed_of_sound=343.9)
>>> edc = ra.schroeder_integration(rir)
>>> pf.plot.time(rir/np.abs(rir.time).max(), dB=True, label='RIR')
>>> ax = pf.plot.time(
...     edc/edc.time[..., 0], dB=True, log_prefix=10, label='EDC')
>>> ax.set_ylim(-65, 5)
>>> ax.legend()

(Source code, png, hires.png, pdf)

pyrato.time_shift(signal, shift, circular_shift=True, unit='samples')

Apply a time-shift to a signal.

By default, the shift is performed as a cyclic shift on the time axis, potentially resulting in non-causal signals for negative shift values. Use the option circular_shift=False to pad with nan values instead, note that in this case the return type will be a pyfar.TimeData.

Parameters:

• signal (Signal) – The signal to be shifted

• shift (int, float) – The time-shift value. A positive value will result in right shift on the time axis (delaying of the signal), whereas a negative value yields a left shift on the time axis (non-causal shift to a earlier time). If a single value is given, the same time shift will be applied to each channel of the signal. Individual time shifts for each channel can be performed by passing an array matching the signals channel dimensions cshape.

• unit (str, optional) – Unit of the shift variable, this can be either 'samples' or 's' for seconds. By default 'samples' is used. Note that in the case of specifying the shift time in seconds, the value is rounded to the next integer sample value to perform the shift.

• circular_shift (bool, True) – Perform a circular or non-circular shift. If a non-circular shift is performed, the data will be padded with nan values at the respective beginning or ending of the data, corresponding to the number of samples the data is shifted. In this case, a pyfar.TimeData object is returned.

Returns:

The time-shifted signal. If a circular shift is performed, the return value will be a pyfar.Signal, in case of a non-circular shift, its type will be pyfar.TimeData.

Return type:

pyfar.Signal, pyfar.TimeData

Notes

This function is primarily intended to be used when processing room impulse responses. When circular_shift=True, the function input is passed into pyfar.dsp.time_shift.

Examples

Perform a circular shift a set of ideal impulses stored in three different channels and plot the resulting signals

>>> import pyfar as pf
>>> import pyrato as ra
>>> import matplotlib.pyplot as plt
>>> impulse = pf.signals.impulse(
...     32, amplitude=(1, 1.5, 1), delay=(14, 15, 16))
>>> shifted = ra.time_shift(impulse, [-2, 0, 2])
>>> pf.plot.use('light')
>>> _, axs = plt.subplots(2, 1)
>>> pf.plot.time(impulse, ax=axs[0])
>>> pf.plot.time(shifted, ax=axs[1])
>>> axs[0].set_title('Original signals')
>>> axs[1].set_title('Shifted signals')
>>> plt.tight_layout()

(Source code, png, hires.png, pdf)

Perform a non-circular shift a single impulse and plot the results.

>>> import pyfar as pf
>>> import pyrato as ra
>>> import matplotlib.pyplot as plt
>>> impulse = pf.signals.impulse(32, delay=15)
>>> shifted = ra.time_shift(impulse, -10, circular_shift=False)
>>> pf.plot.use('light')
>>> _, axs = plt.subplots(2, 1)
>>> pf.plot.time(impulse, ax=axs[0])
>>> pf.plot.time(shifted, ax=axs[1])
>>> axs[0].set_title('Original signal')
>>> axs[1].set_title('Shifted signal')
>>> plt.tight_layout()

(Source code, png, hires.png, pdf)