Page  407 ï~~A Simplified Approach to Modeling Dispersion Caused by Stiffness in Strings and Plates Scott A. Van Duyne CCRMA, Stanford University savd(ccrma. stanford.edu Julius 0. Smith, III CCRMA, Stanford University j os(ccrma.stanford. edu Abstract Adding allpass filtering to a delay line loop can introduce inharmonicity into the resonant loop modes. We propose a simplified filter design method to match loop modes to measured stiff string partial frequencies with a set of cascaded multiply-free one-pole allpass filters. This digital structure corresponds physically to terminating a distributed string model with a non-uniform lumped mass and spring transmission line chain. Promising results are demonstrated in re-tuning modes on the 2D digital waveguide mesh as well. 1 Current Methods The use of allpass filtering to introduce dispersion into feedback delay loops has been known for some time [Jaffe & Smith, 1983]. Recently, [Paladin & Rocchesso, 1992] have been able to generate excellent stiff string sounds on the MARS workstation in real-time using a group delay based analysis approach to the stiff string equation (1), and using one of the allpass filter design methods in the literature such as that proposed by [Yegnanarayana, 1982] to find allpass filter coefficients for a desired group delay response. However, the order of allpass filters found by these methods can become large, requiring a fair amount of computation; and finding the coefficients requires a sophisticated filter design algorithm. We propose here a simplified approach to the filter design, building up allpass filters out of one-pole blocks. These filters are designed to be multiplyfree when implemented directly in hardware by virtue of choosing a coefficient which may be computed with a binary shift and at most one add. 2 The Stiff String The equation of the stiff string may be written as [Morse & Ingard, 1968], this equation, only appreciate how stiffness affects the tuning of the string modes. In the ideal flexible string case, the modes are harmonically related. We can say, p(n) =nfo (2) where fo is the fundamental frequency, n is the harmonic number, and p(n) is the frequency of the nth harmonic. For the stiff string, however, higher spatial frequencies travel faster than lower spatial frequencies, causing the higher partials to become farther apart than the lower. The partial frequencies may be approximated as, p(n) n (a +,On2) (3) [Morse & Ingard, 1968], where /3 is a small inharmonicity factor, and the fundamental frequency is p(1) = a + /3. As n gets large, and higher order terms come into play, this approximation breaks down. From Equation (3), we can say that the frequency separation between the partials, S(po), at some point in the spectrum, Po, is approximately, S(po) _ p(no +.5) - p(no -.5) ti p'(no) (4) where po = p(no). The rate at which the partial frequency separation changes per unit of spectrum height, i. e., the rate of stretching at a given spectral frequency, Po, may be defined as the derivative of S with respect to Po, S'(p S == -_p"(no) S'(O- dpo u d d pfl '(no) (5) T-ffj - QS,Â~-4y= PS-2 y (9X2 Ox4 Ox2 (1) This equation differs from the ideal flexible string equation in that a fourth-order term is introduced causing the string to behave more like a bar as it gets stiffer. For our purposes, we need not solve Figure 1 shows a plot of the rate of stretching, S', for a hypothetical low AO piano string with a reasonable inharmonicity factor of /3 =.001 assumed, and partial frequencies defined as, pAo(n) = n (27.499 +.001n2) (6) 1CMC Proceedings 1994 407 Sound Synthesis Techniques

Page  408 ï~~If the value of S' were constant at.01, the separation between adjacent partials would increase at the rate of 10Hz per kHz of spectrum. At the fundamental, the partials would be about 27.5Hz apart; at 1000Hz they would be about 37.5Hz apart, and so on. Actually, the plot shows that the rate of stretching is not flat. Rather, it ramps up from zero (rate of zero meaning the partials are exactly harmonic), reaches a maximum, and drops slowly as we get higher into the spectrum. It is this qualitative shape of the rate of stretching curve which leads us to model dispersion in stiff strings with cascaded one-pole allpass filters. OaT -o to1).01.005 so (po) Figure 2: Phase Response of One-Pole Allpass unwrap the values of the usual computer generated arc tangent function by subtracting ir whenever the arctan function is positive. Figure 2 shows overlaid plots of q(wT) for values of aG ranging from -0.8 to 0.8. 4 The Dispersive Delay Loop Figure 3: Dispersive Delay Loop System A resonant loop consisting of a length N delay line and a block of M cascaded one-pole allpass filters is shown in Figure 3. The transfer function one time around the loop is, i 10kHz 20kHz Figure 1: Rate of Stretching for Hypothetical AO 3 The One-Pole Allpass The one-pole digital allpass filter may be defined by its difference equation, y(n) = aox(n) + x(n - 1) - aoy(n - 1) (7) where x(n) is its input signal and y(n) is its output signal, or by its Z-transform, a ao+z-' = +ao1 1 (8) We take the coefficient, ao, to be a real number between -1 and 1. The frequency response of this filter is obtained by replacing z with ejiwT, where w is the radian frequency and T 1/f, is the sampling interval. The magnitude response, IH(eJwT)I, is always 1; therefore, we may re-write the frequency response in the following polar form, H(edwT) = ejO((wT) (9) where the phase response is 45(wT) =_ LH(e'wT). We may derive an analytic expression for the phase response, q!(wT) = arctan Im{H(eiwT)}(10) Re{ H(eiwoT)} (0 = arctan2(ao2 -1)wsin tTT1 2ao + (ao 2 + 1) cos woT (11) where the arc tangent function has been unwrapped. Since we restrict w to positive frequencies, we can L(z) zNH(z)M (12) We compute the phase response of L to be, 6(wT) LL(eJwT) = -NwT + MO(wT) (13) For frequency WkT to be supported on this loop, O(WkT) must equal -k27r, for some positive integer k. In fact, k = 1 implies wkT is the fundamental, k = 2 implies wkT is the second partial, and so forth. Figure 4 compares the supported modes for a loop of the form z-( N+ M) with a loop of the form z-NH(z)M, where N = 8, M = 1, and ao = -0.8. The wiT are the harmonics of the pure delay loop, and the wiT are the inharmonic partials of the delay plus allpass loop. Comparing the wiT with the wiT, it seems clear that C.iT is less than w1T, that the spacing between w; T and &2T is about the same as that between 1T 1T S2TOT o)T jI T I 44 )4T -2x. F(geT) e n Figure 4: Modes on z-NH vs. Modes on z-(N+l) Sound SynthesisTechniques 408 ICMC Proceedings 1994

Page  409 ï~~wiT and w2T, and that the spacing between w2T and w3 T is greater than that between w2T and w3T. This means that the wiT are stretching partials. In general, the choice of a0 negative will produce a stretching in the resonant modes of the loop. The more negative it is, the more the stretching will be, at least in the lower frequency range. We can make a qualitative comparison between the kind of stretching in the delay line loop with that found in stiff strings by computing a rate of stretching function for the digital loop. We derive So(wkT) = -21r/O'(wkT) to be the separation between partials on the loop at a particular point in the normalized frequency spectrum, wkT, where wk is the kth partial frequency in radians per second. The rate of partial stretching per unit of spectrum on the digital loop is the derivative of So with respect to wkT, = '(wkT) (14) So'(kT)= -2r O(wkT )Z Figure 5 shows a plot of Se' with N = 8, M = 1, and ao = -0.25. Comparing the physical rate of stretching, S', in Figure 1 with the digital loop rate of stretching, So', one can see a qualitative correspondence: each plot ramps up from zero, reaches a maximum, then falls off. If we choose N, M, and a0 well, we can approximate partials of a real stiff string..04.02 Vd2 n (OkT Figure 5: Stretching of Partials on a Delay Loop 5 Fitting Loop Parameters to Measured Data To design a resonant loop of the form z-N HM, we must find a good a0, M, and N. We reduce the field of search by requlring that the filter coefficient, a0, be a power of 2, or 1 minus a power of 2, so that the multiply may be implemented efficiently in hardware as a binary shift plus at most one addition. Coefficients closer to zero provide the most even stretching of partials, but the number of ailpasses needed will be greater; whereas, coefficients closer to -1 will necessitate less filters in the loop, but the effect will be biased in favor of stretching at the lower frequencies. We can eliminate N as a variable by requiring that one of the partial frequencies on the loop be held constant. The necessary value for N may be computed as, 27rno + Mc(27rp(no)T) 2wrp(no )T (15) Choosing no = 1 holds the fundamental fixed. In the case of missing data, or in the case of very low piano strings, where the perceived pitch is more a function of the partials an octave or two above the fundamental, it may be desirable to hold some other partial fixed. Assume that we have a set of partial frequency measurements, pi = p(ni), for a set of partial indices, ni. Given ao, we can compute the apparent partial number, no(ni), on the delay line loop for each of the measured partial frequencies, p(ni), no(ni) = O(2irp(ni)T)/(-2ir) (16) and compare it to the actual partial number for that frequency, ni. If the digital loop is a perfect match for the data, then no(ni) will equal ni. If not, we can nonetheless minimize a least sum of squares error function, e(M) = > (no(ni)- hi)2 (17) As an example, we measured the first 22 partial frequencies of a recorded E2 piano tone. Choosing ao = -.75, we obtained the best fit parameters, M = 2 and N = 255. Figure 6 plots the difference between the apparent partial index on the digital loop, no(n:), and the actual partial index for each of the measured partial frequencies, n, for values of M from 0 (no allpass filters) to 5 (five allpass filters). It is fairly clear from the figure that M = 2 is a good choice. That is, we need only two allpass filters, with coefficient a0 = -.75. If we choose a0 = -.5, similar analysis leads to a choice of N = 221 and M = 16. In this case, the 16 allpass filters can be implemented in parallel bucketbrigade hardware if z221'H'6 is re-formulated as z-205 (z-' H) 16, double buffering the allpass filters. no(n) - n 1 0.75 0.5 0.25 -0.25 -0.5 -0.75 -1 i t 7 i f 1 1 I 7 I A i = t i i, M=O M=1 M=2 M=3 M=4 M=5 n--- Figure 6: Approximating Measured Data ICMC Proceedings 1994 409 Sound SynthesisTechniques

Page  410 ï~~6 Physical Interpretation 7 Dispersion in Plates,l M=2 7 M=a M=I1 Figure 7: Mass/Spring Chain String Terminations A delay line loop with a one-pole allpass filter in it may be interpreted as an ideal string terminated with an ideal spring, as is shown in Figure 7 for the M = 1 case [Van Duyne, Pierce & Smith, 1994]. The loop computes traveling force waves, where the allpass filter is a transfer function from the right-going force wave, F+, to the left-going force wave, F-, at the spring termination. Given an arbitrary termination impedance, F(s)/V(s) --R(s), the force wave transfer function is, F- _R(s) - Ro F =R(s) +R1o (18) where R0 is the wave impedance of the string, which is the geometric mean of its tension and mass density. When two transfer functions corresponding to two termination impedances, R1 (s) and R2(s), are cascaded in the loop, we find, F- (R1 - R (R2 - Ro\ _ -Ro 1 F =k, = +R, +RO R2+Ro =R+R0 (19) where, F =h=R1R2 + Rol V n +(20) Hence, comparing Equations (18) and (20), we see that the physical termination impedance actually modeled by the two cascaded wave transfer functions is R(s). Using a similar analysis procedure, we may derive the physical string termination impedances for the cascaded digital allpass filters, HM, described earlier. A detailed treatment of allpass modeling of dispersion in the two-dimensional case is not possible here. However, much as in stiff strings, the higher frequencies in stiff membranes travel faster than lower frequencies, causing a "stretching" of the modes. Figure 8 illustrates the spectral effect of terminating a 10 x 10 2D digital waveguide mesh [Van Duyne & Smith, 1993] with single allpass filters around the rim. The six plots in the figure have been scaled horizontally to align the de-tuned fundamental mode in each case so that the mode stretching effect for various choices of allpass coefficient, a0, is clarified. Since the 2D mesh is designed for parallel multiply-free hardware implementation, our simplified multiply-free allpass approach is ideally suited to it. 0 v/4 x/2 3,r/4 u I............... ao = -125..... ao _.2: I.hi i tttIIiif 1I IIt. a LIIi~i.., ii --......_____._..___ttll~ilbldiiri....t:~ ~....:.:.......11 f-11hh, ii::1.ulrbuW,,, h ".-- 0 /4x/2 3/4 g Figure 8: Stretching Modes on the Mesh References [Jaffe & Smith, 1983] D. Jaffe and J. Smith. Extension of the Karplus-Strong plucked string algorithm. CMJ 7 (2), 43-45. [Paladin & Rocchesso, 1992] A. Paladin and D. Rocchesso. A dispersive resonator in real time on MARS Workstation. Proc. ICMC, San Jose. [Morse & Ingard, 1968] P. Morse and K. Ingard. Theoretical Acoustics. McGraw-Hill, New York. [Van Duyne, Pierce & Smith, 1994] S. Van Duyne, J. Pierce, and J. Smith. Traveling wave implementation of a lossless mode-coupling filter and the wave digital hammer. Elsewhere in these proceedings. [Van Duyne & Smith, 1993] 5. Van Duyne and J. Smith. Physical modeling with the 2-D digital waveguide mesh. Proc. ICMC, Tokyo. [Yegnanarayana, 1982] B. Yegnanarayana. Design of recursive group-delay filters by autoregressive modeling. IEEE Trans. of A coust. Sp. and Sig. Proc. ASSP-30, 632-637. H' = ki/s H2 ms+k2/s H3 (k3/s)fJ(m2s + k4/s) H4 = m3s + ((k,/s)11(m4s + k6/s)) (21) (22) (23) (24) where allb ab/(a + b) is a parallel connection of two elements, and all the spring stiffness constants, k=, and masses, mi,, are dependent on kl and the string wave impedance, R0. The physical meanings of these termination impedances are illustrated in Figure 7. We conclude that a delay loop with cascaded one-pole allpass filters corresponds to a distributed transmission line string model extended by a non-uniform lumped transmission line model of chained masses and springs. Sound SynthesisTechniques 410 ICMC Proceedings 1994