rt60

Estimate decay time of room impulse response

Since R2025a

Syntax

rtsummary = rt60(rir,fs)

rtsummary = rt60(rir,fs,Name=Value)

rt60(___)

Description

rtsummary = rt60(rir,fs) returns a table describing the amount of time for the reverberation energy of the room impulse response to decay 60 dB (RT-60). The function estimates the decay time in each octave band according to the ISO-3382 standard [2][3]. The returned table also contains metrics described in ISO-3382 such as the curvature, C50 (speech clarity), C80 (music clarity), D50 (definition), and Ts (center time).

rtsummary = rt60(rir,fs,Name=Value) specifies options using one or more name-value arguments.

rt60(___) without any output arguments plots the reverberation time metrics for each octave band. You can select the octave bands in the plot to see the impulse response power and energy decay curve in each band.

Examples

Measure RT-60 from Room Impulse Response

Read in a room impulse response for the interior of a church from a file.

[rir,fs] = audioread("ChurchImpulseResponse-16-44p1-mono-5secs.wav");

Call rt60 with the impulse response to get a table containing the RT-60 metrics for each octave center frequency.

metrics = rt60(rir,fs)

metrics=6×11 table
    CenterFrequencies     EDT       T20       T30       Topt              DegreeOfNonLinearity              Curvature      C50          C80        D50         Ts   
    _________________    ______    ______    ______    ______    _______________________________________    _________    ________    _________    ______    ________

           125           6.0856       NaN       NaN    6.2365    7.3997      2.438      2.438     1.7044         NaN      -6.0616      -4.8223    19.849     0.44521
           250           6.0207    6.4347    6.7053    6.4347    5.1429    0.58978     1.0775    0.58978      4.2059      -6.9081      -4.8076    16.929     0.40977
           500           5.3181    6.0276    6.0413    5.9859    3.0841    0.77807    0.49503     0.3772     0.22695      -2.9446       -2.438     33.67       0.319
          1000           4.7896    5.3202    5.5298    5.3276    5.4525    0.26611    0.89927     0.2395      3.9388     -0.66439    -0.039779    46.183     0.22816
          2000           4.0822    4.4473    4.5732    4.3312    14.253    0.64798    0.53592    0.30591      2.8315       1.6968       2.2301    59.645     0.15107
          4000           2.9998    3.3613    3.5392    3.2139    53.436     1.5334     1.6486    0.40026      5.2938       3.8169       4.4364    70.659    0.087573

Call rt60 with no output arguments to plot the RT-60 metrics for each octave band. You can select the different octave bands to view the energy decay curve and linear fits in each band.

rt60(rir,fs)

Specify Octave Filter Bank Parameters for Computing RT-60

Load a room impulse response and call rt60 to compute the RT-60 metrics for one third octave bands on the range from 100 to 16,000 Hz.

[rir,fs] = audioread("ChurchImpulseResponse-16-44p1-mono-5secs.wav");
rt60(rir,fs,Bandwidth="1/3 octave",FrequencyRange=[100,16000])

Call rt60 again on the same impulse response, this time modifying the octave filter bank parameters by using a higher filter order and filtering in the backward direction. See how the modifications affect the RT-60 metrics. The rt60 function throws a warning if the specified octave filter bank does not comply with the IEC 61260-1 standard.

rt60(rir,fs, ...
    Bandwidth="1/3 octave",FilterOrder=48, ...
    FilterDirection="backward",FrequencyRange=[100,16000])

Check Reliability of RT-60 Values for Forward Filtering

Compute the RT-60 values of a room impulse response using the default forward filtering.

[rir,fs] = audioread("ChurchImpulseResponse-16-44p1-mono-5secs.wav");
bw = "1/3 octave";
frequencyRange = [50,500];
metrics = rt60(rir,fs,Bandwidth=bw,FrequencyRange=frequencyRange);

According to ISO 3382-1:2009, section 7.3, for results to be reliable, $$BT$$ must be greater than 16, where $$B$$ is the filter bandwidth in Hz, and $$T$$ is the measured reverberation time in seconds.

To check if the results are reliable, first get the filter bandwidths. Use the octaveFilterBank getBandedgeFrequencies for a precise calculation.

bwhz = diff(getBandedgeFrequencies(octaveFilterBank(bw,FrequencyRange=frequencyRange)));

Check if the measurements are reliable.

T20_isreliable = metrics.T20 .* bwhz(:) > 16;
T30_isreliable = metrics.T30 .* bwhz(:) > 16;
Topt_isreliable = metrics.Topt .* bwhz(:) > 16;

Create a table indicating which results are reliable. See how the results are not reliable for the lower bandwidths for this room impulse response.

reliableResults = table(T20_isreliable,T30_isreliable,Topt_isreliable, ...
    VariableNames=["T20","T30","Topt"], ...
    RowNames=arrayfun(@(x)num2str(x),metrics.CenterFrequencies,UniformOutput=false))

reliableResults=10×3 table
            T20      T30     Topt 
           _____    _____    _____

    50     false    false    false
    63     false    false    false
    80     false    false    true 
    100    false    false    true 
    125    false    false    true 
    160    true     false    true 
    200    true     false    true 
    250    true     true     true 
    315    true     true     true 
    400    true     true     true 

Analytical and Physical Modeling of RT-60

In this example, you calculate reverberation time (RT-60) of a room using both the Sabine equation and a more physically accurate room modeling approach. In the first section, you use the Sabine equation to estimate RT-60 from room dimensions and absorption coefficients. Then you simulate the room impulse response (RIR) using physical modeling to obtain a more realistic RT-60. In the second section, you invert the Sabine equation to find absorption coefficients required to achieve a desired RT-60 in each frequency band. Finally, you compare the Sabine-predicted RT-60 values to physical modeling.

Sabine Equation

The Sabine equation is a classic analytical formula to estimate the reverberation time of an enclosed space based on its geometry and average surface absorption. The Sabine equation was published by Wallace Clement Sabine in 1900 after extensive experiments in the Fogg Lecture Hall at Harvard University. Sabine discovered that the reverberation time is proportional to the room's volume and inversely proportional to the total effective absorption:

$${RT}_{60}=\frac{\left(55.25/c\right)V}{A}$$

where:

• ${RT}_{60}$ is the reverberation time in seconds.

• $c$ is the speed of sound in meters per second.

• $V$ is the room volume in cubic meters.

• $A$ is the total absorption of room in metric sabins.

The Sabine equation is still in use today and is a valuable tool for its simplicity. However, it has several limitations, such as:

• Assumes the sound field is perfectly diffuse--meaning a uniform distribution of sound energy throughout the room

• Assumes absorption is evenly distributed on all surfaces

• Is less accurate in small and irregular rooms

• Does not account for air absorption, scattering, or diffraction effects

• Can be inaccurate when the mean absorption coefficients are greater than 0.25

For the purposes of this example, the physical room and simulation are contrived to simulate many of the suppositions of the Sabine equation.

Room Modeling for RT-60 Estimation

Define room dimensions in meters. Choose dimensions that are not integer multiples of each other to avoid modal degeneracies in the physical simulation.

roomdims = [9.5,7.1,3.1]; % length, width, height

Specify absorption coefficients for six frequency bands: 125, 250, 500, 1000, 2000, and 4000 Hz. Use the same absorption coefficients on each surface of the room to follow Sabine's assumption of a diffuse acoustic field. The values are moderate and representative of mostly absorptive surfaces that one could find in a theater.

alpha = [0.2045,0.2100,0.1900,0.2097,0.2495,0.2835];

Use the Sabine equation to estimate the RT-60 for each frequency band. The supporting function iRT60model does this.

t60_sabine = iRT60model(roomdims,alpha)

t60_sabine = 1×6

    0.6925    0.6744    0.7454    0.6754    0.5676    0.4996

Next, use a physical modeling approach to simulate the room's impulse response and estimate the RT-60 from the simulated data.

Define a sampling frequency for the simulation.

fs = 24e3;

Define source and receiver positions for the simulation. The height, 1.55 m, is roughly ear level and a common choice. To reduce the dominance of early reflections and to avoid strong room modes, the positions are away from the walls and corners.

sourceloc = [2.1,2.3,1.55];
receiverloc = [7.0,4.9,1.55];

Use high scattering coefficients to simulate a diffuse acoustic field. Use the same six frequency-dependent scattering coefficients for all surfaces.

scatteringCoeffs = [0.8,0.5,0.5,0.5,0.5,0.5];

To simulate the impulse response using a hybrid image-source and ray-tracing algorithm, use acousticRoomResponse. Use the default parameters of the stochastic ray tracing algorithm (MaxNumRayReflections=10 and NumStochasticRays=1280) to model the reverberant tail. Increase the image-source order for a more accurate model of the sound reflections.

ir = acousticRoomResponse(roomdims,sourceloc,receiverloc, ...
    BandCenterFrequencies=[125,250,500,1000,2000,4000], ...
    AirAbsorption=0, ...
    MaterialAbsorption=alpha, ...
    MaterialScattering=scatteringCoeffs, ...
    Algorithm="hybrid", ...
    MaxNumRayReflections=10, ...
    NumStochasticRays=1280, ...
    ImageSourceOrder=45);

To simulate a realistic measurement, add a noise floor to the impulse response. The supporting function iAddNoiseFloor does this.

noiseFloor = -80; % dB
rir = iAddNoiseFloor(ir,fs,noiseFloor);

Estimate the RT-60 from the simulated room impulse response using T30. RT-60 is the time it takes to the energy of an impulse response to decrease 60 dB from its peak. However, actually achieving a 60 dB dynamic range in a measurement is unlikely due to noise floors in most real measurements. Generally, a 20 dB range or 30 dB range (T20 and T30, respectively) is measured and the slope is extrapolated to 60 dB.

t60 = rt60(rir,fs,FilterOrder=64);
t60_physicalmodel = t60.T30'

t60_physicalmodel = 1×6

    0.7015    0.6875    0.7357    0.6725    0.5711    0.5031

Compare the analytical (Sabine) and simulated (physical model) RT-60s by calculating their root mean square (RMS) difference. Display the RT-60 values and RMS difference.

rms_error = rms(t60_physicalmodel - t60_sabine);

figure
bar(["125","250","500","1000","2000","4000"],[t60_sabine(:),t60_physicalmodel(:)],"grouped");
xlabel("Octave Band Center Frequency (Hz)");
ylabel("RT-60 (s)");
legend("Sabine","Physical Model",Location="best");
title("RT-60 Comparison","Total RMS Difference = " + round(rms_error,4) + " s");
grid on

Invert the Sabine Equation: Find Absorption Coefficients from Target RT-60

Assume you want to design the room to achieve specific target RT-60 values in each frequency band.

t60_target = [1,0.8,0.85,0.8,0.7,0.7];

Invert the Sabine equation to solve for the average absorption coefficient needed for each band. The supporting function iInverseRT60model does this.

alpha_target = iInverseRT60model(roomdims,t60_target)

alpha_target = 1×6

    0.1416    0.1770    0.1666    0.1770    0.2023    0.2023

Verify the resulting RT-60 values match the target values using the Sabine equation with the newly computed absorption coefficients.

t60_target_sabine = iRT60model(roomdims,alpha_target)

t60_target_sabine = 1×6

    1.0000    0.8000    0.8500    0.8000    0.7000    0.7000

Simulate the room again with these new absorption coefficients. Again, add a noise floor to simulate a realistic measurement.

ir = acousticRoomResponse(roomdims,sourceloc,receiverloc, ...
    BandCenterFrequencies=[125,250,500,1000,2000,4000], ...
    AirAbsorption=0, ...
    MaterialAbsorption=alpha_target, ...
    MaterialScattering=[0.8,0.5,0.5,0.5,0.5,0.5], ...
    Algorithm="hybrid", ...
    MaxNumRayReflections=10, ...
    NumStochasticRays=100, ...
    ImageSourceOrder=45);

noiseFloor = -80;
rir = iAddNoiseFloor(ir,fs,noiseFloor);

Estimate RT-60 from the simulated response.

t60 = rt60(rir,fs,FilterOrder=64);
t60_target_physicalmodel = t60.T30'

t60_target_physicalmodel = 1×6

    0.9326    0.9112    0.8589    0.7824    0.6930    0.6406

Compute the RMS error between the Sabine-predicted and simulated RT-60s and display the results.

rms_error = rms(t60_target_physicalmodel - t60_target_sabine);

figure
bar(["125","250","500","1000","2000","4000"],[t60_target_sabine(:),t60_target_physicalmodel(:)],"grouped");
xlabel("Octave Band Center Frequency (Hz)");
ylabel("RT-60 (s)");
legend("Sabine","Physical Model",Location="best");
title("RT-60 Comparison","Total RMS Difference = " + round(rms_error,4) + " s");
grid on

Supporting Functions

Add Noise Floor

function ir = iAddNoiseFloor(ir,fs,noiseFloor)
% Adds a pink noise floor to the impulse response to simulate measurement noise.

% Extend IR for ISO 3382 compliance (minimum 1.6 s duration)
flength = ceil(2*fs);

% Zero-pad IR to minimum length
if numel(ir) < flength
    ir = resize(ir(:),flength);
else
    ir = ir(:);
end

% Generate pink noise and normalize rms.
rng default
noise = pinknoise(numel(ir),1);
noise = noise./rms(noise);

% Scale pink noise to desired rms.
noise_rms = max(abs(ir(:))) * 10^(noiseFloor/20);
noise = noise*noise_rms;

ir = ir + noise;
end

RT-60 Model (Sabine)

function t60 = iRT60model(roomdims,alpha,options)
% Calculates RT-60 using the Sabine equation for a shoe-box room.
arguments
    roomdims
    alpha
    options.SoundSpeed (1,1) = 343 % (m/s)
end

Lx = roomdims(1); % Length
Ly = roomdims(2); % Width
Lz = roomdims(3); % Height

% Surface areas
Sxy = Lx*Ly;
Syx = Sxy;
Sxz = Lx*Lz;
Szx = Sxz;
Syz = Ly*Lz;
Szy = Syz;

S = [Sxy;Syx;Sxz;Szx;Syz;Szy];

% Total Volume
V = prod(roomdims);

% Compute the frequency-dependent effective absorbing area of the room surfaces.
A = sum(S.*alpha,1);

% Apply Sabine formula
t60 = (55.25/options.SoundSpeed)*V./A;
end

Inverse RT-60 Model (Sabine)

function alpha = iInverseRT60model(roomdims,t60,options)
% Inverts the Sabine equation to estimate the average absorption coefficient
% required to achieve a target RT-60 in each frequency band.
arguments
    roomdims
    t60
    options.SoundSpeed (1,1) = 343 % (m/s)
end

Lx = roomdims(1); % Length
Ly = roomdims(2); % Width
Lz = roomdims(3); % Height

% Surface areas
Sxy = Lx*Ly;
Syx = Sxy;
Sxz = Lx*Lz;
Szx = Sxz;
Syz = Ly*Lz;
Szy = Syz;

S = [Sxy;Syx;Sxz;Szx;Syz;Szy];

% Total Volume
V = prod(roomdims);

% Total surface area
Stot = sum(S);

% Invert Sabine formula:
alpha = (55.25/options.SoundSpeed).*V./ (Stot.*t60);
end

Input Arguments

rir — Room impulse response

vector

Room impulse response used to estimate the reverberation decay time of the acoustic space, specified as a vector of real numbers. Impulse responses are generally in the two to five second range and should include the noise floor in the end by continuing to capture after the reverberation tail has diminished beneath the noise floor. Measurements for larger spaces with lower noise floors, such as churches, require longer time captures.

Data Types: single | double

fs — Sample rate (Hz)

positive scalar

Sample rate of room impulse response, in Hz. The sample rate must be greater than or equal to 2000 Hz.

Data Types: single | double

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Example: rt60(ir,fs,Bandwidth="1/3 octave")

Bandwidth — Bandwidth of filters (octaves)

"1 octave" (default) | "1/3 octave"

Bandwidth of filters used to measure reverberation time, in octaves, specified as "1 octave" or "1/3 octave".

Data Types: char | string

FilterOrder — Order of octave filters

10 (default) | positive even integer

Order of the octave filters, specified as a positive even integer.

Data Types: single | double

FrequencyRange — Frequency range of filter bank (Hz)

two-element row vector

Frequency range of the octave filter bank in Hz, specified as a two-element row vector of increasing values in the range [20,fs/2]. The default frequency range depends on the Bandwidth. If the Bandwidth is "1 octave", the default frequency range is [125,4000]. If the Bandwidth is "1/3 octave", the default frequency range is [100,5012].

Data Types: single | double

FilterDirection — Filter direction

"forward" (default) | "zero-phase" | "backward"

Filter direction used by the octave filter bank, specified as "forward", "zero-phase", or "backward".

• "forward" — The function applies minimum-phase filters in the standard forward direction.

• "zero-phase" — The function applies zero-phase filters in the frequency domain.

• "backward" — The function applies minimum-phase filters on the time-reversed signal.

Data Types: char | string

NoiseCompensation — Noise compensation method

"truncate-correct" (default) | "subtract-truncate-correct"

Noise compensation method, specified as "truncate-correct" or "subtract-truncate-correct".

• "truncate-correct" — The function uses the Lundby algorithm to estimate the noise floor and the intersection between the reverberation tail and the noise floor. The function truncates and compensates the squared impulse response according to ISO 3382-1 [2].

• "subtract-truncate-correct" — In addition to the truncation and compensation, the function subtracts the noise floor from the squared impulse response. This method does not follow ISO 3382-1, but [4] demonstrated that this method introduced the least systematic noise of several evaluated methods.

In both methods, the function assumes the room impulse response has at least the last 10% of its duration correspond to noise.

Data Types: char | string

Output Arguments

rtsummary — Summary of RT-60 metrics

table

Summary of the RT-60 metrics specified by ISO 3382-1 and ISO 3382-2, returned as a table with the following variables:

• CenterFrequencies — Nominal center frequencies of the octave bands in Hz.

• EDT — Early decay time (EDT), in seconds. The rt60 function computes the EDT by making a linear fit of the decay curve from 0 to -10 dB and extrapolating to the 60 dB decay time.

• T20 — Reverberation time for 20 dB decay (T20), in seconds. The rt60 function computes the T20 by making a linear fit of the decay curve from -5 to -25 dB and extrapolating to the 60 dB decay time.

• T30 — Reverberation time for 30 dB decay (T30), in seconds. The rt60 function computes the T20 by making a linear fit of the decay curve from -5 to -35 dB and extrapolating to the 60 dB decay time.

• Topt — Optimal decay time, in seconds, based on the range of the decay curve that yields the best linear fit. The measurement begins at -5 dB and chooses the range which yields the best linear approximation while remaining 10 dB above the detected noise floor.

• DegreeOfNonLinearity — Degree of nonlinearity, as described in ISO 3382-2, returned as an N-by-4 matrix, where N is the number of octave bands. This value describes the correlation of the linear fit to the energy decay curve. The columns of the matrix correspond to the four decay ranges analyzed: EDT, T20, T30, and Topt. According to ISO 3382-2:2008 B.2, typical values are between 0 and 5.

• Curvature — Level of curvature, defined as 100*|T30/T20 - 1|, as defined by ISO 3382-2. According to ISO 3382-2:2008 B.3, typical values are between 0 and 5.

• C50 — Clarity metric in dB using a 50 ms cutoff, as described in ISO 3382-1. Clarity is a measure of early to late arriving sound energy. C50 is typically used to describe speech clarity.

• C80 — Clarity metric in dB using an 80 ms cutoff, as described in ISO 3382-1. C80 is typically used to describe music clarity.

• D50 — Definition metric using a 50 ms cutoff, as described in ISO 3382-1. Definition is a measure of the early to total sound energy.

• TS — Center time in seconds, as described in ISO 3382-1. Center time is the center of gravity of the squared impulse response. Center time avoids discrete division of the impulse response into early and late periods.

For the EDT, T20, T30, and Topt metrics, if the detected noise floor does not allow 10 dB of headroom over the range, the measurement is invalid and returned as NaN.

References

[1] IEC 61260-1:2014. "Electroacoustics — Octave-band and fractional-octave-band filters — Part 1: Specifications." International Electrotechnical Commission.

[2] ISO 3382-1:2009. "Acoustics — Measurement of room acoustic parameters — Part 1: Performance spaces." International Organization for Standardization.

[3] ISO 3382-2:2008. "Acoustics — Measurement of room acoustic parameters — Part 2: Reverberation time in ordinary rooms." International Organization for Standardization.

[4] Guski, Martin. Influences of External Error Sources on Measurements of Room Acoustic Parameters. Aachener Beiträge Zur Technischen Akustik, Band 24. Berlin: Logos Verlag Berlin GmbH, 2015.

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

Version History

Introduced in R2025a