Physical Modeling by Directly Solving Wave PDESkip other details (including permanent urls, DOI, citation information)
This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License. Please contact firstname.lastname@example.org to use this work in a way not covered by the license. :
For more information, read Michigan Publishing's access and usage policy.
Page 00000325 Physical modeling by directly solving wave PDE Marco Palumbi, Lorenzo Seno Centro Ricerche Musicali Via Lamarmora, 18 - 00185 Roma, Italy (email@example.com) Associated to Econa (Interunmiversity Centre on Cognitive Processing in Natural and Artificial Systems) Rome - Italy Abstract We are currently developing a physical modeling method for string instruments based upon the direct solution of the string equation, rather than the simulation of travelling waves - as is instead done in the Waveguide method. Although computationally heavier in comparison with Waveguides, our approach enables any physical parameter to be time-varied during playing, obtaining thus interesting effects from a musical viewpoint Particular attention has been paid in the modeling of the bow, which is the main responsible of the peculiar timbre of bowed strings. The rosin thermal behaviour has been modelled by thresholds. Moreover, the bowing noise is simulated by adding noise to the bow force profile. The whole model runs in real time on today Wintel Workstations for bowed string of the violin or viola class. It produces good sounds not only in sustained, gently bowed articulations, but also in harsh articulations like e.g. "strappato", as well as outside the physically possible parameters zone. The model has been actually used with a primary role by the Italian Composer Michelangelo Lupone in two newly composed electronic and mixed music pieces. Glissandos of string parameters - like fi. the internal dumping - and displacements of the bowing point play a main expressive role in these works. The novelty in the string equation we used is a partial third-order derivative mixed term describing the internal friction of the string, which is responsible of the frequency-dependent dumping in actual strings. The free-motion analytic solutions of-the corresponding PDE is discussed, together to the numerical integration method developed (based on a combination of the sinus transform, and of the well-known Stoermer rule), whose results are analysed from a timbre point of view. We also show a short comparison with Waveguides, and future improvements to continuously varying bowing point and left-hand finger position, and to reduce both asymptotic and actual computational complexity. 1 Introduction Most of people uses Waveguides - delay-lines which simulate the D'Alembert solution of the string PDE in terms of travelling waves  - in order to o synthesize musical instruments by physical modeling. The major role played by this approach is mainly due to the low computational cost of delay lines, and to historical reasons: it is a development of the first method invented - the Karplus-Strong algorithm . Our approach can instead be considered as an evolution of the finite difference method, with the aim to bypass its limitations. Finite difference meth6d is actually a formal implementation of a discrete massspring system, a physical system which inherently shows sets of partials having sub-harmonic ratios. Our string model is instead spatially continuous, thus not suffering from this drawback. Moreover, our personal tendencies and our bonds to contemporary composers pointed our attention to the search of new interesting sounds rather than in the imitation of true, traditionally played, instruments. For this reason we were very interested in developing a synthesis method with parameters time-varying capabilities, and able to simulate physical systems ou, de the "classic" performance parameter domain. Our model (in the form of a computer software for Win tel32 Workstations) was actually used by the Italian composer Michelangelo Lupone to get suggestions about new performing techniques and to compose "Metal string", for string quartet, tape and spatializer (Kronos Quartet, Rome, 1997), and "Canto di Madre" for tape, by commission of the Vatican Radio, Oct. 1998. Some excerpts will be played during the speech. 2 The model 2.1 The string Consider a string of length L, whose free PDE is: Sy(x,t) =c2 -y(x,t) -S -y(x,t) + at2 2 at 1) tSX2 ati at (xlt with the boundary conditions: 2) y(O,t) = y(L,t) =0 Where: c (m/sec) Conservative propagation speed of waves (i.e. for S = S, = 0) S (sec') coefficient of viscous friction with the media. Si (m2/sec) coefficient of viscous internal friction. L (m) length of the string ICMC Proceedings 1999 - 325 -
Page 00000326 A few words about the presence (and the absence) of some terms. The classical dispersive term is lacking: a4 &4y(x, t) This term represents the stiffness of the string, and is responsible of the conservative dispersion. Stiff strings like those of the low section of the piano, have partials in super-harmonic frequency relationship. Strings used in bowed instruments (as in the violin family) have instead very low stiffness, and for this reason we skipped this term. The second, right-hand term, is due to the viscous friction with the air (or the medium). Dissipative terms like this one also introduce disperSion. 2.2 The internal viscous friction term The last, right-hand, third-order mixed term is generally neglected. We couldn't find it in the literature. It represents the effect of the energy losses due to internal, viscous friction, dissipating energy when the curvature changes. This phenomenon is responsible of a main behaviour of actual strings: the higher the frequency partial, the faster the dumping. For instance, if you pluck a string, during the transient you can hear a lot of partials - a rich sound. Whereas in the release part of the sound you can hear just the fundamental Putting S = 0, the pure effect of this term can be investigated. The solution of the resulting PDE is: -t 3) y(x,t)=a-sin(k-x~c- t)-er i.e. travelling forward and backward dumped waves. Superimposing the two components, and taking into account the boundary conditions, one can write: k=n n-i with n=1,2...oo L 1 n 2;2 c 2 L2 4) w = * 2S 4*- -1 2 iL n2 *.r2.Si2 -1 5) Yn (x,t)= 2 e a * e sen(k, x).cos(j|w *t) 2 c L * for n <-, i.e. c, real, and -t 6) y, (x,t) = 2 aer'l -sin(k, -x) - cosh(I), 0) for n >, i.e. o, imaginary. (S.k,2) Modes 5) are oscillatory. Higher modes 6) are overdumped, not oscillatory, as you can see besides in de tail. This is a quite important point: this means that higher modes give contributions only to transients, not to steady-state oscillations. On the other hand, high oscillating modes have large dispersion (4)), but are quickly dumped (7)). Let's take a look to the solutions for modes beyond the critical dumping condition. You can rewrite 5) and 6): y, (x,t)=On(x) -, (t) where: on(x) = sin(k, -x) and -r Vn (t) = cos(In, - t) - eT" for oscillating modes -r yIn (t) = cosh(i,2,. t) - e'" for over-dumped modes. The last one brings to: 1 1 8) W,(=t)=a( (e +e ) Taking into account 6), both exponentials are decreasing. So, the decaying of an over-dumped mode is due to the superimposition of two different time constants. 2.3 The stimulus In our model we consider the coupling of the bow to the string, taking into account the pressure and the speed of the bow, the location of the bowing point xp., and the "stick" / "slip" condition of this point. Our bow is in some way different from that usually implemented by means of continuous speed/force functions. During stick, the bowing point xp. (one of the sampling points) is glued to the bow, and we force the speed of that point to be equal to the bow speed. During the slip state, the bowing force F is added to the forces balance in that point. The bowing force depends on the relative speed, in the way shown in fig. 1 underneath. As you can see, the "stick"-f"slip" transition has a threshold that is twice the "slip"-+"stick" transition (which happens for zero relative speed). The effect of the bowing pressure is to linearly scale the force (both threshold and continuous curve). We don't take explicitly into consideration the temperature of the point - a variable affecting the fluidity and thus the friction behaviour of the rosin (a natural polymer with a complex behaviour). The shown threshold emulates however the cooling effect of the rosin during the stick time. The noise due to the rubbing of the bow over the string is modelled by means of white noise added to the F curve, in order to preserve the curve as a "maximum value" of the random value. This is a true "noise model", and adds Darametric noise to the model, which imparts to the string motion a "micro-chaotic" behaviour due to the random interaction between bow and string. 3 The sound The sound is taken as the displacement of the spatial sampling point nearest to the bridge. One can think of this point as a "secant" approximation of the tangent to the string in the bridge position, Le. to the strain against the bridge. -326 - ICMC Proceedings 1999
Page 00000327 F Sv(t)-by(x,t) fig. 1 - State diagram of the bowing point Force versus point-bow relative speed. 4 The method of calculus 4.1 Space The idea consists in the use of a sinus transform of the string shape, truncated to the first N components, which gives us an analytical expression of the shape, with N degrees of freedom. The N coefficients of the sinus series can be thus derived from just N points on the string. Sinus transform also implicitly satisfies the boundary conditions 2). Truncation to N components, band-limits the vibration of the string, but the considerations about the particular contribution of higher partials developed in 2.2 tell us that also actual strings are, in some way, "band-limited". Let's write the sinus series, being, at a given instant, x,,x2... Xy the abscissas and y1,y2... y the corresponding vector of ordinates of N points: N y, = a, -sen(n x- x) nsl The series can be in a more compact way rewritten: y = Ma where: y = [y, a = [aO] and M = iw,,n J being M,, = sen(n, r * x,) Thus, at any given instant: a = M-'y The matrix M is inherently invertible: the orthonormality of the sinus functions gives us (in principle) a complete freedom of choice of the points x,. They can be unevenly spaced, provided that there are no coincident couples. In practice, they must not be too much near each other, in order to avoid precision losses in numerical calculations. The calculation of the curvature of the shape from the coefficients a is a straightforward task: --y(x) = a.n2 2 -sin(n.- r x) ax a2 So, being C = [C,, with C, = y(x), the vector of the curvatures in the "sampling points x: C=M'a where MA, =n2 *r2. sin(n. r xi) Please note, by the way, that M' depends on the choice of the Npoints abscissas. So, the task of computing the curvatures starting from the knowledge of the coordinates of the "sampling points" can be summarised in this way: C = MM-xy, ie. C= My where M, = M'M- is the "curvature matrix". Several optimizations can be made in the vectormatrix product. If the points are evenly spaced, M, shows symmetries that can be exploited to reduce the number of mul-adds. We obtained in this way a reduction by a factor of two, for a 16 x 16 matrix. Another way is the use of FFTs to calculate the direct and inverse sinus transform. For small N's like 16 we could not find any improvement using an FFT, but further researches can be made for ad-hoc simplifying schemes. 4.2 Time Once calculated as shown in 4.1, the curvatures can be used to decouple the N points, getting from 1) N separated motion equations of N dumped oscillators of which you know the recalling force (the curvature). The special nature of the internal dumping factor requires the time derivative of the curvature in each ICMC Proceedings 1999 - 327 -
Page 00000328 point, a task easily performed calculating a finite difference using a previously stored value. To integrate the motion equation, we used the Stoermer rule, known and used for conservative systems since 1913. This second-order method of integration brings to a small error whose effect is to overestimate the frequencies of the partials, bringing to a sort of "numerical stiffness". To reduce this effect we used over-sampling. A dumped string is not conservative, of course, so in the attempt to reduce the over-sampling factor we tried the use of higher order, prediction-correction methods. The following conclusion are to be considered yet provisional, but the results showed - as expected - better frequency precision, introducing however an error adding further, purely-numeric dumping. This makes long, sustained sounds impossible. Stoermer rule, on the contrary, being conceived for conservative systems, gives perfectly sustained everlasting sounds, the dumping acting as a small correction to the steadystate conservative motion. As a temporary conclusion, it seems that finer methods of integration are good for more damped systems than for musical applications. S Results The model has been widely tested and used with N=16 and a four times over-sampling factor. In these conditions, it easily runs in real time on a Pentium II 350 system @ 44.1 KHz sampling frequency. Further researches are in progress in order to make feasible a continuous placement of the bowing point by means of interpolation via polyphase filtering, requiring thus memory but no further computational resources. In a similar way can also be implemented supplementary bonds along the string f.i. slight fingerings like those used to produce harmonics. Simple models of pinch and percussion brought very satisfactory results from the standpoint of imitation - thanks mainly, in our opinion, to the correct modeling of the frequency-dependent dumping due to the internal viscous friction losses. Bowed emissions, with the specific goal to imitate actual instruments gave good results for little variations of the pitch, particularly as to the bow attack transients and the evolution of Ionglasting sounds. The auditory results remain good even in the case of varying pressure, velocity and position in non-standard conditions of playing. 6 Waveguide comparison The main difference between the model shown here and Waveguides is the parameters time-variance capability. Waveguides mimic sampled travelling waves, obtaining computation reductions via the commutativity of the string segments transfer functions, accumulating them into just one (or few) global transfer func tion in one (or few) point In this way, linearity and time constancy of the discrete subsystems are necessary conditions. No time-variance string parameters nor glissando's are allowed with Waveguides. Moreover, the effect of the internal friction is implemented by means of a commuted low pass filter. Easy-toimplement,.IIR filters does not have a good auditory effect, because of their phase shift whose effect is an audible dispersion of partials, quite different from the one described in 4). In our opinion in the bowed strings the most important role is played by the bow, rather than the string. Applying our bow to a Waveguide will likely bring to good results, provided that you are not interested in time-variance. The computational complexity of Waveguide models still remains better that the direct method here shown, whose complexity is O(IN) for direct matrix methods, or O(N. Id(N)) for FFT-based sinus transforms, even pushing Waveguide over and over to implement finer auditory features. On the other hand, computational complexity can be less and less important as the power of GSS and DSPs increases, and if you consider the time-variance an important expressive, musical requisite. 7 References  Raman, C.V. (Ramaseshan S. Editor), "Scientific Papers ofC.V. Raman: Acoustics" 1989. MlTPress. [21 Cremer L., Allen J.S., "The Physics of the Violin" 1985. AT Press. [31 Hiller H., Ruiz P., 1971. "Synthesizing Musical Sounds by Solving the Wave Equation for Vibrating Objects: Part 1" - "Part 2" JAES 12 July (6), Agust (7)  Schelleng J.C. "The Bowed String and the Player" 1973. JASA 53 (1). 5] McIntyre ME. et al., "Fundamentals of Bowed String Dynamics", 1979. Acustica 43.  Schumacher R.T., "Self-Sustained Oscillations of the Bowed String." 1979. Acustica 43.  McIntyre ME. et at, "Aperiodicity in Bowed String Motion", 1981. Acustica 49.  McIntyre M.E., Schumacher R.T., Woodhouse J. "On the Oscillations of Musical Instruments", 1983. JASA [91 Karpus, K., Strong, A. 1983. "Digital Synthesis of Plucked String and Drum Timbres." C/ 7,2 pp.43-55 [101 Jaffe, D.A, Smith, J.O. I "Extension of the KarplusStrong Plucked-String Algorithm." Summer 1983, CJ 7,2 [I1] Chafe, C. "Simulating Performance on a Bowed Instrument" 1989 Current Directions in Computer Music Research - Edited by Matew, M. and Pierce, J. MIT Press.  Woodhouse, J., 1992. "Physical Modeling of Bowed Strings." CM2J, 16,4, pp. 43-56 [t3] Smith, J.O. i 1996. "Physical Modeling Synthesis Update." Ct/. 20,2, Summer 1996 pp.44-56 (141 Smith, J.O. m, 1996. "Discrete Time Modeling of Acoustic Systems with Applications to Sound Synthesis of Musical Instruments." Proceedings of the Nordic Acoustical Meeting, Helsinki, June 12-14 [(151 Valimaki V., Huopaniemi J., Karjalinen M., Janosy Z. "Physical Modeling of Plucked String Instruments with Application to Real-Time Sound Synthesis" 1996. JAES 44 (5)  Mastropietro A., "Ritratto di musicista da sperimentatore: Michelangelo Lupone", Sonus n. 18, Sonus Edizioni Musicali, Potenza 1998.  Lupone M., "Spazio, arco e metallo" in Musica tecnologia interzioni, Quaderi della Civica Scuola di Musica, anno XIV N. 26 - Comune di Milano, 1999. - 328 - ICMC Proceedings 1999