Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 10 additions & 7 deletions examples/energy_decay_curves_and_reverberation_time.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -190,10 +190,13 @@
" The room impulse response with dimension [..., n_samples]\n",
"sampling_rate: integer\n",
" The sampling rate of the room impulse response.\n",
"freq: integer OR string\n",
" The frequency band. If set to 'broadband',\n",
" the time window of the Lundeby-algorithm will not be set in dependence\n",
" of frequency.\n",
"smoothing_parameter : int or array_like of int or {'broadband'}\n",
" Used to determine the smoothing time window in the Lundeby\n",
" algorithm. It should represent the center frequency (in Hz) of the\n",
" frequency band(s) in which the RIR data was computed.\n",
" If set to 'broadband', the smoothing time window will not be set\n",
" in dependence of frequency and a fixed time window of 30 ms\n",
" is used.\n",
"noise_level: ndarray, double OR string\n",
" If not specified, the noise level is calculated based on the last 10\n",
" percent of the RIR. Otherwise specify manually for each channel\n",
Expand Down Expand Up @@ -276,7 +279,7 @@
"source": [
"intersection_time, late_reveberation_time, noise_level_lundeby = \\\n",
" ra.edc.intersection_time_lundeby(\n",
" rir_1_noise, freq='broadband', is_energy=False,\n",
" rir_1_noise, smoothing_parameter='broadband', is_energy=False,\n",
" time_shift=True, channel_independent=False, plot=True)\n",
"\n",
"output_string = \"The estimated intersection time is {}s, the late reverberation time is {}s and the estimated noise is {}dB.\".format(intersection_time, late_reveberation_time, 10*np.log10(np.abs(noise_level_lundeby)))\n",
Expand All @@ -297,7 +300,7 @@
"outputs": [],
"source": [
"edc_truncation = ra.edc.energy_decay_curve_truncation(\n",
" rir_1_noise, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_1_noise, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True, plot=False)\n",
"\n",
"plt.figure(figsize=(10, 4))\n",
Expand Down Expand Up @@ -422,7 +425,7 @@
"outputs": [],
"source": [
"edc_chu_lundeby = ra.edc.energy_decay_curve_chu_lundeby(\n",
" rir_1_noise, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_1_noise, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True, plot=False)\n",
"\n",
"plt.figure(figsize=(10, 4))\n",
Expand Down
35 changes: 24 additions & 11 deletions pyrato/dsp.py
Original file line number Diff line number Diff line change
Expand Up @@ -417,7 +417,7 @@ def _smooth_rir(
The room impulse response with dimension ``(..., n_samples)``.
sampling_rate: integer
Defines the sampling rate of the room impulse response.
smooth_block_length : double
smooth_block_length : double or array-like of double
Defines the block-length of the smoothing algorithm in seconds.

Returns
Expand All @@ -430,25 +430,38 @@ def _smooth_rir(
The time vector fitting the original data.

"""
data = np.atleast_2d(data)
n_samples = data.shape[-1]
data = data.reshape(-1, n_samples)
smooth_block_length = (np.atleast_1d(smooth_block_length)).flatten()
n_channels = len(smooth_block_length)
n_samples_nan = np.count_nonzero(np.isnan(data), axis=-1)

n_samples_per_block = int(np.round(smooth_block_length * sampling_rate, 0))
n_samples_per_block = (np.round(
smooth_block_length * sampling_rate, 0)).astype(int)
n_blocks = np.asarray(
np.floor((n_samples-n_samples_nan)/n_samples_per_block),
dtype=int)

n_blocks_min = int(np.min(n_blocks))
n_samples_actual = int(n_blocks_min*n_samples_per_block)
reshaped_array = np.reshape(
data[..., :n_samples_actual],
(-1, n_blocks_min, n_samples_per_block))
time_window_data = np.mean(reshaped_array, axis=-1)

n_blocks_min = (np.min(n_blocks)).astype(int)
n_samples_actual = (n_blocks_min*n_samples_per_block).astype(int)
time_window_data = np.empty((n_channels, n_blocks_min))
# Use average time instances corresponding to the average energy level
# instead of time for the first sample of the block
time_vector_window = \
if n_channels > 1:
for ch in range(n_channels):
reshaped_array = np.reshape(
data[ch, :n_samples_actual[ch]],
(n_blocks_min, n_samples_per_block[ch]))
time_window_data[ch, :] = np.mean(reshaped_array, axis=-1)
time_vector_window = (
(0.5+np.arange(0, n_blocks_min)).reshape(1, -1) * (
n_samples_per_block/sampling_rate).reshape(-1, 1))
else:
reshaped_array = np.reshape(
Comment on lines +450 to +460
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not entirely satisfied with this cumbersome if statement, so I would appreciate a suggestion for improvement

data[..., :n_samples_actual[0]],
(-1, n_blocks_min, n_samples_per_block[0]))
time_window_data = np.mean(reshaped_array, axis=-1)
time_vector_window = \
((0.5+np.arange(0, n_blocks_min)) * n_samples_per_block/sampling_rate)

# Use the time corresponding to the sampling of the original data
Expand Down
87 changes: 56 additions & 31 deletions pyrato/edc.py
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,7 @@ def _schroeder_integration(impulse_response, is_energy=False):

def energy_decay_curve_truncation(
data,
freq='broadband',
smoothing_parameter='broadband',
noise_level='auto',
is_energy=False,
time_shift=True,
Expand All @@ -202,10 +202,13 @@ def energy_decay_curve_truncation(
----------
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.
smoothing_parameter : int or array_like of int or {'broadband'}
Used to determine the smoothing time window in the Lundeby
algorithm. It should represent the center frequency (in Hz) of the
frequency band(s) in which the RIR data was computed.
If set to 'broadband', the smoothing time window will not be set
in dependence of frequency and a fixed time window of 30 ms
is used.
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
Expand Down Expand Up @@ -284,7 +287,7 @@ def energy_decay_curve_truncation(

intersection_time = intersection_time_lundeby(
energy_data,
freq=freq,
smoothing_parameter=smoothing_parameter,
initial_noise_power=noise_level,
is_energy=True,
time_shift=False,
Expand Down Expand Up @@ -335,7 +338,7 @@ def energy_decay_curve_truncation(

def energy_decay_curve_lundeby(
data,
freq='broadband',
smoothing_parameter='broadband',
noise_level='auto',
is_energy=False,
time_shift=True,
Expand All @@ -357,10 +360,13 @@ def energy_decay_curve_lundeby(
----------
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.
smoothing_parameter : int or array_like of int or {'broadband'}
Used to determine the smoothing time window in the Lundeby
algorithm. It should represent the center frequency (in Hz) of the
frequency band(s) in which the RIR data was computed.
If set to 'broadband', the smoothing time window will not be set
in dependence of frequency and a fixed time window of 30 ms
is used.
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
Expand Down Expand Up @@ -435,7 +441,7 @@ def energy_decay_curve_lundeby(
intersection_time, late_reverberation_time, noise_estimation = \
intersection_time_lundeby(
energy_data,
freq=freq,
smoothing_parameter=smoothing_parameter,
initial_noise_power=noise_level,
is_energy=True,
time_shift=False,
Expand Down Expand Up @@ -620,7 +626,7 @@ def energy_decay_curve_chu(

def energy_decay_curve_chu_lundeby(
data,
freq='broadband',
smoothing_parameter='broadband',
noise_level='auto',
is_energy=False,
time_shift=True,
Expand All @@ -642,10 +648,13 @@ def energy_decay_curve_chu_lundeby(
----------
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.
smoothing_parameter : int or array_like of int or {'broadband'}
Used to determine the smoothing time window in the Lundeby
algorithm. It should represent the center frequency (in Hz) of the
frequency band(s) in which the RIR data was computed.
If set to 'broadband', the smoothing time window will not be set
in dependence of frequency and a fixed time window of 30 ms
is used.
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
Expand Down Expand Up @@ -730,7 +739,7 @@ def energy_decay_curve_chu_lundeby(
intersection_time, late_reverberation_time, noise_level = \
intersection_time_lundeby(
energy_data,
freq=freq,
smoothing_parameter=smoothing_parameter,
initial_noise_power=noise_level,
is_energy=True,
time_shift=False,
Expand Down Expand Up @@ -786,7 +795,7 @@ def energy_decay_curve_chu_lundeby(

def intersection_time_lundeby(
data,
freq='broadband',
smoothing_parameter='broadband',
initial_noise_power='auto',
is_energy=False,
time_shift=False,
Expand All @@ -803,10 +812,13 @@ def intersection_time_lundeby(
----------
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.
smoothing_parameter : int or array_like of int or {'broadband'}
Used to determine the smoothing time window in the Lundeby
algorithm. It should represent the center frequency (in Hz) of the
frequency band(s) in which the RIR data was computed.
If set to 'broadband', the smoothing time window will not be set
in dependence of frequency and a fixed time window of 30 ms
is used.
initial_noise_power: ndarray, double OR string
If ``'auto'``, the noise level is calculated based on the last 10
percent of the RIR. Otherwise specify manually for each channel
Expand Down Expand Up @@ -900,11 +912,24 @@ def intersection_time_lundeby(
sampling_rate = np.round(1/np.diff(data.times).mean(), decimals=4)
energy_data = energy_data.time

if freq == "broadband":
# number of frequency bands given by first channel axis
n_bands = np.prod(data.cshape)
if smoothing_parameter == "broadband":
# broadband: use 30 ms windows sizes
freq_dependent_window_time = 0.03
freq_dependent_window_time = np.atleast_1d([0.03] * n_bands)
elif isinstance(smoothing_parameter, (int, list, tuple, np.ndarray)):
smoothing_parameter = np.atleast_1d(smoothing_parameter)
if smoothing_parameter.size == 1:
smoothing_parameter = np.tile(smoothing_parameter, n_bands)
elif smoothing_parameter.size != n_bands:
raise ValueError(
"The size of smoothing_parameter must match the number of "
"frequency bands.")
freq_dependent_window_time = (800 / smoothing_parameter + 10) / 1000
else:
freq_dependent_window_time = (800/freq+10) / 1000
raise TypeError(
"smoothing_parameter must be an int or array_like of int "
"or {'broadband'}")

# (1) SMOOTH
time_window_data, time_vector_window, time_vector = dsp._smooth_rir(
Expand All @@ -923,12 +948,12 @@ def intersection_time_lundeby(
noise_peak_level = np.zeros(data.cshape, data.time.dtype)

for ch in np.ndindex(data.cshape):

output = _intersection_time_lundby(
time_window_data[ch], noise_estimation[ch], energy_data[ch],
time_vector_window, dB_above_noise, n_intervals_per_10dB,
use_dyn_range_for_regression, sampling_rate, ch, failure_policy)

time_window_data[ch], noise_estimation[ch], energy_data[ch],
np.squeeze(np.atleast_2d(time_vector_window)[ch, :]),
dB_above_noise, n_intervals_per_10dB,
use_dyn_range_for_regression, sampling_rate,
ch, failure_policy)
if output is None:
reverberation_time[ch] = np.nan
noise_level[ch] = np.nan
Expand Down
16 changes: 8 additions & 8 deletions tests/test_data/generate_test_data.py
Original file line number Diff line number Diff line change
Expand Up @@ -100,24 +100,24 @@
substracted_2D = pyrato.edc.subtract_noise_from_squared_rir(rir_array**2)

edc_truncation_1D = pyrato.energy_decay_curve_truncation(
rir_array[0], freq='broadband', is_energy=False, time_shift=True,
rir_array[0], smoothing_parameter='broadband', is_energy=False, time_shift=True,
channel_independent=False, normalize=True)
edc_truncation_2D = pyrato.energy_decay_curve_truncation(
rir_array, freq='broadband', is_energy=False, time_shift=True,
rir_array, smoothing_parameter='broadband', is_energy=False, time_shift=True,
channel_independent=False, normalize=True)

edc_lundeby_1D = pyrato.energy_decay_curve_lundeby(
rir_array[0], freq='broadband', is_energy=False, time_shift=True,
rir_array[0], smoothing_parameter='broadband', is_energy=False, time_shift=True,
channel_independent=False, normalize=True, plot=False)
edc_lundeby_2D = pyrato.energy_decay_curve_lundeby(
rir_array, freq='broadband', is_energy=False, time_shift=True,
rir_array, smoothing_parameter='broadband', is_energy=False, time_shift=True,
channel_independent=False, normalize=True, plot=False)

edc_lundeby_chu_1D = pyrato.energy_decay_curve_chu_lundeby(
rir_array[0], freq='broadband', is_energy=False, time_shift=True,
rir_array[0], smoothing_parameter='broadband', is_energy=False, time_shift=True,
channel_independent=False, normalize=True, plot=False)
edc_lundeby_chu_2D = pyrato.energy_decay_curve_chu_lundeby(
rir_array, freq='broadband', is_energy=False, time_shift=True,
rir_array, smoothing_parameter='broadband', is_energy=False, time_shift=True,
channel_independent=False, normalize=True, plot=False)

edc_chu_1D = pyrato.energy_decay_curve_chu(
Expand All @@ -128,10 +128,10 @@
channel_independent=False, normalize=True, plot=False)

intersection_time_1D = pyrato.intersection_time_lundeby(
rir_array[0], freq='broadband', is_energy=False, time_shift=True,
rir_array[0], smoothing_parameter='broadband', is_energy=False, time_shift=True,
channel_independent=False, plot=False)
intersection_time_2D = pyrato.intersection_time_lundeby(
rir_array, freq='broadband', is_energy=False, time_shift=True,
rir_array, smoothing_parameter='broadband', is_energy=False, time_shift=True,
channel_independent=False, plot=False)

# noise_energy_from_edc_1D = pyrato.edc.estimate_noise_energy_from_edc(
Expand Down
20 changes: 10 additions & 10 deletions tests/test_data/test_notebook.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -168,38 +168,38 @@
" substracted_2D = nh.subtract_noise_from_squared_rir(rir_array**2)\n",
"\n",
" edc_truncation_1D = nh.energy_decay_curve_truncation(\n",
" rir_array[0], sampling_rate, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_array[0], sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True)\n",
" edc_truncation_2D = nh.energy_decay_curve_truncation(\n",
" rir_array, sampling_rate, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_array, sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True)\n",
"\n",
" edc_lundeby_1D = nh.energy_decay_curve_lundeby(\n",
" rir_array[0], sampling_rate, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_array[0], sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True, plot=False)\n",
" edc_lundeby_2D = nh.energy_decay_curve_lundeby(\n",
" rir_array, sampling_rate, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_array, sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True, plot=False)\n",
"\n",
" edc_lundeby_chu_1D = nh.energy_decay_curve_chu_lundeby(\n",
" rir_array[0], sampling_rate, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_array[0], sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True, plot=False)\n",
" edc_lundeby_chu_2D = nh.energy_decay_curve_chu_lundeby(\n",
" rir_array, sampling_rate, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_array, sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True, plot=False)\n",
"\n",
" edc_chu_1D = nh.energy_decay_curve_chu(\n",
" rir_array[0], sampling_rate, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_array[0], sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True, plot=False)\n",
" edc_chu_2D = nh.energy_decay_curve_chu(\n",
" rir_array, sampling_rate, freq='broadband', is_energy=False, time_shift=True,\n",
" rir_array, sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=True,\n",
" channel_independent=False, normalize=True, plot=False)\n",
"\n",
" intersection_time_1D = nh.intersection_time_lundeby(\n",
" rir_array[0], sampling_rate, freq='broadband', is_energy=False, time_shift=False,\n",
" rir_array[0], sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=False,\n",
" channel_independent=False, plot=False)\n",
" intersection_time_2D = nh.intersection_time_lundeby(\n",
" rir_array, sampling_rate, freq='broadband', is_energy=False, time_shift=False,\n",
" rir_array, sampling_rate, smoothing_parameter='broadband', is_energy=False, time_shift=False,\n",
" channel_independent=False, plot=False)\n",
"\n",
" noise_energy_from_edc_1D = nh.estimate_noise_energy_from_edc(\n",
Expand Down
Loading
Loading