Open an interactive online version by clicking the badge

or download the notebook.

Room Acoustics: Measurement of Reverberation Time#

Marco Berzborn

Eindhoven University of Technology, Building Acoustics Group

Contact: m.berzborn@tue.nl

This assignment is part of the course “Architectural Acoustics” at the Eindhoven University of Technology. The goal of the assignment is to estimate the reverberation time of a room based on a measured room impulse response. The assignment covers the band-pass filtering of the room impulse response into octave bands, the calculation of the energy decay curve using Schroeder integration, and the estimation of the reverberation time using the T20 method. Finally, you will be asked to listen to auralizations of speech and music in the room and to evaluate listening experience in the room.

Duration: 90 Minutes

Requirements: Basic knowledge of room acoustics and reverberation time estimation. Familiarity with Python and the used libraries is helpful but not required.

References

1. Kuttruff, Room acoustics, Taylor & Francis.

Dependencies

If you are running this notebook locally or in Google Colab, make sure to install the required dependencies. You can do this by executing the following command in an arbitrary code cell:

!pip install numpy matplotlib pyfar ipympl pyrato

[ ]:

import pyfar as pf
import numpy as np
import matplotlib.pyplot as plt
import pyrato
from matplotlib.ticker import EngFormatter, NullFormatter
from pyfar.plot.ticker import FractionalOctaveLocator
from IPython.display import Audio as AudioPlayer
%matplotlib ipympl

sampling_rate = 48000

Load the data#

For this task you will use a measured room impulse response which is provided in the following code. The function downloads the data and returns it as a pyfar.Signal object, which encapsulate the data as well as relevant meta data. For more information on the pyfar.Signal class, please refer to the pyfar documentation.

[ ]:

room_impulse_response = pf.signals.files.room_impulse_response(sampling_rate)

You can visualize the room impulse response in the time domain by executing the following cell. In order to switch between the time and frequency domain, you can use the hotkeys Shift + t for time domain and Shift + f for frequency domain. For the hotkeys to work, make sure to click on the plot first.

[ ]:

plt.figure()
pf.plot.time(room_impulse_response, dB=True)

Fractional octave filtering#

Task: In order to analyze the reverberation time in different frequency bands, you are required to apply fractional octave band filters to the room impulse response. You can use the function pf.dsp.filter.fractional_octave_bands to apply the filters. The function takes the input signal, the number of fractions per octave, and the frequency range as input arguments. Use the frequency range already specified in the cell, i.e. frequency_range.

Make sure to name the result room_impulse_response_octave to ensure that the following cells work properly.

[ ]:

frequency_range = [125, 8e3]

# YOUR CODE HERE
raise NotImplementedError()

To check your solution, you can visualize the filtered room impulse response by executing the following cell. You can use the same hotkeys as before to switch between the time and frequency domain.

Hint: You can use the keys [ and ] to switch between the different channels in the plot. Alternatively, , and . work as well. To go back to showing all channels press Shift + a. Make sure to click on the plot first for the hotkeys to work.

[ ]:

bands = pf.dsp.filter.fractional_octave_frequencies(1, frequency_range=frequency_range)[0]

plt.figure()
pf.plot.time(room_impulse_response_octave, dB=True, label=[f"{band:.0f} Hz" for band in bands])
plt.legend()

Schroeder integration#

The code in the following cell calculates the energy decay curve for each frequency band using Schroeder integration. The function further ensures that measurement noise influences are removed from the energy decay curve. For the remainder of the assignment you can assume that the returned energy decay curve is free of noise influences and valid for the calculation of the reverberation time.

[ ]:

energy_decay_curve = pyrato.energy_decay_curve_lundeby(room_impulse_response_octave[:, 0], sampling_rate)

Visualize the energy decay curve by executing the following cell.

[ ]:

plt.figure()
pf.plot.time(energy_decay_curve, dB=True, log_prefix=10, label=[f"{band:.0f} Hz" for band in bands])
plt.legend()

Task: Your next task is to calculate the reverberation time for each frequency band using the energy decay curves obtained from the Schroeder integration. For this you can complete the function calculate_reverberation_time in the next cell. Assume that the energy decay curve passed into the function always only contains a single channel. The function should calculate the reverberation time using the T20 method, which which corresponds to the dynamic range of 20 dB, i.e. the range between -5 dB and -25 dB. The function should return the reverberation time in seconds. You can use the function linregress from the scipy.stats module to perform a linear regression on the energy decay curve corresponding range. You can verify the solution by executing the following cell, which calculates the reverberation time for each frequency band and prints the results.

[17]:

from scipy.stats import linregress

def calculate_reverberation_time(energy_decay_curve):
    """
    This method calculated the reverberation time, more specifically T20, from a given energy decay curve.
    The method uses linear regression to fit a line to the energy decay curve in the interval between -5 dB and -25 dB.

    Parameters
    ----------
    energy_decay_curve: pyfar.Signal
        The energy decay curve for which the reverberation time should be calculated.

    Returns
    -------
    reverberation_time: float
        The calculated reverberation time in seconds.

    """
    # The energy decay curve on logarithmic scale
    log_energy_decay_curve = 10 * np.log10(np.squeeze(energy_decay_curve.time))
    # The time values for which the energy decay curve is measured
    t = energy_decay_curve.times

    # YOUR CODE HERE
    raise NotImplementedError()
    return reverberation_time

[ ]:

reverberation_time = np.zeros(len(bands), dtype=float)
for bdx, band in enumerate(bands):
    reverberation_time[bdx] = calculate_reverberation_time(energy_decay_curve[bdx])
    print(f"{band:.0f} Hz: {reverberation_time[bdx]:.2f} seconds")

Finally, plot the reverberation time as a function of frequency by executing the following cell.

[ ]:

plt.figure()
ax = plt.axes()
ax.semilogx(bands, reverberation_time, marker='o')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Reverberation Time (s)')
ax.grid(True)
ax.xaxis.set_major_locator(FractionalOctaveLocator(1))
ax.xaxis.set_major_formatter(EngFormatter(unit='', sep=''))
ax.xaxis.set_minor_locator(FractionalOctaveLocator(3))
ax.xaxis.set_minor_formatter(NullFormatter())

Auralization#

Your final task is to auralize some audio signals using the room impulse response and listen to the results. For this, you can firstly download some audio signals using the pf.signals.files using the code provided in the following cell.

[ ]:

drums = pf.signals.files.drums(sampling_rate)
speech = pf.signals.files.speech(voice='male',sampling_rate=sampling_rate)
guitar = pf.signals.files.guitar(sampling_rate)

You can listen to the dry signals by executing the two following cells and pressing the respective play button.

[ ]:

AudioPlayer(np.squeeze(guitar.time), rate=sampling_rate)

[ ]:

AudioPlayer(np.squeeze(speech.time), rate=sampling_rate)

Task: Now use the function pf.dsp.convolve to convolve the audio signals with the room impulse response to create the auralizations for the two audio signals using the room impulse response. Hint: Make sure to name the resulting auralizations auralization_speech and auralization_music to ensure that the following cells work properly.

[ ]:

# YOUR CODE HERE
raise NotImplementedError()

You can again execute the following cells to listen to the auralized signals.

[ ]:

AudioPlayer(np.squeeze(auralized_guitar.time), rate=sampling_rate)

[ ]:

AudioPlayer(np.squeeze(auralized_speech.time), rate=sampling_rate)

Task: Listen to the auralized signals and compare them to the dry signals. Would you consider the room is well suitable for speech reproduction? Justify your answer. You can also consider the reverberation times previously calculated to support your answer. Give your answer in the following cell.

YOUR ANSWER HERE

License notice#

Watermark#

The following watermark might help others to install specific package versions that might be required to run the notebook. Please give at least the versions of Python, IPython, numpy , and scipy, major third party packagers (e.g., pytorch), and all used pyfar packages.

[ ]:

%load_ext watermark
%watermark -v -m -p numpy,scipy,pyfar,sofar,watermark