Digital Waveguide Modeling of Flared Acoustical TubesSkip 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 00000001 Digital Waveguide Modeling of Flared Acoustical Tubes Maarten van Walstijn (1) Vesa Vilimaki (2) (1) Faculty of Music, University of Edinburgh maarten@music. ed. ac.uk, http://www.ph.ed.ac.uk/fluids/acoustics/ (2) Helsinki University of Technology, Laboratory of Acoustics and Audio Signal Processing email@example.com, http://www.hut.fi/HUT/Acoustics/ Abstract The distributed reflection of the flared tube is calculated using the Schroedinger equation, including the effects ofviscothermal and open end radiation losses. The reflection is approximated with digital filters, which are designed to be efficient and accurate in a digital waveguide model of a wind instrument. Weighted Least-Squares design techniques are applied to model the reflectance of a Bessel horn with an IIR filter. 1 Introduction Digital Waveguide Filters (WGF) are widely used for simulation of wave propagation in musical instruments. The particular efficiency of WGF implementations is based on the separation between unfiltered propagation (digitally represented with delay lines), and additional physical phenomena, such as reflections, losses and dispersion (approximated with digital filters). WGF models have been successfully applied for wind instruments, where the main bore of the instrument is often represented by a single elemental tube shape (cylindrical, conical). Flared tube sections (bell, mouthpiece) are usually approximated with a series of small elemental segments. (see, f.e., , ). In the case of modeling a strongly flared tube (such as brass horns), using piecewise cylindrical segments often requires a large number of segments to enable sufficiently accurate results . A much smaller number of conical segments would be enough, but unfortunately the time-domain filter representations of conical junctions have structural inherent instabilities . Alternatively, the complete reflection and transmission characteristics of the horn may be calculated by applying Schroedinger's equation , and approximated with a single digital filter. This paper addresses the design of such a filter in a WGF Model. In section 2, the use of Schroedinger's equation is discussed, and a formulation with scattering and propagation matrices is presented. Further extensions are made in section 3 by incorporating the effects of viscothermal losses end open-end radiation. The design of an appropriate digital filter is discussed in section 4. Under the assumption that frequencies above the horn cut-off are mainly transmitted out of the horn, we aim to obtain a digital filter that has its main discrepancy in the phase response at the higher frequencies. The phase of the highpass transmitted wave energy is considered to be of minor perceptual importance, whereas the phase of the lowpass reflectance is essential to the sound generation mechanism in wind instruments, and should therefore be approximated accurately. The idea of 'guiding' the main discrepancy into the phase of the higher frequencies is to obtain a horn waveguide structure containing low-order digital filters that will be capable of sufficiently accurate real-time simulation. With these considerations in mind, both FIR and IIR filter design techniques were applied to approximate the analogue filter characteristics. As an example we applied our calculations to a WGF model of a Bessel horn, which is used in a trumpet simulation by Federici and Borin . 2 Schroedinger's Equation Wave propagation in tubes of varying cross-section is modeled by Websters's equation: 1 2p d2_ p 1 dS dp c 8t2 z2 S dz dz (1) where z is the coordinate that corresponds to the position where the spherical wavefront surface S cuts the tube axis (see figure 1). In  a convenient coordinate transformation is used to rewrite (1) into the Schroedinger equation: +2( +(k-F,(z))p= 0 dz (2)
Page 00000002 Here ip = p-S is an alternative wave parameter, k = w/c is the free space wave number, and Fh(Z) = o2 8r/2 is the 'horn function'. + + Pi \ P2 -.................... S....... Pi A-P P2 - 2z Figure 1: Wave propagation in the flared tube In order to find the filter characteristics of the wave propagation-path of p-waves in flared horns, we assume that the horn acts as a one-dimensional waveguide: i.e., p is a sum of right- and left-going waves: Where q0 = qN+1 = k. The horn reflection coefficient is computed as: Rh (w) = G/E. 3 Incorporation of Losses The effects of viscothermal losses can be included by using a 'lossy' wave number that varies with angular frequency c: k(wc) = (o / v()) - ja(C)) (8) with sound velocity v(o) and attenuation factor a(w), for which appropriate values can be found in most literature on acoustics (see, f.e., ). The losses due to the open-end radiation can be represented with scattering (termination) matrix T: p = lp + +p- = Ae-jqz + Beqz (3) 1 1 -R(w) +R() R() +2R( 1 + R(o) R(wo ) 1+ 2R(co) (9) All attenuation and dispersion effects in the horn are expressed in the locally defined propagation constant q(z), which is found by substituting (3) into (2): q(z) = k2 - Fh(z) (4) According to (4), the characteristics of lossless wave propagation are spatially varying (as opposed to wave propagation in cylindrical and conical tubes), and are determined by the horn function. Since this function can be of arbitrary shape, no exact solutions of (2) can be found, and spatial discretization of the wave propagation-path is necessary. Assuming that at junctions between tube sections the pressure is constant and the total sum of incoming volume velocity flow equals zero, we found the following scattering matrices: A suitable expression for the open-end reflection coefficient R(o) is given for spherical wavefronts by Causse et al. . The final EFGH-matrix is now computed by F T. 4 Digital Filter Approximation It is necessary to approximate the reflection function by a low-order digital filter to render it feasible for realtime sound synthesis. It would be easiest to use an FIR filter, but since the impulse response corresponding to the reflectance of a flared tube is usually quite long, a large number of FIR filter coefficients would be needed and this would lead to an inefficient implementation. IIR filters are more complicated to design but they allow effective realization of long impulse responses with a small number of arithmetic operations. An IIR filter design method that allows separate weighting of specified magnitude and phase (group delay) has been developed by Deczky . The method optimizes the poles (rp,i,0p,i) and zeros (ro,,i0o,i0) of a cascade of second-order sections: N 1 - 2ro cos( oi)z + ro 1 H(z) = kON1- oi oi s' - (10) I - 2rpcos( Op)z + rpi 2 Iterative gradient-based techniques are applied to minimize a weighted sum of least-p errors over a discrete set of frequencies cok. The main disadvantage of this method is that good initial pole positions are required to avoid the algorithm from halting at local minima . In practice, we found that restricting the parameter optimization entirely to either magnitude ejqln 0 1 + qn+l 1- Iq Pn =1 J0 n =1 q qI 0 e- jql ' - qn+_1 + q+ L qn n (5) where Pn is the propagation matrix that represents the traveling of waves in section n of length In, and Jn is the scattering matrix for the junction (n,n + 1). Note that P, represents filtered wave propagation, and can not be easily simulated in the digital domain. The final matrix F that relates p-waves on the right side of the tube to ip -waves on the left side, is calculated as a product of N cascading sections: E F=F=Jo'Pl'Jl --.'--Pn'-J, --PNJN (7) G H r
Page 00000003 response or group delay strongly simplifies the choice of the initial parameters and improves the performance of the method. Hence a more robust procedure is to approximate the specified magnitude As and group delay Ts with two separate filters. Once a suitable magnitude approximation filter Hm is found, the desired group delay rd for the second filter becomes: Td(0) = Ts() - Tm(0) (11) where 0q = CkT, with sampling period T. In order to ensure that the second filter does not boost or attenuate any frequencies, (10) is modified into the allpass filter: N r - 2r cos(i)z' + Z2 H(z)=Iil - 2ri cos(i)z' + ri2-2 (12) The following Least-square (LS) error norm was used to determine the coefficients of this filter: K L =,wk[Ta(k)+ TO Td(0k)] k=l Figure 3 shows the impulse response of the horn. The response decays quite slowly and thus an FIR approximation would not be effective. Three secondorder IIR sections were required to approximate the magnitude of the analogue response of the Bessel horn with an average accuracy of 0.25 dB. Here we used equal weighting of a LS error norm in a band of 2,5 KHz, with a sampling frequency of 44.1 KHz (see figure 4a). Using more sections did not improve the accuracy substantially. A larger number of sections (N = 5) is needed for the allpass group delay design (see figure 4b). Fortunately the second-order allpasssections are computationally efficient structures, which can be implemented using just 2 multiplications per section: y(n) = a2,i[x(n) - y(n - 2)] (16) + ai [x(n - ) - y(n - 1)] + x(n -2) with coefficients ali = -2ri cos(oi) and a2,i = r2 The total number of multiplications for this filter is 1+ 3 x 4 + 5 x 2 = 22. This result is comparable to a conical section approximation with similar accuracy. With strategically chosen junction locations, about 10 sections are required for a strongly flared tube , and its WGF implementation would have 10 x 2= 20 multiplications, thereby ignoring the extra fractional delay filters  that will be needed. p + --H---------H.(z) main bore S--Z------JHa (Z) Figure 2: WGF model of wind instrument bore An LS FIR approximation corresponds to truncating the beginning of the impulse response and using that as the filter approximation. In order to approximate the impulse response of the Bessel horn an FIR filter would have about 300 coefficients corresponding to a duration of 6.8 ms. (13) with weighting function Wk. The partial derivatives 8L,I/dr and 8L,/8oi can be expressed in simple closed form and are applied to find the direction towards the minimum error. The stability of the filter is guaranteed by restricting magnitude of each pole (0 < ri < 1). At each iteration the nominal group delay To is optimized and forced to be an (even) integer: / K \ I/K \ To = even wkaW k dk(0k) - Td 0 1 k=1 ) /(=1 / (14) This way the procedure determines the number of delays that might have to be added or subtracted from the delay lines that represent the main bore (see figure 2). The fractionality of the delay is then also incorporated in the allpass filter. To avoid extra filters for the viscothermal losses of the main bore, these may be included into the bell reflectance computation. 5 Results and Implementation The method was applied to a 54cm long Bessel horn with a bore profile described by: r(x) = b(x + a)e (15) For a = 1.05cm, b = 12.7cm,and e = 0.6, the radius r as a function of tube-axis position x forms a good approximation to a trumpet bell. The corresponding horn function is calculated with cubic interpolation . 0 2 4 6 8 10 time (ms) Figure 3: Calculated impulse response of the Bessel horn
Page 00000004 O 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 frequecy (Hz) Figure 4a: Comparison between the ideal magnitude response (solid) and the magnitude approximation filter (dashed) 25o,, I 200o >150 _10 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 frequency (Hz) Figure 4b: Comparison between the desired group delay (solid) and the group delay of the allpass filter (dashed) 6 Conclusion A digital waveguide model of strongly flared acoustical tubes has been developed, including the effects of viscothermal and open-end radiation losses. All these aspects are 'lumped' in our model into two digital filters (one for magnitude and another for phase approximation). The system is designed to have good magnitude characteristics over the entire Nyquist band and good phase properties for those frequencies which are most strongly reflected. Formerly similar models have been realized using cylindrical and conical waveguides and fractional delay filters. The main advantages of the proposed modeling scheme are that (1) the model can be guaranteed to be stable by careful design of the digital filters (whereas conical waveguide filters are unstable in some cases) and (2) separate fractional delay filters that are computationally expensive are not needed. A minor disadvantage is the lack of parametrization since new digital filters need to be designed when the properties of the horn are changed. Fortunately, the geometry of the aircolumn of most musical horns is not changed during playing so the filters can be designed beforehand. Strongly flared tubes, such as the Bessel horn that we used as an example, have a long impulse response that is more efficiently modeled with IIR filters than FIR filters. By separating the magnitude and phase aspects of the IIR design procedure the phase approximation error can be guided into the less important high frequencies. This approach constitutes a structural improvement of the balance between accuracy, efficiency, and perceptual relevance of WGF models of horns and wind instrument bells. References 1] Cook, P.R. 1988. "Implementation of Single Reed Instruments with Arbitrary Bore Shapes using Digital Waveguide Filters." CCRMA Tech. Report no. STAN-M-50, Stanford University.  Vailimaki, V., and M. Karjalainen. 1994. "Digital Waveguide Modeling of Wind Instrument Bores of Truncated Cones." Proc. ICMC, Aarhus, pp. 423-430.  Causse, R., J. Kergomard, and X. Lurton. 1984. "Input impedance of brass musical instruments - Comparison between Experiment and Numerical Models." J. Acoust. Soc. Am., Vol. 75, no. 1, pp. 241-254.  Berners, D.P., and J.O. Smith. 1994. "On the Use of Schrodinger's Equation in the Analytic Determination of Horn Reflectance." Proc. ICMC. pp. 419-422  Federici, R., and Borin, G. 1997."Synthesis of the Trumpet Tone Based on Physical Models." Proc. ICMC, Thessaloniki.  Benade, A.H., and E.V. Jansson. 1974. "On Plane and Spherical Waves in Tubes with Nonuniform Flare I. Theory of Radiation, Resonance Frequencies, and Mode Conversion." Acustica Vol. 31, pp. 80-89  Deczky, A.G. 1972. "Synthesis of Recursive Digital Filters Using the Minimum p-Error Criterion." IEEE Transactions on Audio and Electroacoustics, Vol. 20, no. 4, pp. 257-263.  Smith, J.O. 1983. "Techniques for Digital Filter Design and System Identification, with Application to the Violin." PhD dissertation, CCRMA Tech. Report no. STAN-M-14, Stanford University.  Fletcher, N. H., and T. D. Rossing. 1991. "The Physics of Musical Instruments." New York: Springer Verlag.