Page  463 ï~~Simulation of Room Acoustics with a 3-D Finite Difference Mesh Lauri Savioja las@cs.hut. fi Timo J. Rinne tri@cs.hut.fi Tapio Takala tta~cs.hut. fi Helsinki University of Technology, Department of Computer Science Otakaari 1, FIN-02150Espoo, Finland Abstract Room acoustics has been digitally modeled with impulse responses calculated by my-tracing or image-source methods. Both of them are based on geometrical acoustics, which is well applicable only to high frequencies, where diffraction can be neglected. Based on previous research on waveguide meshes, we have implemented a finite difference scheme of arbitrary dimension and used it to model sound propagation in three-dimensional room spaces. The method has shown promising results, particularly clearly demonstrating the diffraction effect. Future research includes combining the finite difference scheme to the geometrical acoustics as a low frequency counterpart. 1 Introduction Different methods for simulation of room acoustics have been developed for a long time. In current models there are still some deficiencies, like neglection of diffraction. In this paper we show that 3-D finite difference mesh is good for modeling low frequency behaviour of a room and with suitable computing time and memory requirements. This method is based on digital waveguides used previously to model musical instruments [Smith, 1992]. 1.1 Simulation of room acoustics There are two main approaches to model room acoustics with computers. One is based on the wave equation, the other on geometrical acoustics [Vian, 1986]. Solutions based on the wave equation Finite element method (FEM) and integral equation technique are the main analytical methods. The greatest problems with these are the enormous memory and computing time requirements. The 3-D Finite Difference Mesh presented here is a discrete approximation of the wave equation and is somewhat near to the FEM. Geometrical Acoustics Geometrical acoustics is based on the theory of sound rays. A ray propagates homogeneously in the space and is reflected from a surface with the angle of reflection equal to the angle of incidence. Geometrical acoustics is similar to geometrical optics, which is widely used technique in modeling the behaviour of light. Most of the computer gener ated illumination models are based on ray optics [Glassner, 1990]. There are two main methods in geometrical acoustics. Virtual image source method considers each reflection as a new imaginary source [Kirsztenstein, 1984]. This method gives good results at high frequencies if walls are specularly reflective. It is computationally heavy and for that reason applicable to first few reflections only. The ray tracing method does not give as accurate results as the virtual image source method [Kulowski, 1985]. However, it is possible to calculate responses for longer duration. Another advantage is that with ray tracing also diffuse reflections can be modeled [Naylor, 1993]. There are some weaknesses with geometrical acoustics, generated by the wave nature of sound. For example diffraction and interference cause problems at low frequencies. At high frequencies geometrical acoustics is a valid approach. Hybrid Methods Hybrid methods combine different methods into one model. Most hybrids are based on geometrical acoustics [Vorlander, 1989], [Naylor, 1993]. Due to its limitations low frequencies are poorly modeled. Analytic methods should be used to overcome this problem. When limited to low frequencies the memory requirements of analytic methods are quite feasible. 1.2 Digital waveguide mesh One-dimensional digital waveguides are a discrete numerical method widely used to model musical ICMC Proceedings 1994 463 Acoustics

Page  464 ï~~instruments. Van Duyne and Smith [1993] have presented a two-dimensional extension of digital waveguides. The idea of 2-D waveguide is perpendicular interconnection of 1-D waveguides. The 2-D waveguide mesh is constructed of lossless junctions of equal impedances. There are two main constraints which the junction must satisfy. 1. Sum of the inputs to the junction equals the sum of the outputs, i.e. the flows must sum to zero. 2. Pressures in each waveguide attached to the junction are equal, because of equal impedances. 1-D waveguides have been developed for onedimensional vibrating systems such as strings or pipes. 2-D waveguides give good results for plates, drums, etc. For modeling room acoustics with waveguides, we need a 3-dimensional mesh. 2 3-D Finite Difference Mesh Smith and Van Duyne [1993] derived the lossless scattering equations for the interconnection of several 1-D waveguides: The previous equation is based on the waveguide model, but it is equivalent to a first-order difference equation. Our model is implemented based on that equation. Thus we use the term finite difference mesh, instead of 3-D waveguide mesh. 2.2 Termination at the boundaries In our implementation the termination is carried out by adding one node to the end of each waveguide. Nodes at the boundaries have only one neighbor and they behave exactly as in 1-D waveguides. All other nodes have exactly 2N neighbors Fixed end node At a fixed end (phase reversing reflection) terminat ing node is bound to zero: 1) =0. + D = 1)-i (1) (2) = (n)= " (n-l) (3), where u. represents the signal (pressure or velocity) at a junction at time step n, vi are the incoming waves at the junction, and v{ are the outgoing waves. The reason for equation (3) is that the junctions are interconnected by unit delay elements. 2.1 Nodes in N-dimensional mesh In an N-dimensional mesh each junction has 2N neighbors. If the mesh is considered to be homogenous the lines have equal impedances. So equation (1) reduces to the form Free end node At a free end (phase preserving reflection) the node impedance is zero and the node always assumes the value of its neighbor: UJ, k (n) = u1j,1(n) Matched impedance node Unechoic walls are simulated with nodes of matched impedance. This corresponds to a waveguide that continues to infinity. According to (3) and (4), the equation therefore becomes: uJ, k(n) = ufl(n-1) This equation is useful for simulation of unechoic walls, although it generates minor reflections at high frequencies. 3 Implementation Our generic N-dimensional system is implemented with an object-oriented approach and C++ language. The main goal was to make easily modifiable and flexible software. In the future the execution speed will be optimized and the computation parallelized. Modeling objects with different materials in a space is possible. Material parameters determine the absorption of signal and the behaviour of its phase. Also the mesh itself has a parameter for lowpass filtering during propagation. The mesh can be excited with arbitrary signals. Functions used to calculate signal values at junctions are easy to change. For example the functions can. refer to as many previous time steps as required. This makes higher-order approximations of the wave equation possible. Program development and test runs were done in Sun Sparcstation and Silicon Graphics Onyx workstations. N (4) By using equations (2) and (3), equation (4) can be modified to contain no references to u4 or 1): XD1.(n- 1) vJ,k(n) = t N uJk(n-2) where k is the position of the junction to be calculated and 1 represents all the neighbors of k. Acoustics 464 ICMC Proceedings 1994

Page  465 ï~~4 Empirical Results In the following figures there are demonstrations of different situations with a mesh. The test configuration is a 3-dimensional mesh of 80x80x20 nodes. This corresponds to 8x8x2 meters, if one mesh unit is 10 cm. Three rooms were modeled, separated by unechoic walls and two open doorways. The ceiling is parametrized to produce phase preserving reflection. All the other surfaces make phase reversing reflections. Note that the figures are 2-D slices of the 3-D mesh at the height of the excitation signal. Excitation signal has a Gaussian shape located near to a corner of the largest room. diffracted waves Figure 4:Diffraction from the hole Below are figures of lowpass filtered impulse response at point P1 and the corresponding Schroeder integral, which can be used to estimate reverberation times [Schroeder, 1965]. Point P1 Figure J:Room configuration and excitation signal primary wavefront reflection from the ceiling Figure 5:Impulse response at point P1 Figure 2:First reflection from the ceiling e.7 o.a e.e Figure 6:Schroeder integral for estimation of reverberation time In the following examples the test configuration is formed of a corridor with the sound source near one end and a square-shaped pillar at the other end. The long sidewalls are unechoic. phase reversed reflections Figure 3:Phase reversing reflection from the sidewall ICMC Proceedings 1994 465 Acoustics

Page  466 ï~~Figure 7:Primary wavefront and first reflection from the left wall reflected wavefront _/ Figure 8:Reflection from the pillar diffracted primary wavefront Figure 9:Diffraction effect of the pillar 6 Summary We have derived equations for a 3-dimensional waveguide mesh and implemented an N-dimensional system using object-oriented approach. The model has been shown in practice to be correct. It can be used as a counterpart for geometrical acoustics at low frequencies, where dispersion error and computing time requirements of difference scheme are not disturbing. Acknowledgements The authors would like to thank Matti Karjalainen, Antti Jlrvinen, Vesa VAlimlki, Juha Backman and Jyri Huopaniemi at the Acoustics Laboratory at the Helsinki University of Technology for their help in preparing this paper. This work is financially supported by the Academy of Finland. References [Glassner, 1990] Andrew S. Glassner, An Introduction to Ray Tracing, 3. printing, London, Academic Press Limited, 323 p. [Kirszenstein, 1984] J. Kirszenstein, An Image Source Computer Model For Room Acoustics Analysis And Electroacoustic Simulation, Applied Acoustics, Vol. 17, No. 4, Pp. 275-295. [Kulowski, 1985] A. Kulowski, Algorithmic Representation of the Ray Tracing Technique, Applied Acoustics, Vol. 18, No. 6, pp. 449-469. [Naylor, 1993] Naylor G. M., ODEON - Another Hybrid Room Acoustical Model, Applied Acoustics, Vol.38, Special Issue on Computer Modelling and Auralisation of Sound Fields In Rooms, pp. 131-143. [Schroeder, 1965] M. R. Schroeder, New Method of Measuring Reverberation Time, J. Acoust. Soc. Am., Vol37, pp. 4090-412. [Smith, 1992] Julius 0. Smith, Physical modeling using digital waveguides, Computer Music Journal, VoL 16, No. 4, pp. 75-87, Winter 1992. [Van Duyne and Smith, 1993] Scott A. Van Duyne and Julius 0. Smith, Physical Modeling with the 2-D Digital Waveguide Mesh, Proceedings of the 1993 ICMC, Tokyo. [Vian, 1986] J. P. Vian, Different Computer Modelling Methods - Their Merits and Their Applications, 12th LC.A., Toronto, July 1986, E4-10. [Vorlinder, 1989] M. Vorlflnder, Simulation of the Transient and Steady-state Sound Propagation in Rooms Using a New Combined Ray-tracing/ Image-source Algorithm, J. Acoust. Soc. Am., Vol. 86, No. 1, pp. 172-178. 5 Result Analysis In a 3-D mesh, waves propagate one half of the diagonal unit per time sample. Thus the speed is Jfxd C = 2xT ' where d is the distance between two nodes in the mesh and T stands for one time step in the model. The speed of sound in the air is about 343 m/s at room temperature. The upper limit frequency of the mesh is then: 1 _ 2x343 2x= 2xf/3xd For example, if d equals 10 cm, the limit frequency is approximately 2 kHz. About 8000 time steps are needed to calculate the impulse response for 2 seconds with the previous configuration. The main problem of the finite difference mesh is the dispersion of wavefronts [Van Duyne and Smith, 1993]. The largest dispersion is along the coordinate axes at high frequencies, whereas waves propagate diagonally with no dispersion. For this reason valid frequencies are much lower than the limit frequency given by the previous equation. Ways to reduce the dispersion would be taking a high-order approximation of wave equation or using a denser mesh. Acoustics 466 ICMC Proceedings 1994