# Physical Modeling with the 2-D Digital Waveguide Mesh

Skip other details (including permanent urls, DOI, citation information)This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License. Please contact mpub-help@umich.edu to use this work in a way not covered by the license. :

For more information, read Michigan Publishing's access and usage policy.

Page 40 ï~~Physical Modeling with the 2-D Digital Waveguide Mesh Scott A. Van Duyne savd@ccrma.stanford.edu Julius 0. Smith III jos@ccrma.stanford.edu Center for Computer Research in Music and Acoustics (CCRMA) Stanford University, Stanford, CA 94305 Abstract An extremely efficient method for modeling wave propagation in a membrane is provided by the multidimensional extension of the digital waveguide. The 2-D digital waveguide mesh is constructed out of bidirectional delay units and scattering junctions. We show that it coincides with the standard finite difference approximation scheme for the 2-D wave equation, and we derive the dispersion error. Applications may be found in physical models of drums, soundboards, cymbals, gongs, small-box reverberators, and other acoustic constructs where a one-dimensional model is less desirable. 1 Background Theory There are many musical applications of the onedimensional digital waveguide ranging from the generation of wind and string instrument tones, to flanging effects [Van Duyne and Smith, 1992], to reverberation [Smith, 1987]. We review the theoretical derivation of one-dimensional traveling waves as a basis for development of the twodimensional digital waveguide mesh. 1.1 The 1-D Wave Equation The one-dimensional wave equation for displacement of an ideal vibrating string may be written as follows, utt(t,x) = c2Ux,(t,x), where t is time, x is longitudinal position along the string, u(x,t) is transverse displacement of the string as a function of time and position, utt is the second partial time derivative of u corresponding to the transverse acceleration of a point on the string, and ux is the second partial space derivative of u corresponding the "curvature" of the string at a point. The equation says that the force which accelerates a point on the string back toward its rest position is proportional to how tightly the string is curved at that point [Morse and Ingard, 1968]. It is easy to verify by substitution that this equation is solved by the sum of two arbitrary traveling waves, u(t,x) = g+(x - Ct) + g-(x + ct), where g+(x - ct) represents an arbitrary fixed wave shape traveling to the right, and g-(x + ct) represents an arbitrary fixed wave shape traveling to the left. To see that these waves travel at speed c, note that, as t is increased by I in the expression, g+(x - ct), x must be increased by c for the function argument to remain unchanged. The wave speed is given by c = (T/e)0-5, where T is the constant tension on the string and e is the mass per unit length. Intuitively, we may check that increased tension should speed up wave travel and increased mass should slow it down. 1.2 The Digital Waveguide The traveling wave solution to the one-dimensional wave equation may be implemented digitally with a pair of bi-directional delay lines as shown in Figure 1. The upper rail contains a signal traveling to the right and the lower rail contains a signal traveling to the left. This structure is known as the digital waveguide. Two arbitrary traveling waves propagate independently in their respective left and right directions, while the physical wave amplitude at any point may be obtained by summing the left- and right-going waves. --[ N sample delay N sample delay I-* L I# Figure 1. The Digital Waveguide 1.3 Force, Velocity, and Impedance We need not choose displacement as our wave variable. By taking the time derivative of displacement waves, we may obtain velocity waves. In this case, the physical transverse velocity at a point on the string is the sum of its two traveling components, v = v+ + v-. The transverse force component on a string is proportional to the slope of the string at a given point for small displacements. Therefore, force waves may be obtained by properly scaling the position derivative of the displacement waves to obtain, f =f" + f. See [Smith, 1992] and [Smith, 1993] for a full derivation of all these wave relationships and their digital implementation. 5A.2 40 ICMC Proceedings 1993

Page 41 ï~~Force and velocity are a convenient choice of wave variables as there is a well understood impedance relation between force and velocity in mechanical systems which is analogous to the impedance relation between voltage and current in electrical systems. There is also a wave impedance relation between force and velocity waves on strings which is analogous to the wave impedance relation between voltage and current waves on electrical transmission lines. The acoustical system of the vibrating air column is also mathematically equivalent to both the vibrating string and the electrical transmission line. The wave impedance relationship between the traveling components of force and velocity for the string can be written, = Rv+, where R = (Te)0-5. Intuitively, when a force is applied transversely to a string, the resultant transverse velocity should be slower for greater string mass and also slower for greater string tension. 1.4 The Lossless Scattering Junction It is useful to be able to interconnect waveguides of possibly varying wave impedance at junctions which may be lossless, or which may be loaded with impedances of their own, or be driven by external forces. For example, driving a violin string with a pulsed noise signal representing the bow requires a scattering junction on the string where the bow divides it [Chafe, 1990]. Tone holes in wind instrument models may take advantage of scattering junctions. Strings may be coupled together at a bridge via scattering junctions [Smith, 1993]. Scattering junctions may be used to build up acoustic tubes of varying diameter by joining segments of cylindrical tubes [Cook, 1990]. Julius Smith [1987] has developed a reverberation algorithm which depends on interconnecting any number of varying length and varying impedance waveguides into an arbitrarily elaborate network. The membrane model presented in this paper may be viewed as a canonical form of this reverberation structure. When several strings, say N of them, intersect at a single point, or junction, without loss of energy, we have a "series" junction and require two conditions: (1) that the velocities of all the strings at the junction be equal since they are all moving together at that point, and (2) that the forces exerted by all the strings must balance each other at that point, i.e., they must sum to zero, ft+f2 +... +fN =0. Note that the acoustic tube junction is "parallel" and has the dual constraints, i.e., that the pressures must all be equal and the flows must sum to zero. Figure 2 shows a schematic representation of waveguides intersecting in a lossless scattering junction. The line segments with opposing arrows on them represent the bi-directional delay lines of the digital waveguide shown in Figure 1, with their associated wave impedance, Ri. The circumscribed S represents the junction. Figure. Scattering Junction for N=5 case. Combining the two series junction constraints with the wave variable definitions, vi = v+i + v- and Jh= fli +f-, and with the wave impedance relations, f; = Ri v+i and f-i = -R iv-i, we can derive the lossless scattering equations for the interconnection of several strings, vi -(2 Riv+i) / Ri i i - =Vj - v+ i, where vj represents the junction velocity, the v+i are the incoming waves at the junction, and the v-i are the outgoing waves. These equations say that, as a wave is coming into a junction along a string, some portion of the wave reflects off the junction and travels back where it came from, while the rest of it travels into the junction and is divided among the outgoing waves along the other strings. The relative proportions of this scattering effect is dependent only on the relative impedances of the strings and not on their length. 2 The Two-Dimensional Case 2.1 The 2-D Wave Equation The two-dimensional wave equation for displacement of an ideal membrane may be written as follows, utt(t,x,y) = C2 [ux(t,x,v) + uyy(t,x,y)I, where t is time, x and y are spatial coordinates on the membrane, u(t,x,) is transverse displacement of the ICMC Proceedings 1993 41 5A.2

Page 42 ï~~membrane as a function of time and spatial position [Morse and Ingard, 1968]. In the one-dimensional string case, we could solve and implement the wave equation as two bidirectional traveling waves. In the 2-D membrane case, the traveling wave solution involves the integral sum of an infinite number of arbitrary plane waves traveling in all directions, u(t,x, y) = Jgc(x cost + y sina- ct) dx Since assigning one waveguide to each of the infinite plane waves is not feasible, we need an alternative approach. 2.2 The 2-D Digital Waveguide Mesh Proposed in this paper is a formulation of the 2 -dimensional wave equation in terms of a network of bi-directional delay elements and 4-port scattering junctions. This structure can be viewed as a layer of parallel vertical waveguides superimposed on a layer of parallel horizontal waveguides intersecting each other at 4-port scattering junctions between each bidirectional delay unit. Figure 3 shows such a mesh. In the canonical case, the scattering junctions are taken to be equal impedance lossless junctions and the interconnecting waveguides are of unit length. velocity waves traveling in the two-port, bidirectional, delay units. On the other hand, if we view the mesh as a lattice of interconnected acoustic tubes, the pressures at each junction must be equal, and the flows into each junction must sum to zero; in this case, we have parallel scattering junctions with pressure or volume velocity waves traveling through the delay units. 3 Empirical Analysis It is evident that, given an initial excitation at some point on the digital waveguide mesh, that energy from that excitation will tend to spread out from the excitation point more and more as the traveling waves scatter through the junctions. It is not, however, easy to see that the wave propagation on the mesh converges to that on the ideal membrane. Figure 4. Wave Propagation on the Mesh 3.1 Animation of the Mesh A visual verification of the waveguide mesh algorithm can be seen in Figure 4, which shows three separate time frames of an animation computed directly from the algorithm. The top frame shows the initial deflection loaded into the mesh. Each intersecting grid point represents a scattering junction. The next two frames show the circular propagation outward of the initial excitation in a way consistent wave propagation on the ideal membrane. Figure 3. The 2-D Digital Waveguide Mesh If we view the mesh as a lattice of interconnected vibrating strings, the displacement velocities at the four ports of each junction must be equal, and the forces at each junction must sum to zero; in this case, we have series scattering junctions with force or 5A.2 42 ICMC Proceedings 1993

Page 43 ï~~3.2 Sounds from the Mesh As another check of the mesh algorithm, we can compare the expected modal frequencies on an ideal membrane with those generated from the mesh model. The allowed frequencies in a theoretical ideal square membrane with clamped edges are proportional to (m2 + n2)0.5, for m = 1,2,..., and n = 1,2,...[Morse and Ingard, 1968]. These modes may be labeled (m,n) for any given m and n. In Table 1 is computed a list of the normalized frequencies of the first few of these modes given as multiples of the lowest allowed frequency. Table 1. Modes on Ideal Square Membrane (1,1)-- 1.00 (1.5) -- 3.60 (1,2)-3 1.58 (2.5) -- 3.80 (2,2) -- 2.00 (4.4) --- 4.00 (1,3) - 2.24 (3,5) -- 4.12 (2,3) - 2.55 (1,6) -- 4.30 (1,4) - 2.92 (2.6) --- 4.47 (3,3) -43.00 (4,5) -- 4.50 (2,4) -- 3.16 (3.6) --- 4.74 (3,4) - 3.54 (5,5) --, 5.00 Figure 5 is a spectral analysis of a sound generated by a square 10 junction by 10 junction digital waveguide mesh reflectively terminated at the boundaries. A careful inspection of the plot will reveal that the theoretical modal frequencies listed in Table I are all present and accounted for in the sound generated by the model. This would indicate that the mesh is doing the right thing. Figure 5. Measured Modes on the Square Mesh Frequencies up to half the sampling rate are shown in Figure 5. However, notice that the spectrum mirrors around one quarter of the sampling rate. This symmetry, which also occurs in the one dimensional waveguide case, is a result of the fact that, when the waveguide or waveguide mesh is reflectively terminated, all the unit delays on the upper rails can be commuted down to the lower rails making the system a function of z"2, in effect, over-sampling the system by a factor of two. 4 Mathematical Analysis Given the the function f(x), one may approximate its first derivative by the difference, f '(x) = ftx + Ax)-.fx) Ax By applying this expression to itself we may arrive at the standard difference scheme approximation for the second derivative, f"(x) Ax + Ax) - 2Aix) +f(x - Ax) Ax2 Computable and numerically stable difference schemes can be found for many partial differential equations by substituting approximations of this kind into the equation. The digital waveguide mesh algorithm, in fact, may be interpreted as a difference scheme for computation of the two-dimensional wave equation. 4.1 The Mesh as a Difference Scheme In the digital waveguide mesh, we require that the impedances in all directions be equal, as we would desire for the isotropic membrane case. Setting the wave impedances of the four unit waveguides attached to each junction point in the mesh equal, i.e., R1 = R2 = R3 = R4, the scattering equations then reduce to, V+ I mn(n)+v+2l nm(n)+v+3l m(n)+v+4l m(n) vj,,(n) = 2 V'it~1( n) = vjim(n) - v+itm(n), where the l,m indices represent the spatial position of the junction in the mesh and the n index represents the current time sample. vj,.n(n) represents the velocity of the junction at position l,m at time n; v+il.m(n) and vil,1(n) represent the four input and output waves to that junction, respectively. As an implementation note, observe that this junction computation may be performed with 7 adds (or subtracts), 1 shift (to divide by 2), and no multiplies. In addition to these scattering equations, we may also note that the sum of the inputs to a junction equals the sum of the outputs, I v~i,~1(n)= X.,,(n), i i and that, since the junctions are interconnected by unit delay elements, the input at one port of one junction is equal to the output at the opposing port of the adjacent junction at the previous time sample, v+ ops(n)= v-o.st(n-). With a little perseverance, one may manipulate all these relations algebraically into the following difference equation, Vjl~1(n) - 2Vj, (n-!l) + Vjl,,1(n-2) = O.5 [v J/.,1 (n-I1) - 2Vjl, n(n-l1) + Vjl, 1,I(n-I) ] + O.5[vjl1l~(n-l) - 2vj,~,(n-l) + vlimnl] ICMC Proceedings 1993 43 5A.2

Page 44 ï~~A comparison of this difference scheme with the two-dimensional wave equation, utt(t,x,y) = c2(uxx{(t,x, y) + u,(t x,y)), reveals that it is the standard second-order difference scheme for the hyperbolic partial differential wave equation for the ideal membrane, with wave propagation speed c = 2-.5 --0.7, and the time and spatial sampling intervals (X=Ax, Y=Av, T=At) taken to be equal to each other. The mesh implements a wave propagation speed of one-half unit diagonal distance per time sample. Intuitively, when we superimposed the perpendicular layer of parallel waveguide strings to form the mesh, we doubled the mass density per unit area, thereby reducing the wave speed by one over the square root of two. Defining, X = TX-' = 1, we observe that the CourantFriedrichs-Lewy stability condition, lcAI 2-0-5, is satisfied by this difference scheme [Strikwerda, 1989]. This condition says that for a difference scheme to track the solution of a hyperbolic equation with two space dimensions, the cone of dependence for each point of the continuous solution must lie within the pyramid of dependence for each point of the difference scheme solution. Since the condition is satisfied in the equality, the lowest possible dissipation and dispersion error for this particular scheme is obtained. The numerical approximation schemes for initial value problems involving second-order hyperbolic partial difference equations usually require a multistep time scheme which retains values for at least two previous time frames. This is to cope with the second partial time derivative in the equation. The waveguide mesh reduces this structure to a one-step time scheme where each new time frame may be computed wholly from the previous time frame. This is made possible by the use of traveling wave components in place of physical wave variables. 4.2 Von Neumann Error Analysis Von Neumann analysis of finite difference scheme approximations of partial differential equations uses Fourier transform theory to compare the evolution over time of the spatial spectrum in the continuous time solution to that in the discrete time approximation [Strikwerda, 1989]. Recall that to solve an ordinary linear differential equation, we may reduce the problem to a polynomial in s by taking the Laplace transform and replacing orders of derivatives with powers of s. Similarly, we can take a spatial transform of a partial differential equation with independent time and space variables to obtain an ordinary differential equation describing the evolution of spatial spectra over time. From here, we can check how the spatial spectrum evolves after a time delay of T seconds. The ratio of the spatial spectrum at time t + T to the spatial spectrum at time t is known as the spectral amplifi cation factor. Recall that in the discrete case, to solve a timeindexed difference equation, we may reduce the problem to a polynomial in z'1 by taking the Z - transform and replacing samples of time delay in the index with powers of z". If we have a time- and space-indexed difference scheme approximation to compare with our partial differential equation, we may perform a similar spectral evolution analysis on the difference equation using discrete Fourier transforms to obtain a discrete spectral amplification factor. In the case of the ideal membrane equation, we already know the solution is an integral sum of plane waves all moving at constant speed c. The speed of wave travel is independent of spatial frequency, i.e., smooth low spatial frequency waves travel at the same speed as jagged or ripply-looking high spatial frequency waves. This means that there is no frequency dispersion in the ideal membrane. After T seconds, the position of a plane wave of spatial frequency 4, traveling in direction a, would have moved forward a distance of cT, corresponding to a spatial phase shift of --ctT. The spectral amplification factor would then be e-CT. In the difference scheme derived in Section 4.1, computation of a discrete spectral amplification factor is a little messier. Unfortunately, the speed of plane wave travel on the digital waveguide mesh is dependent on both the direction of travel, a, and on the spatial frequency of the wave, 4. The discrete spectral amplification factor may be written in the form e-JC(t),T, where c'(at,4) represents the direction and frequency dependent speed of plane wave travel on the waveguide mesh. When making two-dimensional spatial transforms, we take x *- 4 and y H. The coordinate frequencies, t and 2 are hard to understand conceptually, but viewing the transform in polar coordinates, we see that the point ( l in the twodimensional frequency space is referring to a plane wave of spatial frequency = ( _2 + 22)0.5, oriented in the radial direction, at = tan~l1. Using the procedure outlined above, a closed-form 5A.2 44 ICMC Proceedings 1993

Page 45 ï~~expression for the normalized speed of plane wave travel in the waveguide mesh may be found, c'( l, ) / c = 20.5 ( T)-1 tan-1 (4 - b2)05 b-1, where b = cos 41T + cos 4.T. 2 - 2T Figure 6. Wave Travel Speed vs. Frequency Figure 6 shows a plot of the normalized wave travel speed on the mesh versus spatial frequency. The center region of the plot corresponds to low spatial frequencies; the outer regions of the plot correspond to higher spatial frequencies. The angular position on the plot, as seen from the frequency plane origin, corresponds exactly to the direction of plane wave travel on the mesh. Notice that near the center of the plot, corresponding to smooth, low frequency plane waves, the c'I c ratio is fairly close to 1. Also, c/c = I exactly along the diagonals, corresponding to no dispersion at any spatial frequencies when traveling in a direction diagonal to the mesh coordinate system. In waves traveling along the coordinate axes of the mesh, we see a fall off in travel speed in the higher spatial frequencies. A Figure 7 shows three time frames of the mesh initialized with a deflection containing high spatial frequencies. Notice how the wavefront smooths out along the mesh coordinate directions, corresponding to high frequency dispersion, while it remains sharp along the diagonal directions, corresponding to no dispersion. In a bounded mesh, speed distortion results in a mistuning of resonant modes. This distortion can be reduced by allpass filtering and/or warping of the membrane boundary in a compensating manner. Oversampling the mesh and low-passing can eliminate the effect to arbitrary accuracy. We note that the high frequency modes of a membrane become so dense that, in audio contexts, this error may not be important. 5 Implementation Features The digital waveguide mesh may be computed in parallel, and without multiplies. In addition numerical round off loss may be redistributed back into the mesh to create a zero-loss system. 5.1 Two Pass Parallel Computation The network elements in the 2-D digital waveguide mesh are of two types: 4-port scattering junctions and 2-port bi-directional unit delays. If the unit delays are double buffered, so that each delay has its own input and output buffers, the computation of all the elements in the mesh can be segregated and computed in any arbitrary order or in parallel, according to the following two pass computation scheme: (1) The scattering junction outputs are computed from their known inputs and placed at the junction outputs. This constitutes the scattering pass. (2) The outputs from each scattering junction are placed at the inputs of the adjacent scattering junctions, thereby implementing the bi-directional delay units. This constitutes the delay pass. Due to the possibility of arbitrary ordering of the scattering computations, implementation on a parallel computing architecture with local four-sided connectivity between processors is ideal for the mesh algorithm. In this implementation, the junction equations are computed in the processors; and then the data transfer cycle is used to transfer data from the outputs of each processor to the appropriate inputs of the adjacent processors. Since the equal impedance 4-port lossless scattering junction is multiply-free, as pointed out in Section Figure 7. Wavefront Dispersion on Mesh ICMC Proceedings 1993 45 5A.2

Page 46 ï~~4.2, a VLSI implementation may be constructed with a handful of gates with no need for hardware multipliers. Since the junctions may be computed in parallel, the whole mesh may be computed in the time it takes to do 7 adds (i.e., 3 adds and 4 subtracts) and one shift. In fact, the four subtracts may be performed in parallel. 5.2 Energy Preserving Junctions When performing the multiply-free junction computation, one divide by two is required. If a simple sign-preserving right shift and truncation is used for this operation, the junction value is rounded toward zero in the case of positive numbers and rounded away from zero in the case of negative numbers. This is a round down in both cases, which could introduce a negative offset into the values of the mesh which may eventually lead to numerical instability or reduced dynamic range, if there is no loss in the system somewhere else. The usual solution to such a problem would be to make a conservative rounding toward zero in both the positive and negative cases. This way no energy and no DC drift will be introduced into the system. This method is known to work quite well in onedimensional feedback loops. Unfortunately, in the two-dimensional mesh case, there are so many junctions that the cumulative losses in all the junctions add up to a noticeable amount. An energy preserving method of junction computation may be constructed as follows. When shifting a binary number to the right, there are exactly two cases: (1) the low-order bit which is shifted off the end of the word is zero, and the computation is exact, or (2) the low-order bit is one, and the error is exactly 0.5. When the junction inputs are subtracted off vj in computing the scattering junction outputs, this error is magnified by 4 and the 0.5 error propagates into the four output signals equally. Note that the error is in the same direction in all four cases so the total error is 4. 0.5 = 2 full bits. To preserve energy in the mesh, round two output signals down (i.e., just truncate) and round the other two output signals up (i.e., add the low-order bit back in after truncating). This re-distribution of the error produces a numerically exact lossless scattering junction. in effect, the slight numerical error has been converted into a slight scattering dispersion error. Whereas the numerical error was problematic, the dispersion adjustment is vanishingly small. 6 Modeling with the Mesh To build a model of a drum membrane, we need to clamp down the boundary of the mesh, corresponding to terminating the mesh reflectively, inverting the traveling waves at the edges. Since the 2-D digital waveguide mesh is just a big linear system, filters representing loss in the system are easily interconnected, and may be consolidated as desired around the rim due to commutativity of linear systems. Modeling a stiff plate might be accomplished by letting the waves reflect off the boundary with out inversion (for an unclamped plate). To help out with the greater spacing of higher modes in the plate caused by stiffness, some appropriate allpass filtering might be introduced. Figure 8 suggests the possibility of modeling a harpsichord by connecting an array of waveguide plucked string models to a waveguide mesh representing the sounding board via 5-port junctions. The soundboard mesh would have appropriate boundary filters with low pass characteristics to represent loss and allpass characteristics to represent stiffness effects in the board. In this model the mesh is used as a resonant coupling connection which both reverberates the string outputs and scatters energy into strings which have not been plucked, thereby inducing sympathetic vibrations. Z1 NJ Mn N, M~z1) ZNF \ - -1 Z-N \._. o(n) Figure 8. The 8-String Harpsichord. 6.1 Extensions of the Mesh The mesh algorithm assumes nothing about its boundary conditions. There is nothing to stop one from connecting one edge of the mesh to the opposite edge to produce a cylindrical topology, as shown in Figure 9. Furthermore, it is straight forward to extend the algorithm into three dimensions by layering several 2-D meshes above each other, replacing all the 4-port junctions with 6-ports and 5A.2 46 ICMC Proceedings 1993

Page 47 ï~~connecting up the layers. This topology is shown in Figure 10. Figure 9. The Cylindrical Mesh Figure 10. The 3-D Digital Waveguide Mesh cylindrical mesh modeling the sides of the drum and 3-D mesh inside the drum modeling the air cushion. Figure 12 shows a guitar model with 2-D mesh material modeling the bridge and body shell and 3-D mesh material modeling the resonant body cavity. 7 Summary Although finite element and difference scheme approximation methods are known which can help with the numerical solution of the 2-D wave equation, these methods have two drawbacks: (1) their heavy computational time is orders of magnitude beyond reach of real time, and (2) traditional problem formulations fit poorly into the physical model arena of linear systems, filters, and network interactions. On the other hand, the 2-D digital waveguide mesh formulation proposed in this paper, while corresponding exactly with the standard difference scheme approximation, may be implemented in a fully parallel, multiply-free formulation; the energy preserving, digitally exact round-off method eliminates numerical problems; the mesh extends simply to 3- or N-dimensions; and, finally, the algorithm is a linear network which connects up easily to other physical models. References [Chafe, 1990] Chris Chafe. "Pulsed Noise in SelfSustained Oscillations of Musical Instruments," Proceedings IEEE International Conf. Acoust. Speech and Signal Proc., Albuquerque, NM. [Cook, 1990] Perry Cook. Identification of Control Parameters in an Articulatorv Vocal Tract Model, with Applications to the Synthesis of Singing. Ph.D. Diss., Elec. Eng. Dept., Stanford University. [Morse and Ingard, 1968] P. M. Morse and K. U. Ingard. Theoretical Acoustics. McGraw-Hill, New York. [Smith, 1993] Julius O. Smith. "Efficient Synthesis of Stringed Musical Instruments," Proceedings ICMC, Tokyo. [Smith, 1992] Julius 0. Smith. "Physical Modeling Using Digital Waveguides," Computer Music Journal, 16:4. [Smith, 1987] Julius O. Smith. Music Applications of Digital Waveguides. CCRMA, Stanford Univ., Stanford, CA, Tech. Rep. STAN-M-39. [Strikwerda, 1989] J. Strikwerda. Finite Difference Schemes and Partial Differential Equations. Wadsworth & Brooks, Pacific Grove, CA. [Van Duyne and Smith, 1992] 5. A. Van Duyne and J. 0. Smith. "Implementation of a Variable Pick-Up Point on a Waveguide String Model with FM/AM Applications," Proceedings ICMC, San Jose, CA. Figure 11. The 3-D Drum Model. M2 Figure 12. The 3-D, One-String Guitar. With constructs such as these, fully physical models of musical instruments can me made. For example, Figure 11 shows the 3-D drum model with a mesh modeling the drum head, connected to a stiff ICMC Proceedings 1993 47 5A.2