# On the Use of Weighted Sample Methods in Digitizing the Clarinet Equations

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 00000001 ON THE USE OF WEIGHTED SAMPLE METHODS IN DIGITIZING THE CLARINET EQUATIONS Federico Avanzini, DEI-CSC University of Padova. Via Gradenigo, 6/A - 35131 Padova, Italy. Email: avanzini@dei.unipd.it Abstract- A general scheme for time domain simulations of single reed wind instruments is proposed, which makes use of accurate, stable and efficient numerical techniques. First, the K method is used in order to solve non-computability problems in the discrete model, then the differential equations are discretized using the 1-step Weighted Sample method. Frequency and time analysis shows that the physical properties of the continuos model are mapped with good accuracy in the discrete domain. Time domain simulations show that the model reacts consistently to changes in significant physical parameters such as the reed resonance and damping. Such quantities are related to player's embouchure, and can therefore be used together with the blowing pressure as musically meaningful control parameters in real-time simulations. Keywords- Physical Models, Weighted Sample Methods, Non-linear Systems, Wind Instruments. I. INTRODUCTION Once a set of differential equations has been obtained which accurately describes the physical functioning of a musical instrument, the continuos time model has to be discretized in order to develop numerical simulations. The problem of finding a suitable discretization technique for the continuous system is often neglected in literature, and the equations are usually discretized using simple numerical methods such as the Euler method, the bilinear transformation or the impulse response invariance. However, each of these methods exhibits some drawbacks: frequency warping, artificial damping, aliasing, as well as poor accuracy in time. In this work we present a discrete time model for the single reed, which makes use of accurate numerical techniques. First, the K method [Borin et al., 2000] is used in order to solve a non-computable loop in the reed equations, then the linear differential equation of the system is discretized using the 1-step Weighted Sample (WS) method [Wan and Schneider, 1997]. Analysis both in frequency and in time domains shows that the resulting model well reproduces the physical properties of the continuous system. Sec. II briefly describes the model and the numerical techniques; in Sec. III the properties of the "digital reed" are analyzed; Sec. IV illustrates some significant results obtained from simulations. II. DISCRETE TIME MODEL Table I summarizes the meaning and values for the variables and constants that will be used throughout the paper (except for Sec. IV-C); they are the ones used by Schumacher [Schumacher, 1981]. TABLE I VARIABLES AND CONSTANTS IN THE REED MODEL. Quantity Sampling period Sampling rate Reed equil. position Reed displacement Reed relative displ. Reed mass/area Effective reed area Reed resonance freq. Reed damping Amplitude parameter Blowing pressure Mouthpiece pressure Pressure drop Mouthpiece flow Flow through the slit Wave impedance Speed of sound Length of the resonator Pr. wave from the bore Pr. wave to the bore Symbol Value Ts Fs H Xr Yr Pr Sr Wr gr A Po P Ap Ap u Uf Zo c L p P+ [s] [s-1] 0.4 10-3 [m] [m] [m] 0.0231 [Kg/m2] 1.46. 10-4 [m2] 23250 [rad/s] 3000 [s-1] 0.0797 [m3/(N2/3s)] [N/m2] [N/m2] [N/m2] [m3/s] [m3/s] 2290133 [Kg/m4s] 355 [m/s] 0.3 [m] [N/m2] [N/m2] A. Equations of the model The bore is described by a basic and efficient waveguide model [Smith, 1998], which simulates the propagation by means of two delay lines of length mL = F8L/c. Losses in the wave propagation, as well as wave reflection at the bell, are taken into account by inserting a digital low-pass filter Rd(z) in the delay loop. The cutoff frequency is chosen to be fc = 1.5 kHz. Therefore, the digital bore takes the incoming pressure wave p+ from the exciter, and returns to it an outgoing pressure wave p- given by P-(z) = -Rd(z) -2mL P+(z). The reed model is the one proposed by Schumacher [Schumacher, 1981], and is described by the set of

Page 00000002 equations: tr [yr(t) + grr(t) + Wyr(t) = -wp(t), Uf(t) = u(t) + S-rU(t), Ap = A-3/2sgn(uf) Uf 3/2/X, III. PROPERTIES OF THE DIGITAL REED Consider the transfer function H,(s) of the relative displacement yr versus pressure drop Ap; this is a simple harmonic oscillator: where Ap = po - p and /r = Xr - H. The reed is supposed to be a rigid body, so that it can be modeled as a second order harmonic oscillator. The flow uf through the slit is the sum of the flow u at the mouthpiece and the flow Srýv produced by the reed motion. The third equation describes the relation between the pressure drop Ap and the flow uf, and is derived by Bernoulli's law for perfect fluids. B. Numerical techniques When the linear differential equation in system (1) is discretized, a delay-free non-computable loop is produced, which connects the variables Xr and Ap; due to the presence of a non-linearity, solving this loop is not trivial. Some authors (see [Gazengel et al., 1995]) restore computability by inserting a fictitious delay in the loop; however, this introduces an error in the numerical scheme, and can lead to instability. A general solution to this problem is given by the K method, recently proposed in [Borin et al., 2000]. The method is independent on the chosen discretization technique; it operates a linear variable transformation in such a way that a new system can be defined, in which the noncomputable loop is solved without introducing errors in the computation (see also [Avanzini, 1999]). Concerning the discretization technique, usual choices in literature are the Euler method, the bilinear transformation or the impulse response invariance. In this work we use the 1-step Weighted Sample (WS) method. The method, introduced by Wan and Schneider [Wan and Schneider, 1997], discretizes a linear system ~i = Aw + bu by performing a linear interpolation of the input u. Details about the method are given in [Wan and Schneider, 1997], [Avanzini, 1999], here we stress some of its noticeable qualities: * Stability is preserved for any choice of the sampling rate F,. * The discrete time equations produced by the method are linear equations with constant coefficients, which can be implemented with low computational costs. * The order of is preserved, i.e. discretizing a N-order differential equation produces a N-order difference equation. Previous works [Avanzini, 1999], [Avanzini and Rocchesso, 2000] have been devoted to comparison of this method with other techniques (namely, the bilinear transformation, Adams-Moulton methods, higher order WS methods), and analysis both in time and frequency domains have shown that the 1-step WS method seems to be the technique that best preserves the properties of the continuous system. 1 1 II(S) = Pr s2 + grs + w2 Ar 8 r (2) The discrete transfer function Hwsi (z), obtained by applying to Eq. (2) the 1-step WS method, is given by Hwsi(z) = [1, 0]. [zI - (T)]-1 W [1, z - 1]T, (3) where b(t) = exp 0-g -2 1 t) and calculation of W involves calculation of the integrals 4fT 1(T _-t)0t'dt (see [Wan and Schneider, 1997], [Avanzini, 1999] for details). Again we stress the fact that the computation of these matrices is done offline. Following some ideas developed in [Gazengel et al., 1995], we focus on three physically meaningful parameters: the resonance frequency Wr, the oscillator stiffness (defined as |Hc(0)I-1 = urW), and the damping coefficient g, = WI - w2 (w1,2 are the 3 dB cut-off frequencies). From (3) one can see that the two poles of Hwsi (z) are pi = exp(AiTs) (i = 1, 2), Ai being the poles of He(s). Therefore, if Im(Ai) < 0, then |pi| < 1, so that stability is preserved for any choice of T,. Moreover, Im(Ai) is mapped into F, -arg(pi), so that the oscillator resonance is preserved. Fig. 1 shows the behavior of pi for different sampling rates (p2 is its complex conjugate). From this it can be seen that the system is stable, as predicted by the theory; moreover, a numerical damping can be observed, which is more sensible at low F,. 0. 9 F =20 kHz 0.8 0.7 -:0.6 - S0-5 S0.4 - 0.3 - 0.2 0.1 - F =100 kHz 01 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 Re[p (F)] 0.4 0.6 0.8 1 Fig. 1. Behavior of the pi's for different sampling rates (F, = 20...100 kHz). Fig. 2 shows the frequency response Hwsi (ejw), derived from Eq. (3) with F, = 20 kHz; this is compared with the continuous response, HIc(jw). Again, it can be seen that the resonance is preserved, and a small

Page 00000003 x 10 B' -4 a) O 2. C -a ca 2. IV. RESULTS FROM SIMULATIONS A. Comparison with the quasi-static model A comparison between our dynamic reed model and the well known quasi-static approximation [Hirschberg et al., 1996] can be obtained by plotting the corresponding uf vs. Ap phase diagrams for the steady state signals. In Fig. 4 we show an example of such phase diagrams; the dynamic model exhibits an interesting behavior: uf and Ap move along a hysteretical (a) path. This phenomenon comes from to the presence of memory in the dynamic reed model, while the quasistatic approximation results in a memoryless model; we can say that in the dynamic model the attractor is not a curve in the plane anymore: instead, it is a closed orbit in a higher dimensional phase space. 3000 4000 f [Hz] 7000 (b) E x 10 4.5 B L 0 1000 2000 3000 4000 5000 60 Ap [Pa] Fig. 2. Frequency responses (F8 - id line) and |Hwsi| (dotted); (b) arg(Hwsi) (dotted). 20 kHz); (a) He (solarg(Hc) (solid line) and 00 numerical damping is introduced; stiffness is well preserved. The same results are obtained by comparing the impulse response hwsi with the continuous one, hc. Fig 3 shows an example of such responses, with F, = 20 kHz. x 10 - EL 0.5 o -o E - -1.5 0 0.5 1 1.5 2 2.5 t [sl x1o-3 Fig. 3. Impulse responses (Fs = 20 kHz); he (solid line) and hwsi (dotted). Summarizing, the digital reed preserves the physical properties of the continuos system with good accuracy, even at low sampling rates (F8 = 20 kHz). Less accurate methods do not give the same results; as an example (see [Gazengel et al., 1995]), the Euler method produces a discrete frequency response which does not fit at all with the continuous one. Fig. 4. Phase diagram for Uf vs. Ap (F8 = 44.1 kHz). Quasistatic model (solid line) and dynamic model (dotted line). B. Threshold pressure A second study on the simulations concerns the threshold pressure pt, i.e. the value for blowing pressure above which oscillations can take place. A rough estimate for pt is given in the quasi-static approximation, from which the expression pt - H /PrW /3 is derived [Hirschberg et al., 1996]; with the values listed in Table I we find pt - 1664 Pa. However, as observed by Keefe [Keefe, 1990], this value underestimates the true Pt. This estimate has been compared with experimental results from the simulations. Results are summarized in Fig. 5; from this we can see that, for increasing F8, Pt converges asymptotically towards a "true" value which is found to be close to 1802 Pa. Moreover, the system exhibits a good robustness with respect to F8. C. Transitions to upper registers As many authors have pointed out (see for instance [Thompson, 1979]), both the resonance wr and the damping g, play a key role in helping transition to the second register (clarion register). Experiments with artificial lips on real clarinets have proved that, if the reed resonance matches a low harmonic of the playing

Page 00000004 -1780 Cý (a) 2 3 4 5 6 7 8 9 10 Fs [Hz] xlo4 Fig. 5. Dependence of the threshold pressure pt on Fs. frequency and the damping is small enough, the clarion register can be produced. Moreover, an exceedingly low damping causes the "reed regime" (squeaks) to be produced, i.e. the oscillation is governed by the reed resonance. From a musical standpoint, squeaks are often explained as a consequence of insufficient breathing, while the fundamental register comes in as the blowing pressure is increased. All these effects are well reproduced by the model. Fig. 6a shows a transition to the clarion register, obtained by matching w, to the seventh harmonic of the playing frequency and lowering gr down to 980 s-1 In Fig. 6b a transition to the reed regime is showed. Squeaks are also observed in simulations by injecting a blowing pressures which is just below the threshold. V. CONCLUSIONS We have shown that applying the 1-step Weighted Sample method to the single reed equations leads to a discrete time system which accurately preserves the physical properties of the continuous one. Time simulations obtained from the model show that it reacts consistently to changes in physical parameters such as wr and g., which are related to player's embouchure. Moreover, the model results in low computational costs, so that it can be used for accurate real-time simulations of clarinet-like instruments, in which both the blowing pressure po and embouchure-related quantities are taken as musically meaningful control parameters. We also stress the fact that this digital reed can be coupled to more accurate bore models, even with different geometries, thus leading to a general tool for synthesizing single reed wind instruments. Finally, we notice that the numerical techniques developed in this work can be used in simulations of higher dimensional, non-linear fluid mechanical systems, such as the double reed and the lip reed, or multiplemass models [Ishizaka and Flanagan, 1972] of the vocal cords. -500 - -1000 _W-1500 -2500 -3000 -3500 - 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 t [s] (b) Fig. 6. Upper registers (F, = 44.1 kHz); (a) clarion register (Wr = 12580 rad/s, gr = 980 s-1, po = 1800 Pa); (b) reed regime (wr = 1.9770 rad/s, gr = 300 s-1, po = 1900 Pa). REFERENCES [Borin et al., 2000] Borin, G., De Poli, G. and Rocchesso, D. 2000. "Elimination of Delay-free Loops in Discrete-Time Models of Nonlinear Acoustic Systems." IEEE Trans. on Speech and Audio Process., in press. [Wan and Schneider, 1997] Wan, C. and Schneider, A. M. 1997. "Further Improvements in Digitizing Continuous-Time Filters." IEEE Trans. Sig. Process., 45(3):533-542. [Schumacher, 1981] Schumacher, R. T. 1981. "Ab Initio Calculations of the Oscillations of a Clarinet." Acustica, 48(2):71-85. [Smith, 1998] Smith, J. O. 1998. "Principles of Digital Waveguide Models of Musical Instruments." Kahrs, M. and Brandenburg, K., eds., Applications of DSP to Audio and Acoustics, pp. 417-466. Kluwer Academic Publishers. [Avanzini, 1999] Avanzini, F. 1999. "Stable and Accurate Numerical Methods in Physical Models of Musical Instruments. A Study on the Clarinet." Proceedings of the Diderot Forum on Mathematics and Music, pp. 11-19, Vienna. [Avanzini and Rocchesso, 2000] Avanzini, F. and Rocchesso, D. 2000. "Efficiency, Accuracy, and Stability Issues in Discrete Time Simulations of Single Reed Wind Instruments." Submitted. [Gazengel et al., 1995] Gazengel, B., Gilbert, J. and Amir, N. 1995. "Time Domain Simulation of Single Reed Wind Instruments. From the Measured Impedance to the Synthesis Signal. Where are the Traps?" Acta Acustica, 3:445-472. [Hirschberg et al., 1996] Hirschberg, A., Pelorson, X. and Gilbert, J. 1996. "Aeroacoustics of Musical Instruments." Meccanica, 31:131-141. [Keefe, 1990] Keefe, D. H. 1990. "On Sound Production in Reed-Driven Wind Instruments." Technical report, (No. 9003), Systematic Musicology Program. School of Music, Univ. of Washington, Seattle. [Thompson, 1979] Thompson, S. C. 1979. "The Effect of the Reed Resonance on Woodwind Tone Production." J. Acoust. Soc. Am., 66(5):1299-1307. [Ishizaka and Flanagan, 1972] Ishizaka, K. and Flanagan, J. L. 1972. "Synthesis of Voiced Sounds from a Two-mass Model of the Vocal Cords." Bell Syst Tech J., 51:1233-1268.