t60.m

Returns an estimate of t60, a measure of reverberation time.

Description

Pass a recorded noise burst or impulse, and this script returns an estimate of t60. Pass only the part of the recording after the offset of the noise. It uses Schroeder's integrated impulse decay method, and then performs a cumulative r-squared to find the straight part of the slope. eg, if you recorded a noise segment that is 8192 samples long, syntax would be:

rt = t60(noise(8193:length(noise)),srate);

The script will also plot the decay function if desired.

• MATLAB release: MATLAB 5.3 (R11)

[Olivier Deille]

It seems there are some mistakes in this program, thus giving wrong results.

[Matthias Geier]

Hi Christopher.

Thanks for this really useful function!

I have 2 questions:

1) Can you please provide the original formula for the "r-squared" you used?
It seems to be related to the squared "sample Pearson correlation coefficient" but I think it could be off by a division by "n".
I searched in some statistics books and on the web and the problem is there are many different definitions of "r-squared".

2) By default, your function discards the first 50 ms of the result of the "cumulative r-squared". Is this a value you found by experience? In which situations would you suggest different values?

TIA,
Matthias

[Matthias Geier]

It's me again ...

Finally, I found the problem in your code!

First, I was wrong with my speculation that the r-squared formula would be off by a division by "n".
The formula is OK (it's just the squared version of Pearson's product-moment correlation coefficient).
The following regression is also correct.

However, there is an error when calculating the intersection between the regression line with the -60 dB line.

Your code is:

rt = round(abs(60/slope(1))-yintercept);

This works if yintercept is (very close to) zero (as it is in the nice picture above. But in the general case (and in my case), yintercept is a non-negligible negative value and the resulting reverberation time is wrong.

The correct line would be:

rt = -(yintercept + 60)/slope;

As I said, the regression is correct, however it could be written a little more compactly like this:

xx = xtime(1:stop);
yy = y(1:stop);
slope = (mean(xx.*yy) - mean(xx)*mean(yy))/(mean(xx.^2) - mean(xx)^2);
yintercept = mean(yy) - slope * mean(xx);

That's it about my first point, I'm still searching for an answer regarding my second question.

cheers,
Matthias