Page  411 ï~~Traveling Wave Implementation of A Lossless Mode-Coupling Filter and The Wave Digital Hammer Scott A. Van Duyne CCRMA, Stanford University savd(ccrma. stanford. edu John R. Pierce CCRMA, Stanford University jrp~ccrma. stanford. edu Julius 0. Smith, III CCRMA, Stanford University jos(ccrma. stanford. edu Abstract Recent work has led to physical models of strings and acoustic tubes based on the digital waveguide, and to physical models of membranes and plates based on the 2D digital waveguide mesh. This paper introduces two new developments for these models which extend their usefulness: first, a mode-coupling filter structure based on a passive nonlinear impedance circuit; and second, a physical model for a piano hammer, or felt mallet, based on the wave decomposition of a nonlinear mass and spring system. 1 Motivation Since the development of the 2D digital waveguide mesh to model wave propagation in membranes and plates [Van Duyne & Smith, 1993], it has been evident that a physically accurate excitation means was necessary, as well as a way to convert linear plates into gongs and cymbals. 1.1 Memoryless Nonlinearities Nonlinearities, small or large, favorably affect the sounds of many musical instruments. In gongs and cymbals, nonlinearities cause the transfer of energy from lower frequency modes to higher frequency modes after the instrument has been struck. These nonlinearities do not generate new energy, only transfer it. While memoryless square-law and look-up table nonlinearities may be incorporated in computer generation of sounds, these means often cause system energy loss or gain, and are difficult to control when a range of large and small effects are desired. Almost all memoryless nonlinearities involving only one input signal may be put in a form where the nonlinear effect is handled by a look-up table, where the output signal is taken from a table indexed by the input signal. If the table contains a simple line of slope 1 (45 degrees), the output equals the input. The nonlinear effect is produced in accordance with how much the table deviates from a line of slope 1. While many spectral modifications can be achieved by memoryless nonlinearities [Rodet, 1993], energy conservation is out of control without additional amplitude tracking and/or scaling elsewhere in the loop. If the slope of the nonlinear look-up table near index zero is less than 1, low amplitude signals are driven to zero in a feedback loop. If the slope of the table near index zero is greater that 1, then low amplitude signals are driven to some higher and higher amplitude until the indexing signal values reach a point in the table where table values drop below the line of slope 1. In addition, the greater the nonlinear effect, i. e., the further the table deviates from a line of slope 1, the more troublesome the energy control becomes. Yet, for gongs and cymbals, a large nonlinear effect is required. To resolve this dilemma, we have approached nonlinearity through physical modeling. We have derived a nonlinear filter based on a real physical system constructed from passive, lossless elements only, namely two springs of differing stiffness. When such filters are connected to a 2D digital waveguide mesh membrane model, a natural, controllable spreading of spectral energy occurs, and fine gong sounds may be produced. 1.2 Striking Musical Instruments The musical effect of a struck string or drum is greatly determined by the nature of the hammer or mallet which strikes it. The attack transient of a struck instrument can be approximated by the injection of an appropriate, pre-computed excitation signal into the resonant system. However, this excitation method is not sufficient to cope with the complexities of certain real musical situations. For ICMC Proceedings 1994 411 Sound SynthesisTechniques

Page  412 ï~~example, when a mallet strikes an ideal membrane, it sinks down into it, feeling a pure resistive impedance. The depression induces a circular traveling wave outward. If the membrane were infinite, the waves would never return, and the mallet would come to rest, losing all its energy into the membrane. If the membrane is bounded, however, reflected waves return to the strike point to throw the mallet away from the membrane. The first reflected wave to reach the mallet may not be sufficiently powerful to throw the mallet all the way dear, or may only slow down its motion; and later reflected waves may finally provide the energy to finish the job. This complex mallet/membrane interaction can have very different and difficult to predict acoustical effects, particularly when a second or third strike occurs while the membrane is still in motion, such as in a drum roll. We have obviated these problems of system excitation by constructing a physically accurate model of the piano hammer, or felt mallet, which can be coupled into a waveguide string or membrane model. As a result, no complex set of excitation signals need be designed and stored. Only the correct choice of physical hammer parameters must be specified: hammer mass, felt stiffness response, felt hysteresis loss, and so on. The complex interaction between hammer and string, or between mallet and membrane, falls out of the model naturally. 2 Waveguide Review As a theoretical introduction to these two additions to the physical modeling repertoire, we review the traveling wave approach used to model the distributed impedance of the ideal string. The one-dimensional wave equation may be solved as the sum of two waves traveling in opposite directions [Morse & Ingard, 1968]. This solution has led to the modeling of strings and acoustic tubes as a pair of bi-directional delay lines, each delay line corresponding to one of the traveling waves. This structure is known as the digital waveguide [Smith, 1992]. The traveling waves may represent displacement, velocity, slope, force, or other physical variables. Traveling waves are mathematical constructs, the actual physical values at any point on the string being the sum of the right- and leftgoing traveling components at that point. Figure 1 illustrates a physical variable, u, as the sum of its traveling components, u+ and u-, at some point on the wavegulde. There is a wave impedance relationship between the traveling components of force and velocity in the string, corresponding to the wave impedance relation between pressure and flow in acoustic tubes, or KfDELAY UNE + a DELAY UPIE Figure 1: The Digital Waveguide voltage and current in electrical transmission lines, which can be written, f+=Rov+ and f-=-Rov (1) where fÂ~ and vÂ~ are the traveling wave components of force and velocity respectively, and where the wave impedance is R0 = VK, "K being the constant tension on the string and E being the mass density per unit length. Intuitively, these equations say that when a force is applied transversely to a string, the resultant transverse velocity should be less for greater string mass and for greater string tension. The change in sign for the left-going components is due to coordinate system choice. Note that whether we have modeled velocity or force waves in our system, we may compute both physical force and physical velocity at any point along the string from the resident traveling wave components using Equation (1) to make the change of variables, V = V+ + v- = (1/Ro)(f+ - f-) f =fI+f-o(v,-) (2) (3) 3 A Lossless Mode Coupler 3.1 The Physical Impedance Circuit Consider a string terminated by a double spring apparatus as shown in Figure 2. Three states of the system are shown in the figure: First, the lower spring is compressed, while the upper spring is at rest; second, both springs are at rest; and, third, the upper spring is compressed, while the lower spring is at rest. In effect, the spring termination apparatus is equivalent to a single nonlinear spring whose stiffness constant is k1 when the displacement is positive and k2 when the displacement is negative. Now consider what is happening to the energy in the system. When the lower spring is compressed, some energy from the string is converted to potential energy stored in the spring. When the lower spring returns to its rest state, the stored spring energy is entirely returned to the string, and the spring contains no stored energy. When the upper spring is then compressed, exactly the same kind of Sound SynthesisTechniques 412 ICMC Proceedings 1994

Page  413 ï~~1 at rest k compressed k2 attrest k at rest k compressed k1 at rest kk Figure 2: String Terminated with Double Spring energy exchange occurs. This ideal system is physically passive and lossless since energy is neither created nor destroyed. If the spring stiffness constant had changed while stored potential energy was still in the spring, i. e., when one of the springs was still compressed, the stored energy would be scaled by the new relative stiffness of the spring. In this case, the stored energy before the stiffness change would be different from the stored energy after the stiffness change, leading to the creation or loss of energy, possibly resulting in a non-passive system. In our model of this physical system, we must be careful to change the spring stiffness coefficient at the right time, to preserve passivity and losslessness. Passivity is the requirement that no energy be created by the system. When energy is created in a feedback loop, stability problems may ensue. We have specifically tried to discover a nonlinear system which is passive, and which is lossless, so the system loss may be decoupled from the nonlinear effect, and designed separately. 3.2 Deriving the Nonlinear Filter We will first derive the linear spring equations, then terminate a string with a linear spring. From there we determine how to vary the spring stiffness appropriately, and investigate the way in which energy spreads through the modes. Figure 3: Simple Linear Spring System 3.2.1 The Simple Linear Spring The force equation for the ideal linear spring shown in Figure 3 is, f(t) = kx(t) ==. df(t)/dt = kv(t) (4) where f(t) is the force applied on the spring, x(t) is the compression distance of the spring, v(t) is the velocity of compression, and k is the spring stiffness constant. Taking the Laplace transform, and assuming no initial force on the spring, we get, F(s) = (k/s)V(s) (5) Here, k/s is the lumped impedance of the spring. Setting s = jw will give the frequency response of this system. 3.2.2 Terminating a String with a Spring Wave Impedance R0 fr k Figure 4: String Terminated by Simple Spring Figure 4 shows a string terminated by a spring. As was suggested in Equations (2) and (3), the physical force at the string termination is the sum of the transverse force waves on the string at that point, f = f,- + fl, while the physical velocity at the termination is the difference between the force waves scaled by 1/Ro, v = (1/Ro)(fr - fl). We may, therefore, re-formulate Equation (5) as a transfer function from Fr to F, F(s) = (k/s)V(s) F r(s) + F(s) = (k/s) F,(s) -F(s) Fo F (e) - kAs - RO F(s) F()=k/s +R (6) (7) (8) The force wave transfer function is stable allpass since its pole is at s = -k/Ro and its zero is s = k/Ro, where k and R0 are defined to be positive real numbers. To move from the continuous physical system to the discrete digital filter, we use the conformal bilinear transform from the s-plane to the z-plane [Nehari, 1952], 1 +-z-I S 1+ z-1 (9) The bilinear transform maps DC in the continuous system to DC in the digital system, while mapping infinite frequency in the continuous system to ICMC Proceedings 1994 413 Sound SynthesisTechniques

Page  414 ï~~frn) f (n) Figure 5: The Allpass Spring Termination half the sampling rate, or the Nyquist frequency, in the digital system. a is a degree of freedom which may be used to control the frequency warping. It is usual to choose a = 2/T to obtain faithful frequency response at the low end of the frequency range. We apply the bilinear transform to Equation (8) to obtain, Fi(z) = H(z)Fr(z) (10) where, ao +z-' _k- aR0 H(z) = ao z and ao=k+ao (11) Hz=1+aoz- Ro ao ranges from -1 to 1 as k ranges from 0 to oo. Figure 5 shows an implementation of this spring termination filter. 3.2.3 Making the Spring Nonlinear Since the filter coefficient, ao, represents the spring stiffness, we need only change ao to effect a change in the stiffness of the spring termination. We must, however, effect the change just at the point where the spring is at rest to preserve physical and digital passivity and losslessness. Consider the internal filter signal u(n) in the alpass filter implementation of H(z) shown in Figure 5. From the diagram, we may observe, is closest to zero. This is the physically correct time to let the spring stiffness coefficient change to model the nonlinear spring termination system shown in Figure 2. Digital Interpretation The internal filter state signal, u(n), represents the energy stored in the filter. Assuming no input signal, i. e., fr(n) 0=0, and some internal state value at time n = -1, namely, u(-1) -- uo, we may then compute the output of the filter caused only by this internal energy state value as, fg(n) = (1 - ao2)(-ao)"uo (16) If we change the filter coefficient, a0, it is clear that the internal state energy will ring out of the filter with a different decay rate than if the coefficient had not been changed. Such coefficient changes, if made arbitrarily, may lead to instability in a feedback loop. However, if uo is zero or near zero, we change the coefficient with relative impunity, since the resultant discontinuity in the state energy will be minimal or zero. Therefore, we choose to gate the filter coefficient change on the sign of u(n), to maintain passivity in the nonlinear allpass filter. This method of gating the coefficient change is both physically correct and digitally correct. Implementation Figure 6 shows a digital system diagram for the nonlinear string/spring system of Figure 2. Note that, for convenience, we use a look-up table notation to show how u(n) gates the coefficient change; but, in general, this would be implemented in simple conditional logic. fi(n) = aou(n) + u(n - 1) u(n) = fr(n) - aou(n - 1) (12) (13) Physical Interpretation The actual physical force applied to the spring termination is equal to the sum of the input and output force waves, as defined in Equation (3). We may derive an expression of the actual force on the spring, f(n), from (12) and (13), f(n) - fr(n)+ft(n) (14) = (1 + ao)[u(n)-+u(n-1)J (15) Equation (15) says that the actual physical force on the spring is proportional to a linearly interpolated value of signal u at time n -0.5. From Equation (4), displacement of the spring termination is zero when force is zero, and f(n) is zero when u(n) + u(n - 1) is zero. Therefore, when u changes sign between times n-1 and n, the spring displacement Figure 6: Nonlinear String/Spring System Figure 13 in Section 4.2.3 shows how a gong model may be constructed by attaching a series of these passive nonlinear fiters (marked "PNF" in the diagram) to the rim of a 2D digital waveguide mesh. 3.3 Phase Modulation Explanation The time-varying allpass termination filter shown in Figure 6 is essentially lossless and passive. However, for it to be useful it must behave in the desired manner. That is, energy in existing resonant modes of the main system must be caused to spread locally in the spectrum to nearby modes of the system, as Sound SynthesisTechniques 414 ICMC Proceedings 1994

Page  415 ï~~is observed in real musical instruments. This termination filter is somewhat hard to analyze strictly, due to its signal dependent time-variation and due to its inclusion in a feedback loop system. But we may obtain an intuitive understanding of its operation by considering a simpler, but non-passive, form. We have observed that the termination filter in Figure 6 is a one-pole allpass filter with a time-varying coefficient. Consider the frequency response of an allpass filter of the same form with a sinusoidally varying coefficient taking on values between a3 and a2. Intuitively, this filter will be performing a phase modulation on the input signal. We should expect the output of the filter to contain sidebands generated by this phase modulation at multiples of the modulation frequency. In fact, Figure 7 shows the dB magnitude spectrum of the output signal of this filter with an input sine wave "center frequency" of 8000Hz and a coefficient "modulating frequency" of 2000Hz. The resultant sidebands of the output signal are as expected. 75 50 25 -25 -50 -75 I _ - 500 10 00 15000 20000 Figure 8: Evolving Spectrum from Nonlinear Loop a sideband coincides with a supported mode, that mode will be driven by the energy from the appropriate sideband. Energy from sidebands which do not fall on supported modes will not drive any particular mode and will simply be absorbed back into the system. Since the passive nonlinear filter produces sum and difference frequencies of the input signal, we can expect that at least some of the main system modes will be hit and energy spreading will occur. Figure 8 shows a spectrogram of the gradual energy spreading of a waveguide string system terminated with a passive nonlinear filter such as that shown in Figure 6. The rate of spreading, and the spectral region where it is most active, may be controlled by the choice of a1 and a2. 4 The Wave Digital Hammer 4.1 The Physical Hammer Model We follow [Chaigne & Askenfelt, 1994] in viewing the piano hammer, or felt mallet, as a nonlinear mass and spring system, the spring representing the felt portion of the hammer or mallet. Since the felt is very compliant when the mallet is just barely touching the membrane, yet very stiff when fully compressed, we must use a nonlinear spring in the model, whose stiffness "constant" varies with felt compression. 4.2 Derivation of the Digital System In deriving the wave digital hammer filter, we first make a wave decomposition of the linear mass and spring system, then attach it to a waveguide string or membrane model via a scattering junction, find a method to vary felt stiffness with compression, and finally account for hysteresis loss in the felt. 4.2.1 The Mass and Spring System First we consider a linear mass and spring system, as shown in Figure 9. We wish to define the equations of motion for a mass and spring system driven at the spring end by an external force, where the mass has an initial velocity at time zero. Since the mass and spring are in parallel, we may observe that the driving point force applied to the left side Figure 7: Output dB Magnitude Spectrum from Sinusoidally Driven Allpass Filter with Sinusoidally Modulated Coefficient If we now modulate the coefficient, ao(n), with a square wave, instead of a sine wave, the output signal spectrum of the filter will contain greater emphasis in the odd sidebands than a simple sinusoidally modulated filter, due to the odd harmonics in the square wave coefficient modulation signal. In the passive nonlinear filter form of Figure 6, where the coefficient change is controlled by the change in sign of the internal signal, u(n), we have a modulating signal which is quasi-square wave, in that it flips between two distinct values, and it has a "fundamental" frequency related to the input signal to the filter, since u(n) is just a highpassed version of the input signal, fr(n), as indicated in Equation (13). The filter output signal will therefore contain sidebands corresponding to sum and difference frequencies of the input signal, generally grouping near frequencies in the input signal. For coupling to occur in the resonant system, the sidebands produced by the modulated allpass filter must fall on supported modes of the system. When ICMC Proceedings 1994 415 Sound SynthesisTechniques

Page  416 ï~~k Figure 9: Mass and Spring Force Diagram of the spring, f, is equally applied to the spring and to the mass, f =fk=fm (17) Also, the driving point velocity, v, is equal to the sum of the spring compression velocity and the mass velocity, V=Vk+Vm (18) We have already discussed the spring system, and can find an expression for Vk (s) by rewriting Equation (5), vh-(n) = v(n) vh+(n) =vown) LThz va6(n) Figure 10: The Wave Digital Mass and Spring between the force and velocity waves traveling in the same direction determined by the parameters of the medium in which the waves are traveling. In the lumped case, the waves may only be understood to be traveling instantaneously; and the waves are not traveling in any particular medium. However, we may still define an arbitrary reference impedance, Rh, and obtain a wave impedance relation between force and velocity waves, Vk(s) = (s/k)Fk(s) (19) Fin= RhVn Fois = -RhVout (24) We may find an expression for Vm (s) by taking the Laplace transform of, Note that Rh is truly arbitrary. We will make good use of this fact later on. Making the change of variables from F and V to 14,, and Vout, we obtain a velocity wave transfer function from Equation (23), dVm(t) mvo3(t) fm (t) = rndt (20) where mvo3(t) represents an initial momentum impulse setting the mass in motion at velocity v0. Combining Equations (17), (18), and (19) with the Laplace transform of Equation (20), we compute the driving point admittance relation for the mass and spring system in Figure 9, V(s) = Vk(s) + Vm () (21) Vo, = s2 - (k/Rh)S + k/rnV s2 + (k/Rh)s +k/mrn + k/Rh + 2 + (k/nh)S + k/m (25) To move to the digital domain, we apply the bilinear transform (9), and band limit the continuous time impulse function mapping 6(t) -+ 6(n)/T, to obtain a second-order digital transfer function. With a careful choice of the arbitrary reference impedance, 8 1 V0 -F k (s ) + - F ( s ) + - k ms s =_(s2+klm) ks 8 (22) (23) mka Rh = ma2 + k (26) The admittance function, (s2 + k/m)/ks, is the steady state response of this second-order system. Note the admittance zeros are at Â~j /k/fm, and the system oscillates at that frequency. The vo/s term represents the transient effect of the time zero momentum impulse on the mass. 4.2.2 Traveling Wave Decomposition of the Lumped System It is possible to define wave variables for lumped impedances in a manner analogous to distributed systems [Kuo, 1966] [Fettweis, 1986]. We may define f = fin + four, where f,,, is viewed as a force wave traveling into the impedance, and fout is viewed as a force wave traveling out of the impedance. Similarly, we may define v = vin +Vou0. In the distributed case, waves actually travel a distance over time; and there is a wave impedance relation a unit of delay may be factored out of the secondorder transfer function. The result is the efficient digital system shown in Figure 10 defined by, Vout(z) = z-'H(z)Vn(z) + G(z) vo where, ao +z-1 H(z) = 1 + az-1 _ i 1 +2z- 1 + z-Z 1 + aoz-1 (27) (28) (29) with ao defined as, k - ma2 k + ma 2 (3o) Combining Equations (26) and (30), we may arrive at a simple relationship between Rh and a0, ma Rh = - (1+ ao) (31) Sound SynthesisTechniques 416 ICMC Proceedings 1994

Page  417 ï~~Sv0 Figure 11: String Loaded with Mass and Spring 4.2.3 Scattering Junction Connection The next step is to connect the mass and spring system to a string, as shown in Figure 11. Since we have decomposed the lumped mass and spring system into traveling waves, and have picked an associated wave impedance, Rh, we may join the system to the string system at a lossless 3-port scattering junction as shown in Figure 12. Vh Vh+ M V S V+ Figure 12: Attaching the Wave Digital Hammer We define vot0 = vh+ and vi_ = vh-, as was shown in Figure 10, so that traveling wave signals entering the scattering junction are superscripted with a plus sign, and traveling wave signals leaving the junction are superscripted with a minus sign. At the junction point where the mass and spring system connects to the string, the velocity of both sides of the string, vi and v2, and the driving point velocity of the spring, vh, must all equal vj, the junction velocity, V1 = V2 = VhV =-vi(32) v1+ +v1 = V2 + V2- = Vh+ + Vh- (33) In addition, the sum of the forces exerted by the string halves and the mass and spring at that point must be zero, since it is a massless point, fl+f2+fh=fJ=0 (34) fl+4f- +f2++f2-+fh++f =0 (35) Combining these two series junction constraints, (32) and (34), with the wave variable definitions, (2) and (3), and wave impedance relations, (1), we can derive the lossless scattering equations for the interconnection of the wave decomposed mass and spring with the string halves, _ 2 (Rovi+ R-Rov2+ -Rhvh+) (36) Vj-= 2o+ Rh (6 V1- =Vj - Vl (37) v2-= Vj, - v +(38) Vh- = Vj - Vh + (39) These equations say that, as a wave is coming into a junction, 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 other outgoing waves. The relative proportions of this scattering effect is dependent only on the relative wave impedances of the interconnected elements. Figure 13 shows how the mass and spring system may be attached to a 2D digital waveguide mesh [Van Duyne & Smith, 1993]. Junctions marked J are lossless 4-port junctions. The mass and spring system is attached at the 5-port scattering junction marked S. Equations (36) through (39) must be modified appropriately to compute the 5-port junction. Figure 13: Attaching the Wave Digital Hammer to a 2D Digital Waveguide Mesh 4.2.4 Making the Felt Nonlinear Now, to convert our string loaded with a mass and spring to a string being hit with a felt hammer, we must cope with the nonlinearity of the spring "constant". k(Xk) depends on the compression distance, Xk, of the felt. Fortunately we have signals lying around in the mass spring loop which can give us this distance directly. In general, fk = kxk or xk = fk/k for a spring system. If we can find the force being applied to the spring, we may divide it by k to obtain the compression distance. But, the driving point force on the mass and spring system of Figure 9 is equal to the compression force on the spring alone. Hence, Xk = (1/k)(fh- +fh+) = (Rh/k)(v- -Vh+) (40) (41) ICMC Proceedings 1994 417 Sound Synthesis Techniques

Page  418 ï~~Rh/k in Equation (41) may be rewritten in terms of a0 by combining Equations (26) and (30), to obtain 1-a0 Xk= 2a (vj--vh+) (42) in contact with the string. As the hammer felt compresses, Rh (n) ramps up, representing the increased impedance of the compressed felt on the string. 4.2.6 Modeling Hysteresis Loss in the Felt In Figure 14, L computes a small offset of the indexing signal for K which models loss in the felt. It should generally be of the form, With this final piece of the puzzle, we can complete lossless wave digital hammer system as shown in Figure 14. v -(a) r( 'o(") xA(") L Â~ (s) Figure 14: The Wave Digital Hammer 4.2.5 The Complete Digital System Referring to Figure 14, the system loop containing S, H, and one unit of delay, z-1, performs the essential mass and spring modeling. S is the scattering junction defined by Equations (36) through (39). H is the time-varying allpass filter defined by Equation (28). This loop oscillates at an instantaneous frequency of arccos(-ao(n)) and models the second-order mass and spring system. Incidentally, if the computation of S is replaced with, Vh-(n) = --Vh+(n), (43) which is physically equivalent to replacing the string or membrane portion of the system with a rigid wall, then the S -z-'-H loop stands alone as an efficient, one-multiply oscillator. G is the modified integrator defined by Equation (29) which integrates the incoming initial hammer velocity impulse signal taking into account the wave decomposition and the nonlinear effects of the hammer felt. Signal u3(n) is a combination of the velocity waves traveling in the hammer system which represents a normalized force of compression on the hammer felt. X computes the actual felt compression according to Equation (42), taking advantage of the fact that the force of compression on the felt is instantaneously proportional to the felt compression. K is a look-up table for the stiffness coefficient, a0o(n), indexed by felt compression. R computes the effective wave impedance of the hammer system from a0 according to Equation (31). Note that Rh(nl) is zero when the hammer is not touching the string, which means that there is no effect on the rest of the system when the hammer is not us (n) = e(in))(xk(n)- Xk(f - 1)) (44) where xk(n) - xk(n - 1) represents the velocity of felt compression. This may be understood intuitively as a direct implementation of hysteresis by realizing that a resistive loss in the felt will cause the felt to feel stiffer when being compressed and to feel less stiff when decompressing. The springiness of the felt is proportional to the felt compression, whereas the effect of loss is in the direction opposite to the direction of movement. Therefore, Equation (44) offsets the index to the stiffness lookup table, K, such that, when the felt is compressing, the spring "constant" is greater, and when the felt is relaxing, the spring "constant" is smaller. e(xk) should ideally be zero when the hammer is not in contact with the string and should ramp up smoothly from zero as felt compression increases. References [Chaigne & Askenfelt, 1994] Antoine Chaigne and Anders Askenfelt. Numerical simulations of piano strings. I & II. JASA 95 (2) and (3). [Fettweis, 1986] Alfred Fettweis. Wave digital filters: theory and practice. Proc. IEEE 74 (2). [Kuo, 1966] Franklin Kuo. Network Analysis and Synthesis. John Wiley & Sons, New York. [Morse & Ingard, 1968] P. Morse and K. Ingard. Theoretical Acoustics. McGraw-Hill, New York. [Nehari, 1952] Zeev Nehari. Conformal Mapping. McGraw-Hill, New York. [Rodet, 1993] Xavier Rodet. Flexible yet controllable physical models: a nonlinear dynamics approach. Proc. ICMC, Tokyo. [Smith, 1992] Julius 0. Smith III. Physical modeling using digital waveguides. (MJ 16 (4). [Van Duyne & Smith, 1993] Scott A. Van Duyne and Julius 0. Smith III. Physical modeling with the 2-D digital waveguide mesh. Proc. ICMC, Tokyo. Sound Synthesis Techniques 418 ICMC Proceedings 1994