A Model for Room Acoustics
Copyright 2003-2004 Stephen G. McGovern, All rights reserved.

Preface: This paper uses a simple image method to calculate a room impulse response. A mathematical sequence giving the position of the nth virtual source is given. Another sequence giving a coefficient for the collective effect on amplitude of all reflections is given as well. Neither sequence uses a permutation, and both use the bottom-front-left corner as a point of origin. Two Matlab functions are linked. One calculates the room impulse response, and the other implements it. In addition, a demonstration MP3 is linked at the end of the paper. Comments and review are welcomed.

Introduction

There are various schemes for calculating a room impulse response. The following model uses an image method. This article will be broken into three major parts so as to clearly explain how the model works. First we will illustrate an approximately equivalent physical situation that will make easier to visualize the individual echoes that together produce reverberation. Second we will find a unit impulse response for each echo with the proper time delay as it would be heard at a particular position in the room. Third we will calculate the magnitude of each echo's unit impulse response. What we are more or less going to do is find the time and magnitude of each echo as it is heard from a particular position in the room. Finally all this information will be put together into a one-dimensional function of time. This function will be our room impulse response. In the impulse response function time can be made discrete. Discrete time will allow us to use it as an FIR (finite impulse response) filter for simulating reverb.

Visualizing the Process

The figure below shows a rectangular room. Within it is a sound source, shown as a green circle, and the location at which we are trying to calculate the impulse response, depicted as a black star. Let's say that this black star represents the location of a microphone. The line between the star and the circle is the path taken by the sound wave. This is the direct sound.

Another part of the sound wave gets reflected off a wall and then impinges upon the microphone. This reflected wave will be referred to as an echo. The listener perceives this echo as radiating from a point past the wall from which it was reflected. Let's say we make a mirror image of the room and placed it adjacent to the original as is shown below in figure 2. If we were at the location of the star we could see the sound source as well as its mirror image. This mirror image will be referred to as a virtual source, and is shown below as a black circle. It is also the location from which we will perceive the echo to be radiating.

In the above figure (citing Pedrotti & Pedrotti, page 44) the black line represents the actual path of the sound wave, while the blue line represents the perceived path of the sound wave. This process can then be repeated by making a mirror image of the room's mirror image. Each mirror image of the source represents another virtual source. A diagram of the real sound source along with two virtual sources is shown in the figure below.


This process can be extended into two dimensions. A map of the virtual sources is illustrated in the figure below (see Allen & Berkley).


Again the star is the listening location and the black circles are the virtual sound sources. This two dimensional model can again be extended and put into three dimensions. In our analysis we will treat the virtual sources as an individual sound sources and ignore each virtual sources corresponding echo.


The Unit Impulse Response Functions


Locating the Virtual Sources

Now consider a three dimensional extension of figure 4. Our goal is to find the locations of the closest virtual sources. To begin let's first look at the model in just one dimension. This is depicted below in figure 5.


The red part is the origin. The x-coordinate of the virtual sources can be expressed using the sequence below.

(1)

xs is the x-coordinate of the sound source and xr is the length of the room in the x-dimension. The location of the ith virtual source is determined by plugging in an integer for i. If i is negative then the virtual source is located on the negative x-axis. If i = 0 then the virtual source is actually the real source. We can find the distance between the ith virtual sound source and our microphone by subtracting the microphone's x-coordinate, xm, from xi. This is shown below.

(2)

The relative positions of the virtual sources along the y and z axes can be found in a similar fashion using equations 3 and 4.

(3)

(4)


Finding the Distance to Each Virtual Source

To find the distance to each virtual source we simply plug xi, yj, and zk into the Pythagorean theorem, as shown below.

(5)

This equation will represent the distances in the form of a three-dimensional matrix.


Finding the Unit Impulse Response Function of Each Virtual Source

Lets say we're given the equation

(6)

t is the time, di,j,k is the distance given by equation 5, and c is the speed of sound. In the above equation di,j,k/c is the effective time delay of each echo. We now must define a function, ai,j,k(u), such that ai,j,k(u) is 1 when ui,j,k(u)=0 and 0 when ui,j,k(u)&ne0.

(7)

The above equation is the unit impulse response function we have been looking for.


Magnitudes of Unit Impulse Responses


Right now each unit impulse response function has a magnitude of one when ui,j,k(u)=0. We will now proceed to fix this. There are two things that will affect the magnitude of our echoes. The first is the distance it travels to get from the source to the microphone. This is given by the following equation (The relation symbol means proportional to).

(8)

The other is the number of reflections the sound wave makes while in transit. If all the wall reflection coefficients are the same, we can take the wall reflection coefficient, say rw and raise it to the exponent n where n is n=|i|+|j|+|k| and represents the total number of reflections the sound has made. This is shown in equation 9 for the virtual source with the indices i, j, and k.

(9.0)

If each wall has a different reflection coefficient then the situation becomes a bit more complicated. If rx=0 is the reflection coefficient for the wall perpendicular to the x-axis that is closest to the origin, and rx=xr is the reflection coefficient for the wall opposite that, then the combined reflection coefficient for all the reflections made by the ith virtual source along the x-axis can be found by using the following equation.

(9.1)

Similarly the combined reflection coefficient for the jth and kth virtual source along the y and z-axes can be found using the respective equations.

(9.2)

(9.3)

To find the total reflection coefficient of a virtual source with the indices i, j, and k we can simply multiply rxi, ryj, and rzk together as shown below.

(9.4)

We now find the total magnitude of each echo by multiplying equations 8 and 9 together as shown in equation 10.

(10)


Construction of the Impulse Response


Finally to obtain the impulse response we multiply equations 7 and 10 together and sum over all three indices. This can be thought of as the summation of the all the sounds as they stream from all of the virtual sources. This is shown below in equation 11.

(11)


Implementation of the Room Impulse Response


Model

RIR is a program that calculates our room impulse response h(t). This program differs from the model in two ways. First it uses discrete time instead of continues time. Second, it finds ai,j,k and ei,j,k only when ai,j,k=1. This is done to conserve memory.

Simulate

Ordinarily an FIR filter modifies a sound file via a convolution. Our filter h(t), however, is extremely large. An ordinary convolution is much too slow. To use the filter one must perform a fast convolution. To find out what a convolution is and get a superficial explanation of what a fast convolution is go here. A program called FCONV does a fast convolution and is provided at the end of this file.


Discussion


In figure 6 a circle has been drawn concentric with the microphone and so that its edge touches the nearest wall. When extend into three dimensions figure 6 would form a rectangular prism with a sphere inside it. Only the portion of the impulse response that comes from virtual sources with in this sphere are valid. This comprises the first part of the impulse response. The remaining virtual sources that are outside the sphere have counterparts that have been omitted in the Matlab version of RIR. This means that the latter portion of the impulse response will have significant error if the summation indices in equation 11 are set equal to one another. The stand-alone Windows and Linux versions of RIR do not produce this error.

This paper makes a few assumptions that are known to cause error. It assumes that the reflection coefficients are angle and frequency independent. It also assumes that no change in phase takes place upon reflection. Additionally it assumes that air has no effect on the magnitude of the sound wave. There may be other errors as well as I am by no means an expert in the area of acoustics.


Conclusions


I suspect that the most perfect model would have to use the leapfrog method (A recursive method used in computational physics). The leapfrog method could probably simulate small room acoustics on a desktop with less that 500 MB of memory. I can not, however, comment on the time it would take to execute such an algorithm. For larger, say gymnasium size rooms, one would likely need upwards of 30 GB of memory. A solution in the frequency domain may also prove to be more realistic.

The Matlab function RIR.m is equivalent to this model and uses relatively simple code. The program produces data that is virtually identicala to that produced by Allen and Berkley's1 code. The differences are the output data format, a uniform difference is magnitude, and error in the tail. The uniform difference in the magnitude can be eliminated by multiplying all the output data points by the same coefficient. This, however, has not been checked for the possibility of small rounding errors. Likewise the error in the tail can be corrected by finding a longer impulse response and truncating the trailing end.

The stand-alone Windows and Linux versions of RIR are also equivalent to this model. They are faithful to Allen and Berkley's code in that they do not produce any error in the calculation of the tail. The actual algorithm used in these versions, however, has been optimized for speed and differs markedly from the one outlined in this paper.


[a] This also presumes that the highpass filter is removed from the end of Allen and Berkley's code.



RIR - A command line utility that uses the mirror image method to calculate a room impulse response function. The output can be written to the hard drive in the form of a .WAV file. It is available for Windows, Linux, and Matlab.

Matlab download - For documentation put this in your Matlab work folder and type "help rir" at the Matlab command prompt.
Screenshot (in an audio editor)


FCONV - An open source command line utility that uses an FFT to convolve two .WAV files. The output is then written to disk as a third .WAV file. FCONV is means of implementing FIR filters on audio files. It is available for Linux and Matlab.

Matlab download - For documentation put this in your Matlab work folder and type "help fconv" at the Matlab command prompt.


Bibliography

[1] Jont Allen and David Berkley.
Image Method for Efficiently Simulating Small Room Acoustics.
Journal of the Acoustic Society of America,
p912-915, April 1979.

[2] Coppers, Frey, Kinsler, Sanders.
Fundamentals of Acoustics.
Wiley, 3rd edition, 1982.

[3] Eugene Hecht.
Optics.
Addison Wesley, 4th Edition, 2002.

[4] Scott Lehman.
Reverberation.
http://www.harmony-central.com/Effects/Articles/Reverb
Effects Explained.
Harmony Central 1996.

[5] Sanjit Mitra.
Digital Signal Processing,
A Computer-Based Approach.
McGraw-Hill, 2nd edition, 2001.

[6] Pedrotti and Pedrotti.
Optics and Vision.
Prentice Hall, 1st edition, 1998.

[7] James Stewart.
Calculus: Concepts and Contexts.
Gary W. Ostedt, 1998.

[8] Udo Zölzer.
DAFX, Digital Audio Effects
Wiley, 2002.