Page  00000001 TOWARDS THE SYNTHESIS OF WAVEFRONT EVOLUTION IN 2-D Georg Essl* ABSTRACT The dynamics of drums can be described by wavefronts and their wake on a bundle of rays. Using the ray picture in their construction preserves the structure of the wave propagation, which in general is dispersed by traditional meshed methods. The idea is based on emphasizing local dynamical descriptions and impulse responses of continuous solutions in the plane and translating these properties into a computational model. 1. INTRODUCTION This work is part of a larger program trying to understand if we can find efficient computational structure in the 2 -dimensional wave-equation'. The previous results in this program can be found in [4, 6] and discuss geometric features of wave fronts rendered with the algorithm to be described here. Aspects of numerical behavior as well as the nature and impact of wakes - originally to be presented here - have been moved to an alternative publication [5] due to space restrictions. In [4] I described the basic algorithm and related the results of the algorithm to known features of ray dynamics on a circular domain. In [6] the domain was altered to an ellipse and a stadium (two semicircles connected with straight lines) to give an argument relating the observed wave fronts to the whispering gallery phenomenon. Here I would like to bring the argument on the one hand close to the full acoustics of drums in the plane. On the other hand I want to give a first rendering of this algorithm for sound synthesis. For this, a number of pieces are missing in [4, 6] to arrive at a realistic simulation for sound synthesis. Most notably dissipation and radiation has not been modeled so far. Hence we are adding new structure and information to the algorithm to complete the picture. This is however still not an exact numerical integration of the wave equation as wakes have been omitted. See [5] for details. Digital waveguides [14] have proven to be an extraordinarily successful way of simulating the wave equation in one dimension. This has to do with the particular structure of the solution space of the wave equation, when written as traveling waves (d'Alembert's original solution). This form allows for a computational structure to be utilized that is, in a sense, better than a naive estimation of *The author is currently without affiliation. Please use the listed email for correspondence. SThroughout this paper dimensions refer to spatial dimensions of the domain. The additional temporal dimension is implied. a meshed method to integrate the wave equation. By requiring only constant, rather than order of spatial samples, waveguides are optimal with respect to spatial integration. The second main advantage is that the traveling waves are not, or only mildly disturbed by propagation. This allows for copy-propagation rather than arithmetic propagation. Hence rounding errors are avoided or kept limited to few operations thus giving waveguides highly desirable numerical properties as well. At the same time waveguides are exact with respect to the continuous equations on sampling points. We are still far from having the same properties for the 2-dimensional case. A number of methods have been employed to simulate drums. These are all predominantly meshed methods. These include spring-based methods [2], standard finite differencing methods, and waveguide meshes [16, 8] and hybrids between waveguide and finite differencing methods [1, 11]. All these share that the dynamics is discretized on mesh points and different directions of propagation are coupled at these mesh points. The difference either consists in the motivation or implementation. The mesh itself is responsible for inexactness in the solution known as mesh dispersion and various attempts have been proposed to reduce this effect for various situations [12, 9]. Another approach is using functional transforms[15]. This is really a modal methods, in which the modes-strengths are calculated from the equation once transform functions have been picked. As these transform functions usually constitute a support that ranges over the whole domain of solution, local aspects of the problem become approximated by truncated infinite series of these transform functions. Banded waveguides have also been proposed for 2-dimensional structures such as drums and cymbals [7]. This is a modal method which tries to retain some notion of spatial information. However this information is asymptotic and hence incomplete. Secondly, it is very difficult to recover the meaningful spatial positions of various waveguide bands crossing. Recently a related method, which employs details about the circular symmetry have been proposed [17] as well as hybrids with waveguide meshes[13]. 2. APPROACH The method to be described here does not solve the full problem. It rather is an attempt to understand the structure better. It differs from earlier work in that it does not

Page  00000002 employ a static grid, does not assume specific functional forms of the solution nor does it only consider a limited modal set of propagating paths. However, this choice is not a primary one, but rather a consequence of trying to maintain as much of the structure of the solution in numerical simulation. The key approach is to try to model the geometric and functional content of what is known as the fundamental solution of the 2-dimensional wave equation. The fundamental solution can be thought of as the continuous function-theoretic analogue of the discrete impulse response as familiar in digital signal processing. There is no space here to discuss the theory of fundamental solutions of partial differential equations and we refer the curious reader to [3]. Further literature can be found in [5]. The Diracdelta within this theory is described using the functiontheoretic concept of a distribution. These are explained in details in [3] and other sources (see [5]. However fundamental solutions of the wave equations in one to three dimensions were already known before the theory of distributions was fully developed. For relevant information we refer the reader to [10]. The traveling wave solution, that is so successfully employed in waveguide synthesis, is in fact the fundamental solution of the 1-d wave equation (compare [3, p. 144]). The main difference between the traditional simulations, and also traditional analytic solution methods and work with fundamental solutions, is that it described the behavior of the problem directly by answering how a solution behaves locally both in space and time. Hence if a local disturbance is given, the question how this disturbance will evolve can be answered immediately. In comparison, the traditional solution of the circular membrane is written down in terms of functional expansions using Bessel functions arrived at via the method of separation of variables. However, Bessel functions, as are trigonometric functions in case of the 1-dimensional equation, are distributed qualities and hence it requires additional effort to recover local behavior from these description. Often the correct Fourier-Bessel series is not easily found. With respect to meshed methods, meshes are often chosen independent of any knowledge of the structure of the solution. Rather they are chosen to approximate the form of the (partial) derivative. However, similar partial derivative operators can have different resulting solution spaces, as is already visible in the case of the 1-D versus the 2-D wave equation. 3. STRUCTURE OF THE WAVE EQUATION The wave equation in two dimensions without external forces and dissipation reads: H(~ct - x\)ct y(x, t) = 3 2w r ((ct)2 _ j2 H(+fct - |x|) y(x,t) c (ct)2 2 (2) (3) where 2 is the response to an impulsive displacement and 3 is the response to an impulsive velocity. H(.) is a step function also called the Heaviside distribution. These equations have two parts to them. One is the support of the solution. This part is encoded in the Heaviside function of equations (2-3).The other is the functional shape created by an impulse through propagation inside this support. The first part encodes what we will call the wavefront. Here we will discuss this part of the problem. The second part will be called the wake. A full treatment of wakes is future work. A theoretical discussion of the impact of its omission can be found in [5]. Wave fronts are the point of first arrival of a disturbance in response to an impulse (or a convolution of impulses). The argument of the Heaviside functions in (2) and (3) forms a circle for constant time. The radius of the circle expands linearly with time. Hence the wavefront inhabit a cone with the tip at the point of the impulsive excitation. Any section of this cone will give a wavefront at a certain time, which is a line embedded in the plane. Point of the wavefront propagate along straight lines called rays. Hence one can think of rays as the direction of propagation of waves and wave fronts as the points of arrival on them. The crux of the proposed method is to maintain the support of wavefronts on rays. Hence the position will be accurate, up to numerical accuracy on such rays and there is no numerical dispersion by construction. This method is in addition excitation-point centric, hence solves the problem of banded waveguide-like methods in this respect. For space reasons I cannot present the details of the algorithm here. These details can be found in [4, 6]. Using this method also complicated shapes can be simulated. If the domain is sufficiently simpled the reflection code can be replaced by an iterative function of reflection angle and reflection point. For details see [4]. 3.1. Damping, Radiation and Pickup There are two sources of dissipation modeled. One is geometric. If a wave front travels outward it forms an ever increasing circle. The total amplitude will be diluted in proportion of the spreading of the circumference of this circle. As the increase of the circumference is linear with time, i.e. 2-t the reduction of local amplitude will linear as well. On a closed domain the overall amplitude, however, remains constant and hence no energy is dissipated. The second source of dissipation are friction effects of various sorts. In the current implementation we only em 02 2Y 2 Y t2- c X 2 + X22 9t2 \9x12 9x22} 0 (1) The fundamental solution of the wave equation in the plane without boundaries is [10, 3]:

Page  00000003 (9 cl Figure 1. Examples of wave front shapes on a circular domain. Top: Excitation point at 1/3 of the radius. Bottom: Excitation point 9/10 of the radius. Minimum order of reflections are 1, 7 and 13 (from left to right). ploy a very simple heuristic model to implement some dissipation. This model assumes that the amplitude additionally decreases linearly with distance according to some proportionality factor. By choosing a suitable constant, reasonable first-order damping behavior can be achieved. A basic problem for 2-dimensional simulation, is the problem of the pickup. Where should the virtual ear, accepting a one-dimensional audio-stream, be placed? For the purpose of this presentation, a simple conic radiation method was used. The listener is assumed to be placed directly above the center of the drum at some distance h,. The position x, y on the membrane of a disturbance induces a certain travel-length d(x, y) to the position of the listener according to the basic length of a line on the cone. d(x, y) = x2+y2+ h, (4) The delay to the listener can then be calculated using the speed of sound in air and can appropriately be implemented using a delay line, which is loaded at the correct distance according to equation (4). 4. RESULTS 4.1. Geometric Results A detailed discussion of geometric features of this method for the case of the circle can be found in [4]. A few examples of shapes of wave fronts on a circular drum are depicted in Figure 1. One can see that excitations close to the wall yield a pattern of dense wave fronts creeping along the wall. For this reason these types of wave fronts have been associated with whispering gallery effects [4, 6]. Here we will illustrate examples of alternative domain shapes as seen in Figure 2. All simulations were performed with 500 rays. This includes both generated sound files and graphical renderings. 4.2. Sound Synthesis Results In order to get a qualitative comparison of the proposed method with reality, simple experiments were performed Figure 2. Examples of first order (top) and third order (bottom) reflection wave front shapes on (left to right) rectangle, ellipse and stadium domains. Source fih:fo f2:f 0 f3f 0f4h:fo Simulation 2.7857 3.6429 4.6429 6.5000 Recording 2.1333 3.7333 4.8667 6.2000 Theory 2.294 3.598 4.903 6.208 Table 1. Comparison of modes between simulation, recording and theory for a strike at the center of the drum. with a real drum. All the recordings were performed on a Pearl tom drum with a Pearl ProTone drum head. The bottom drum head was removed to eliminate beating due to mistunings and to reduce the impact of the resonator. Strikes were performed with an SX 2B wooden drum stick. We show here two synthesized sound results. First a strike at center. This excitation is particularly interesting, because the theoretical eigenfrequencies of the excitation points can be readily found. All eigenfunctions but the zeroth order Bessel function has a zero at the origin. Hence their eigenfrequencies also vanish. This allows us to use the zeros of the zeroth order Bessel function to compare both the simulated sound, the recorded sound to idealized theoretical predictions. The basic result of the measurement is found in Table 1 and a pictorial comparison of the spectra of the simulation and the recording can be seen in Figure 3(a). The recording also shows the following additional peaks comparable to the one associated with the modes of Jo: 4.4, 5.4, 5.8. As can be seen in the spectra of Figure 3(a), the synthesized sound is much cleaner than the physical recording. This has to do with the idealizing assumption of the model equation (1) as well as with the idealized assumption of perfect impulse. The recordings were performed by hand using a physically extended drum stick. Hence there is also a deviation in position possible. Certainly there is much deviation possible due to the author's severe deficiencies as a drummer. The second case is illustrated in Figure 3(b). Here the drum was hit at a one third of the radius away from the center of the drum. In the recording we see that the second partial is much more dominant than in the simulation. This again may be due to inaccurate striking or due to deviation in dissipation mechanisms. In this case, the deter mination of modes to eigenfrequencies is more laborious and has been omitted.

Page  00000004 Simulation _5 Recording Simulation Recording 35C 45 30 3135. 25 25C 30 |15, 20C 25 - 1E50 15< 15C 20 < 0 15 50 10 50 5 L S200 400 600 0 200 0 600 0 0 600 0 200 400 600 Frequency (Hz) Frequency (Hz) Frequency (Hz) Frequency (Hz) (a) Center of the drum. (b) 1/3 of the radius of the drum. Figure 3. Spectra of simulated and recorded drum strikes. The synthesized result sounds different from the recorded sounds. It appears though that this is mostly due to insufficient dissipation models. The recorded sounds have long ringing higher partials, whereas the simulation damps all partials at the same speed. As it is well known that damping is a very important perceptual factor, it remains to be explored how to improve the damping model to make the sound more realistic. 5. CONCLUSION This paper describes first steps towards sound synthesis in the plane using wave fronts as primary numerical objects. This approach preserves many aspects of the structure of the solution, which are lost in methods proposed so far. In particular it retains accurate positions of first arrival fronts up to floating point precision and can easily handle complicated domain shapes. To complete this program, a detailed dissipation model needs to be developed, and the impact of the wake needs to be investigated in more detail. Adaptive ray densities may be helpful to deal with the ongoing divergence of the wavefront. In the long run, however, this may not be a viable synthesis method because the number of rays for reasonable resolution will remain high. The main advantage of this work, is to highlight previously neglected structures of the solution. Maybe ways can be found to exploit this structure in a way to avoid the need for dense sets of rays and hence make this approach practical. In addition ideas presented here can be used to gauge the success of numerical simulations with respect to their ability to resolve geometric features. Acknowledgments Earlier versions of aspects of this work were presented during a Ph.D. short course at University of Verona in September 2004 with the kind invitation of Davide Rocchesso. I thank Renee Hall for acquiring a number of relevant references for me. Eric Blankinship and Erin Panttaja were supportive during early stages of this work. I appreciate the input of three anonymous reviewers. 6. REFERENCES [1] S. D. Bilbao. Wave and Scattering Methods for the Numerical Integration of Partial Differential Equations. PhD thesis, Stanford University, May 2001. [2] C. Cadoz, A. Luciani, and J.-L. Florens. CORDISANIMA: A Modeling and Simulation System for Sound and Image Synthesis - The General Formalism. Computer Music Journal, 17(1): 19-29, Spring 1993. [3] Y. V. Egorov, A. I. Komech, and M. A. Shubin. Elements of the Modern Theory of Partial Differential Equations. Springer, Berlin, 1999. [4] G. Essl. Computation of Wave Fronts on a Disk I: Numerical Experiments. Submitted to Electronic Notes in Theoretical Computer Science, Proceedings of MFCSIT'04, September 2004. [5] G. Essl. The Computational Structure of Waves on Drums for Sound Synthesis. To appear in Proceedings of Forum Acusticum, Budapest, Hungary, Aug 29 - Sep 2 2005. [6] G. Essl. "Whispering" waves and Bate's ridges in a virtual wineglass and elsewhere. Submitted to Acoustics Research Letters Online, January 2005. [7] G. Essl, S. Serafin, P. R. Cook, and J. O. Smith. Theory of Banded Waveguides. Computer Music Journal, 28(1):37 -50, 2004. [8] F. Fontana and D. Rocchesso. Physical Modeling of Membranes for Percussion Instruments. Acustica united with acta acustica, 84(3):529-542, 1998. [9] F. Fontana and D. Rocchesso. Signal-Theoretic Characterization of Waveguide Mesh Geometries for Models of TwoDimensional Wave Propagation in Elastic Media. IEEE Transactions on Speech and Audio Processing, 9(2):152 -161, 2001. [10] K. F. Graff. Wave Motion in Elastic Solids. Dover, New York, 1991. [11] M. Karjalainen and C. Erkut. Digital Waveguides Versus Finite Difference Structures: Equivalence and Mixed Modeling. EURASIP Journal on Applied Signal Processing, 7:978-989, 2004. [12] L. Savioja and V. Vilimaki. Reducing the Dispersion Error in the Digital Waveguide Mesh Using Interpolation and Frequency Warping Techniques. IEEE Transactions on Speech and Audio Processing, 8(2): 184-194, 2000. [13] S. Serafin, P. Huang, and J. O. Smith. The Banded Digital Waveguide Mesh. In Proc. Workshop on Future Directions of Computer Music (Mosart-01), Barcelona, Spain, November 2001. [14] J. O. Smith. Digital Waveguide Modeling of Musical Instruments. Draft of unpublished online manuscript, available at ~jos/waveguide/, 2003. [15] L. Trautmann, S. Petrausch, and R. Rabenstein. Physical Modeling of Drums by Transfer Function Methods. In Proceedings of the International Conference on Acoustics, Speech & Signal Processing (ICASSP), volume 5, pages 3385-3388, Salt Lake City, Utah, May 2001. IEEE. [16] S. A. Van Duyne and J. O. Smith. The 2-D digital waveguide mesh. In Proceedings of the IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA), New Paltz, New York, Oct. 1993. IEEE Press. [17] S. Zambon and D. Rochesso. Space-time simulations of circular membranes by parallel comb filters. In Proceedings of Understanding and Creating Music, Caserta, Italy, December 2004.