Common methods in acoustics include the finite element method (FEM), the boundary element method (BEM), and finite difference (FD) methods. Some methods seek solutions to the stationary problem in the frequency domain, which is useful for describing modal distributions in low frequencies. However, many interesting and useful properties of sound fields are captured or encoded in the transient, time-domain response. 

In room acoustics and more generally linear acoustics, the computational challenges are not the same as many other disciplines.  The physics of linear wave propagation is relatively straightforward, but the typical size of a problem and necessary boundary conditions present unique challenges.  For instance, the volume of a concert hall can be on the order of 10,000 m³, while the wavelengths of sound in the audible range are on the order of centimeters.  Capturing the wave physics on these scales requires massively parallel computing. Further inhibiting efficient computation, realistic geometries and boundary interactions also require careful design.  Ever-increasing computing power and demand for physics-based simulation make wave-based modeling a continually active area of research.

Finite-Difference Time-Domain Technique

The approach that is currently most promising for wide-band room acoustic simulation is the finite-difference time-domain (FDTD) technique. The first papers to discuss 3D, wave-based simulation in rooms, using FDTD date back to 1994 [1, 2].  Ever since, there has been a division between two separate approaches, and both of them are still valid and in use (see e.g., [3, 4]). The approach presented by Botteldooren [1] is based on the linearized Euler equations and it utilizes two variables to present the sound field, most typically sound pressure and volume velocity. The technique originates from the electro-magnetic literature where Yee presented this technique for simulation of electro-magnetic fields already in 1966 [5]. This approach has come to be known as the standard leapfrog scheme.

The other option is to start from the second-order, scalar wave equation, which requires only one field variable, such as sound pressure [2, 4] or velocity potential. This saves both in memory consumption and computation time, while still providing the same accuracy as the standard leapfrog scheme. The video below is computed with this technique.

A major problem in all the wave-based techniques is their heavy memory consumption at higher frequencies. Finite difference techniques require discretization in both time and space.  In particular, the room under study must be divided into small volumetric elements, called nodes, and the simulation progresses in discrete time steps in which all the nodes are updated according to given equations. The relation of the temporal and spatial discretization step sizes is dictated by the Courant condition [6].

Finite difference solvers can be further divided into explicit and implicit schemes. In the explicit case, the new state of a node can be computed based on its own history and on the past values of the neighboring nodes. This means that it is possible to perform all the updates in parallel thus suggesting it to be a good fit for parallel computation. In the implicit FDTD schemes, the new values of the nodes depend on the new values of their nearest neighbor nodes. Although there are ways to make this efficient, this approach is currently not favored as a general simulation tool. The main advantage of the implicit schemes is that they have less strict conditions on stability and they can be more accurate than explicit schemes, but due to the computational requirements the explicit ones are favored at the moment.

One of the first 3D FDTD implementations for room acoustics

The FDTD simulation for room acoustics has gained a great deal of popularity in recent years, but our earliest papers in the topic date back already to 1994 [2]. That paper was insipred by the 'Digital Waveguide Mesh' approach originally targeting physical modeling of musical instruments [7]. Although the technique is derived from signal processing viewpoint, it is equivalent to an FDTD in which the relation of the temporal sampling to the spatial sampling is set to the Courant limit [6].

Interpolated FDTD schemes

A fundamental problem in the FDTD simulations is dispersion that depends both on the frequency and on the propagation direction. In a rectangular mesh this is seen such that high frequencies get delayed in axial propagation directions in relation to the underlying grid whereas in diagonal directions all the frequencies propagate at the same constant speed. To overcome this limitation of direction-dependent dispersion, we introduced the interpolated digital waveguide mesh structure [8]. The proposed technique results in nearly isotropic propagation. This approach was later developed further by Kowalczyk et al. [4].

Frequency warping with FDTD

While the interpolated mesh improves over the original standard rectilinear scheme, it still has frequency-dependent dispersion error. Frequency warping [9] is a technique that can be used to post-process the obtained results such that this dispersion error is minimal [10, 11]. With this approach the usable frequency band in FDTD simulations can be remarkably increased [12].

Parallelization of a FDTD solver for real-time room acoustic simulation

FDTD is a technique that is easy to parallelize. Introduction of modern GPUs has had a remarkable impact on practical parallel computation in several areas including signal processing and acoustics [13]. Current GPUs are massively parallel processors that have made room acoustic simulation with FDTD a very attractive option. For small rooms it is possible to have even real-time performance for low- and mid-frequency simulation [14].

Frequency band extrapolation

As the computational load in FDTD simulation grows with respect to the fourth power of the frequency, simulation of high frequencies is still a major challenge for current computers. For this reason, we have developed an extrapolation technique that enables full-bandwidth auralization [15]. It is based on analysis of responses such that each reflection is detected from a spatial response, it's direction is estimated, and then extrapolated for full bandwidth. The process included simulation of air absorption that has a big role in sound attenuation at higher frequencies.

Ray Tracing

One of the most widely used modeling technique is ray tracing [1]. A large number of rays are emitted from the sound source. When a ray hits a surface of the room, it is reflected and the ray continues in an other direction. Usually the ray is reflected in the specular reflection direction, which means that the angles between incoming and outgoing rays with respect to the surface normal are equal. Each ray is supposed to carry a small portion of the sound energy emitted by the source. Some of that energy is absorbed by the surface at reflections. The portion of the energy absorbed is determined by material- and frequency-dependent absorption coefficients. Some of the rays eventually hit the receiver, which is modelled as a finite-size volume, such as a sphere. Those rays contrubute to an echogram, which collects the arriving energy in a time-dependent response. The arrival times are calculated by dividing the ray path length from the source to the receiver by the speed of sound.

The ray tracing technique can be extended to model non-specular reflections by generating the outgoing rays at reflections randomly from an angle distribution, instead of using only the specular direction. In this case, the energy is divided between the specular ray and other rays according to a scattering coefficient.

Ray tracing can be seen as a stochastic method, where the space of possible reflection paths is sampled with rays. This also means that some paths might be missed due to low sampling density or their contribution to the echogram can be miscalculated. Increasing the number of rays or splitting them later during the tracing usually improves the quality of the results.

Image Source Method

Another common modeling technique is the image source method [2], [3]. Image sources are created by taking a mirror image of the sound source around the planes at which the surfaces of the room lie. Each mirror image is treated as a source which emits sound like the original source. From the receivers point of view, the sound emitted by these image sources corresponds to the first order reflections from the surfaces which were used in creating those image sources. The ray path can be reconstructed by starting a ray from the receiver point towards an image source, and when the reflecting surface is hit, then continuing the ray towards the original source. Higher order reflections can be modelled by creating image sources of the image sources recursively, and tracing the paths backwards from the receiver in a hierarchy of image sources.

The validity of the reflection paths created as above must be checked before constructing an impulse response out of them. There are two ways a path may by invalid: the ray does not hit the reflecting surface corresponding to an image source or the ray hits another surface before hitting the reflecting surface. It is possible to use ray tracing to test this, and discard false paths.

After the validity of the paths is confirmed, they can be added to the room impulse response. The path length is related to the arrival time as in ray tracing. The amplitude of the peaks corresponding to the reflections is inversely proportional to the path length. Unlike in ray tracing., it is also possible to account for the phase information, since the specular reflection paths are accurately modelled.

One limitation to the image source method is that the number of image sources grows exponentially in relation to the order of reflections. Thus, it is only feasable for low order reflections.

Beam Tracing

The image source method can be accelerated by taking into account the validity of the paths already during the creation of image sources. This can be done with a technique called beam tracing [4]. Instead of starting rays from the sound source, cones or beams are sent forward from it. In an ideal case, the full solid angle around the source is divided between the beams so that none of them overlap or leave gaps between them. A convenient way of doing this is using the source as the apex and the polygons visible to it as bases of the cones. Then the beams are traced further and flipped at the surfaces so that each beam edge is reflected like specular rays in ray tracing. Incidentaly, the apex of the cone is then the image source corresponding to the surface at which the beam is flipped. When a beam hits multiple surfaces, it is split so that the parts of the beam covering different surfaces are put into different beams. Those beams are further flipped and traced. The beams are stored in a hierachical structure called beam tree. At the first level of the beam tree are the initial beams. When the beams are flipped and split, new levels are constructed in the tree under the previous level of beam. This structure remains unchanged for a fixed sound source position, and it can be precalculated and stored.

To find the reflection paths, the first step is to examine in which beams the receiver position lies. When those beams are found, the paths can be generated as in the image source method by tracing backwards from the receiver through the image sources and eventually to the source. The image sources are either explicitly stored in the beam tree or they are quickly calculated from the list of surfaces that are hit during the tracing of a given beam. Although the handling of beams requires extra effort compared to the image source method, the beam tree is much narrower than the hierarchy of image sources of the same maximum reflection order. In addition, if the beam splitting is done accurately, there is no need for validity checks.

The beam tracing technique usually allows higher order of reflections to be modelled with the same computation time when compared to the image source method. The benefit is greatest in densely occluded environments such as office spaces with narrow corridors. In large open spaces, like in the concert halls, large part of the image sources are valid, and thus the performance gain of beam tracing is less significant. On the other hand, beam tracing is efficient in interactive applications, where the sound source is fixed, but the receiver is moving, because the reflection paths can be easily updated by utilizing the precomputed beam tree [5].

Radiance Exchange Methods

While specular reflections are physically correct representations of reflections from rigid flat surfaces with infinite extent, in practice some surfaces produce more diffuse-like reflections. Diffusors in concert halls are an example of such surfaces. They can be modelled either with a modified ray tracing algorithm as mentioned earlier, or with an algorithm specialized in diffuse reflection modeling. They can be used in combination with the other algorithms or as such for room acoustics modeling [6], [7].

The fundamental idea behind this class of modeling algorithms is that the surface of the room is divided into small elements and acoustic radiance exchange between them is modelled. Initially, the sound source radiates onto the surface elements. Then each of these elements reflects the incoming radiance in all directions according to the diffuse reflection model. This radiance arrives at other elements and is reflected further. The radiance exchange between the elements can be explicitly calculated by utilizing their geometric properties [7]. Part of the radiant energy is absorbed at each reflection and the total energy in the system decays as the radiance is iteratively propagated between the elements. When the total energy is below a threshold, the propagation can be terminated. In the end, each element has stored an echogram based on the radiant energy that it has received over time. Then, a final gathering step is performed, where the stored echograms contribute to the echogram at the receiver position according to their geometric relationships.

The results of the radiance exchance methods resembles that of ray tracing, because the final result is an echogram. On the other hand, because the number of elements is fixed in the beginning, these methods do not suffer from the exponential explosion in the number of image sources. Thus, they are well fit for modeling the late reverberation. The main limitation is the reflection model, which is either ideal diffuse or directional diffuse [7]. Early reflections from rigid flat surfaces are better modelled with the other methods.

Diffraction Modeling

As mentioned earlier geometrical room acoustics ignore the wave nature of sound. While at high frequencies the ray approximation is quite good, the low frequencies are most often incorrectly modelled, especially in the presence of diffrating edges. This error can be compensated for by adding a separate diffraction model on top of the geometric model. There are several diffraction models with different levels of accuracy and computational efficiency.

There is an explicit formula for the diffraction from a point source via an infinite edge towards a receiver [8], [9]. This formula can be re-formulated as a line integral which makes it possible to model finite-length edges [10]. Applying this model in practice requires numerical evaluation of the line integral and taking into account special cases [11]. In many cases where interactive performance is required, this is too slow, although the results are accurate.

Another approach is using the so-called Uniform Theory of Diffraction [12] which provides a more efficient way of calculating the diffracted part of the impulse response in the shadow zones.

The Room Acoustic Rendering Equation

We have derived a theory for geometrical room acoustics that covers all the previously mentioned modeling methods as special cases. It relates the outgoing acoustic radiance from any point on the surface of the room to other surface points. This is expressed as a surface integral. The name room acoustic rendering equation comes from its time-independent counterpart in computer graphics. [7]

Acoustic Radiance Transfer Method

We have proposed a novel room acoustics modeling method based on the room acoustic rendering equation. It is a radiance exchange method that uses directionally dependent reflection models. They are represented as bi-directional reflectance distribution functions (BRDFs), which relates the incoming acoustic irradiance to the outgoing acoustic radiance for each incoming and outgoing angle.

The modeling algorithm is called the acoustic radiance transfer (ART) method. The surface of the room is divided into elements and for each of those, the hemisphere above it is further split into directional slots. In the ART method, the acoustic radiance is propagated iteratively from element to element and from directional slot to directional slot, until the propagated energy is neglible. The initialization and final gathering phases are as desribed above in the case of radiance exchange methods in general.

The ART method produces results that are similar to the measured ones when compared with the standard room acoustical parameters. It provides flexibility in modeling different scattering materials through the use of BRDFs. [7]

Accelerated Beam Tracing

When a beam tracer is used as a part of an interactive application, where the reflection paths are calculated for a moving listener, there is coherence between the reflection paths of adjacent listener positions. This coherence can be utilized in the path calculation algorithm.

We have developed a beam tracing algorithm that does not do occlusion tests during the beam tree construction, but postpones the testing as long as possible [5]. It also does not store the clipped polygons, but only references to the original polygons, in the beam tree. This way the beam tree remains simple and it can be fairly deep. The occlusion testing is done during the path construction phase and even then it can be skipped in most cases.

The main idea is to store so-called "fail planes" for each beam. A fail plane is a plane in the hierarchical beam tree against which the path construction has previously failed. This plane is propagated to the lower levels of the beam tree by reflecting it around the surfaces like the beams. Thus, because the path is first tested against the fail plane, the path validation can be quickly terminated already at the lower levels of the beam tree in the case of an invalid paths. Because most paths tend to be invalid, this optimization can tremendeously speed up the calculation. [5]

Another optimization we have proposed for the beam tracing algorithm is the use of a "skip sphere". The beam tree nodes are grouped into "buckets" containing a few nodes. If all of them produce invalid paths, then a skip sphere is calculated around the listener and while the listener stays side that sphere all the nodes in the bucket can be skipped. The radius of the skip sphere is the minimum distance from the listner to the fail planes of the nodes in the bucket. [5]

With these optimizations, the beam tracer algorithm is acapable of updating the early reflection paths at interactive rates for a moving listener in fairly complex models.