Revision 8825a502-e6e9-4e79-a11e-d145c3b2d4dd - Signal Processing Stack Exchange

Encouraged by Hilmar, I've decided to update the answer with all the steps necessary to calculate the Reverberation Time from a scratch. Presumably it will be useful for others interested in this area. Obviously it is a simplest approach because more advanced are definitely beyond a scope.

In the beginning you must obtain the impulse response of a room. It can be done on various ways:

 - **Firing a starter gun**, **popping a balloon**, etc. - basically recording any impulsive-like signal with broad frequency content and omni=directional characteristic. This is the simplest method of obtaining impulse response.
 - **Sweep Sine** measurement with usage of either linear or exponential sweep. This method is my favourite as it allows you to extract many different parameters of your system at the same time.
 - **MLS** (Maximum Length Sequence) measurement with this kind of noise and using Hadamard transform to obtain the impulse response. Some comparison with sweep sine technique can be found [**here**](http://ev.fe.uni-lj.si/3-2011/Policardi.pdf).

After getting the impulse response $h(t)$ of your acoustic space you are able to calculate the reverberation time and other speech related parameters such as $C_{80}$, $STI$, $D_{50}$ (Clarity, Speech Transmission Index, Definition), etc.

Very first step is to filter your impulse response in appropriate frequency band. This is due to fact that reverberation time is a function of position in the room and frequency. Obviously we assume that sound field is diffused and at each position in a room decay rate is the same. Thus we must still consider dependency on frequency band. We are using 1/1 octave or 1/3 octave filters. Its desired parameters are defined by [***EN 61260***](https://law.resource.org/pub/in/bis/S04/is.6964.2001.pdf) standard, preferably you always wish to use the filter of class 0. For most of the time Butterworth 3rd order filter is enough (although for frequencies below 125 Hz you might experience some problems with characteristic not meeting specifications). I did whole implementation on my own with few tweaks, but for usual applications widely used [**MATLAB implementation**](http://www.mathworks.com/matlabcentral/fileexchange/22512-nth-oct-hand-arm---ac-filter-tool-box) is good enough.

Next step after filtering the $h(t)$ and obtaining $h_f(t)$, is to make the response as smooth as possible before conversion to logarithmic scale. For that purpose Hilbert Transform is a widely used tool. The goal is to create the analytic signal:

$$h_A(t)=h(t)+j\tilde{h}(t) $$ 

where $\tilde{h}(t) $ is the HT of your filtered impulse response and $j$ is a complex unit. Now, knowing that your analytic signal is complex, two things are coming from its representation:

$$h_A(t)=A(t)e^{j \psi (t)} $$

where $A(t)$ is the envelope of your signal and $\psi (t)$ is a instantaneous phase. Obviously we are interested mainly in envelope (magnitude of analytic signal). Below you can see overlayed filtered response and its envelope:

![enter image description here][1]

One clue for people speaking MATLAB'ish. The `hilbert` function returns the analytic signal, so no need of multiplying by $j$ and adding original $h(t)$. Just use `abs(hilbert(h))` to get the envelope.

Now we have the smoothed version of our impulse response, but it must be smoothed further. It is optional and one might consider it unnecessary. It is based on Moving Average filter of appropriate length $M$. Everyone knows the moving average filter, it is basically a low-pass FIR filter with frequency response given by:

$$H(f)=\frac{\sin (\pi f M)}{M\sin (\pi f)} $$

Plots for different lengths $M$ are shown below:

![enter image description here][2]

And on following figure you can observe the effect of smoothing of the previously calculated energy curve ($A(t)$ in logarithmic scale becomes $E(t)$) for two different values of $M$. For sampling frequency of $48\mathtt{kHz}$ I've used $M=5001$.

![enter image description here][3]

So now we have the Energy Curve (not a Power) which is simply the smoothed Hilbert envelope in logarithmic scale. Conversion to decibel scale is done by equation:

$$E(t) = 20\log_{10} \frac{A(t)}{\max A(t)}  $$

Curve $A(t)$ can be further smoothed for calculations by using the [**Schroeder Integration method**](http://www.ee.columbia.edu/~dpwe/papers/Schro65-reverb.pdf) of your envelope (also known as inversed time integration). This method gives you maximally flat decay curve and is very easy to implement.


$$L(t)=10 \log_{10} \left[ \dfrac{\int_t^{\infty}h^2(\tau )d\tau}{\int_0^{\infty}h^2(\tau )d\tau} \right] $$

Again for MATLAB people:

    L(td:-1:1)=(cumsum(hA(td:-1:1))/sum(hA(1:td)));

Mind that limit of integration (`td`) is equal to $\infty$. It is true when we do not have any environmental noise. Otherwise you will see that decaying sound 'dives into noise level'. It is depicted below. Red curve is based on the correct limit, where decay crosses noise floor. If you choose your limit of integration to be too short, then as in green curve, the estimate is not long enough, so you will get under-estimation while doing linear interpolation. On the contrary, orange curve has too large integration limit and you will get over-estimated. You always want to find the correct limit. If you didn't obtained curve in such way, then please do so. 

![enter image description here][4]

There are many different methods for estimation of integration limit in Schroeder integral, I think that easiest to implement is the one proposed by [**Lundeby et al.**](http://www.ingentaconnect.com/content/dav/aaua/1995/00000081/00000004/art00009)


Regarding calculation of RT itself. It must be done by performing linear interpolation of your decay curve (or should I say, Schroeder curve) with linear function: $L=A\cdot t+B$ on the correct range (which is described below). When it's done, you calculate the reverberation time from equation. You can see that point of intercept doesn't really matter (per se), you care about gradient of your line.

$$RT = \dfrac{-60}{A} $$

where $A$ is a slope coefficient from interpolated line (in dB/s). While you have that, you can also calculate the correlation coefficient of your linear fit (or r-value), that is telling you how good is the fit.

Finally - limits. You cannot really calculate $T_{60}$ having $60 \mathtt{dB}$ of dynamic range as you are trying to do - you always need more (unless you really know what you are doing). So you are doing linear fit on decay curve at following dynamic ranges, and result is being extrapolated, according to *ISO 3382-2*:

 - $EDT$ (Early Decay Time): upper limit is $0 \mathtt{dB}$ and lower is $-10 \mathtt{dB}$. This parameter correlates well with perceived reverberation time. In practice though beginning for the sake of algorithms, people are using interval of $-1 \mathtt{dB}$ and $-10 \mathtt{dB}$ (i.e. in Norsonic analysers).

 - $T_{10}$: upper limit must start at $-5 \mathtt{dB}$ to remove any fluctuations and then lower limit is taken to be $-15 \mathtt{dB}$, but it always must be at least $10 \mathtt{dB}$ above the noise floor. So in fact you need at least $25 \mathtt{dB}$ of dynamic range (or INR) to be able to calculate $T_{10}$ ($5+10+10$).

 - $T_{20}$: upper limit at $-5 \mathtt{dB}$, lower at $-25 \mathtt{dB}$. Minimum dynamic range needed is $35\mathtt{dB}$

 - $T_{30}$: upper limit of $-5 \mathtt{dB}$, lower at $-35 \mathtt{dB}$, with minimum $45 \mathtt{dB}$ of dynamic range.

Why people are using these values? Well, because you very rarely have $75 \mathtt{dB}$ of dynamic range to be able to estimate $T_{60}$. For very smooth decay these values should be equal. Some more detailed analysis of effects can be found in this publication: 

> [**Measuring Room Impulse Responses: Impact of the Decay Range on
> Derived Room Acoustic
> Parameters**](http://purl.tue.nl/352481346918469.pdf)

If you will not perform Schroeder integration, then most probably they will diverge a lot. For very nasty room you can get something like:

![enter image description here][5]

In theory: $EDT$  = $T_{10}$ = $T_{20}$ = $T_{30}$. Although they are not necessarily going to be. For example if your small room is coupled with large one (let's think of small chapel coupled with cathedral by a small doorway), then you will get very sharp decay for smaller one at the beginning of curve, and very long tail in the end:

![enter image description here][6]

To wrap up. According to standards, you should always estimate your reverberation time $T_{N}$ starting from $-5 \mathtt{dB}$ and down to $(-5-N) \ \mathtt{dB}$. Also ensure that you are at least $10 \mathtt{dB}$ above the noise floor. If not, then you must use other estimate. In your case I suggest to use the $T_{30}$.

Below you have example of $T_{20}$ estimation for a troublesome impulse response. $E(t)$ is the energy decay curve obtained by taking the logarithm of impulse response Hilbert envelope. $L(t)$ is the Schroeder integral. Correlation coefficient between linear fit and $E(t)$ on a given range ($-5 \mathtt{dB}$ - $-25 \mathtt{dB}$) for this case is: $cc=-0.9987$. Calculated RT is $T_{20}=1.15 \mathtt{s}$

![enter image description here][7]

I think that should do. If anyone will follow these steps, then he should get accuracy of $0.1 \mathtt{s} $ against widely used [**Dirac**](http://www.acoustics-engineering.com/html/dirac.html). Anyway it must be remembered that even different softwares tend to give different results (i.e. EASERA and Dirac), so you should be totally fine. For more complicated impulse responses with artifacts at the end of impulse response one might get worse performance, but anyway such measurements are not reliable and reflect wrongly conducted experiment.

  [1]: https://i.sstatic.net/IZo1D.png
  [2]: https://i.sstatic.net/vBszu.png
  [3]: https://i.sstatic.net/AZp7m.png
  [4]: https://i.sstatic.net/15esi.jpg
  [5]: https://i.sstatic.net/1uRiX.jpg
  [6]: https://i.sstatic.net/q3hkJ.png
  [7]: https://i.sstatic.net/474by.png