Copyright © 20032004 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 n^{th} 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 bottomfrontleft 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. IntroductionThere 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 onedimensional 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
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 xcoordinate of the virtual sources can be expressed using the sequence below.  
(1)  
x_{s} is the xcoordinate of the sound source and x_{r} is the length of the room in the xdimension. The location of the i^{th} virtual source is determined by plugging in an integer for i. If i is negative then the virtual source is located on the negative xaxis. If i = 0 then the virtual source is actually the real source. We can find the distance between the i^{th} virtual sound source and our microphone by subtracting the microphone's xcoordinate, x_{m}, from x_{i}. 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 x_{i}, y_{j}, and z_{k} into the Pythagorean theorem, as shown below.  
(5)  
This equation will represent the distances in the form of a threedimensional matrix.  
Finding the Unit Impulse Response Function of Each Virtual Source  
Lets say we're given the equation  
(6)  
t is the time, d_{i,j,k} is the distance given by equation 5, and c is the speed of sound. In the above equation d_{i,j,k}/c is the effective time delay of each echo. We now must define a function, a_{i,j,k}(u), such that a_{i,j,k}(u) is 1 when u_{i,j,k}(u)=0 and 0 when u_{i,j,k}(u)&ne0.  
(7)  
The above equation is the unit impulse response function we have been looking for.  
 
Right now each unit impulse response function has a magnitude of one when u_{i,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 r_{w} 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 r_{x=0} is the reflection coefficient for the wall perpendicular to the xaxis that is closest to the origin, and r_{x=xr} is the reflection coefficient for the wall opposite that, then the combined reflection coefficient for all the reflections made by the i^{th} virtual source along the xaxis can be found by using the following equation.  
(9.1)  
Similarly the combined reflection coefficient for the j^{th} and k^{th} virtual source along the y and zaxes 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 r_{xi}, r_{yj}, and r_{zk} 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 a_{i,j,k} and e_{i,j,k} only when a_{i,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.  
 
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 standalone 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.  
 
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 identical^{a} to that produced by Allen and Berkley's^{1} 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 standalone 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
