Page  384 ï~~Losses in Tubular Acoustic Systems - Theory and Experiment in the Sampled Time and Frequency Domains Noam Amir Giora Rosenhouse EE Dept. Civil Eng. Dept. Technion-IIT Technion-UT Haifa Haifa Israel 32000 Israel 32000 ii Ui Shimony EE Dept. Technion-IIT Haifa Israel 32000 Itai Shelef 49 Hagvura st. Givataim Israel Abstract In this paper we present a discrete model for viscous and thermal losses in a cylindrical acoustic tube. Though good approximations for the frequency response of such a system exist for the continuous case, we show that transforming these approximations to the discrete time domain is not trivial. We show how the discrete impulse response can be calculated theoretically using methods derived from digital filter design, and compare the theoretical calculations to experimental results, with good agreement. Finally, we show how to incorporate this model into a discrete model of a tube with varying cross section, defining a "lossy digital waveguide filter". 1. Introduction The problem of computing the lossy frequency response of an acoustic tube has been dealt with extensively. [Keefe, 1984] gives a good summary of the various approximations for different ratios between radius and wavelength. All of the expressions given by Keefe are valid for the continuous frequency domain. These are complicated functions of frequency though, and therefore it is difficult to inverse-Fourier transform them into the continuous time domain. Furthermore, modem experimental apparatus can deal only with discrete-time measurements. Performing a DFT (Discrete Fourier Transform) on experimental results in order to compare them directly to the continuous-frequency approximations gives seemingly good results, but may be misleading, as will be shown later. This is due to the mathematical nature of transforms from the continuous time and frequency domains to the discrete time and frequency domains. Therefore, performing a precise comparison between theoretical and experimental data requires a more rigorous approach. 2. Continuous frequency calculations A cylindrical acoustic tube can be characterized by the complex wavenumber 'y which determines how a harmonic wave having a planar front and radian frequency W will propagate down the tube. Let us define T as: (1) y=a+ jo/vp, so that a is real, representing the attenuation, and vp is the phase velocity, both of them being functions of frequency. Keefe[ 1984] gives approximations for both of these values as a truncated expansion in r = a(po)/11)112, where a is the radius of the tube, p is the density of air, and i1 is the shear viscosity coef ficient. Let us now define the reflection coefficient F(x) at any point in a tube as the ratio between the pressure waves propagating to the left and to the right. For a tube of length 1 with its left end at the origin and a known termination condition at its right end, we can compute the reflection coefficient at the origin. For instance, assuming a totally reflecting termination and a semi infinite tube to the left, we obtain: (2) r(0) = e2V. From a signal processing viewpoint, this is simply the frequency response of the tube, if we assume that the input is the wave propagating to the right, at the origin, and the output is the wave propagating to the left, at the origin. It follows that the impulse response of the tube is the inverse Fourier transform of r(0) = e-2', though using the series expansions for a and ve it is clearly a formidable task to compute this transform analytically. 3. Discrete time and frequency calculations If we try to carry over the results presented in the previous section to the sampled time and frequency domain, which are the only domains we can refer to in numeric and experimental contexts, we run into difficulties. These difficulties may not be obvious at first glance. A paper by Watkinson et aL [1982], for instance, describes an experiment measuring the pressure and volume velocity at the input to the tube, from which a and ve are calculated. The values obtained show a close correspondence to the theoretically computed values. A similar experiment was performed by the authors, comparing the magnitude and phase of r(0)to its theoretical values. However, we note that the phase is composed of a dominant linear term, which merely represents the 4P.05 384 ICMC Proceedings 1993

Page  385 ï~~delay created by the length of the tube. If we subtract the linear term from the experimental and theoretical graph, the results are quite dissimilar, as shown in Fig. 1. C 10 cc 0o.r a.).15 0 -.15 -.30 -.45 -.60 0 2000 4000 Frequency [Hz] Fig. 1: Discrepancy between theory and experiment The reason behind the discrepancy between the theoretical and experimental results is that the shift from the continuous domain to the discrete domain has been overlooked. Indeed there are many ways to do this, but not all of them are appropriate for the case at hand. This shift in our point of view must in fact be dealt with correctly, and it should be emphasized that the correct formulation of the problem must take into account the details of the experimental apparatus used. When seeking the theoretical value of the frequency response in the discrete domain, we are actually asking what digital filter, when excited by the samples of the pulse used in our experiment, will give the same response as the samples of the response measured in our experiment. This is a known method of designing a digital filter based on a given analog filter, called the input invariance method. This method describes exactly the comparison made between the theoretical and experimental results. Denoting the frequency response of the digital filter (as a function of digital frequency 0) as H(e'Â~), the frequency response of the analog filter as a function of continuous frequency H(jw), and the Fourier transform of the excitation as X(jo), it can be shown that: (3) HeJ) =Z{H(jw). X(jw)} H)0z{xo)} where the operator Z on a Fourier transform denotes performing the inverse transform, sampling, and performing a Fourier transform on the sampled data. The dependence of H(eie), or its inverse Fourier transform h(n), on the specific excitation used in the experiment is now apparent. On examining eq. (3) it also becomes apparent that we have run into numerical difficulties. Computing the inverse Fourier transform of H(jo) is difficult to carry out even in the simple case of the cylindrical tube - not to mention the inverse Fourier transform of X(jwo), for which we will not usually have an analytical expression. 3.1 Finding the equivalent digital filter In this section we demonstrate how to overcome the obstacles presented above. We shall assume in all our computations that a sampling rate of 12 kHz is used, corresponding to a signal bandwidth of 6 kHz, since this was the case in our experimental setup. Method 1- Finite bandwidth input invariance. Assuming that X(j0o) is bandlimited to 6kHz, as it is in our experimental setup, it follows that H(jco)-X( jo) is also bandlimited to 6kHz. Therefore, assuming a sampling rate of 12kHz or more, we obtain: (4) Z{H j)."X(o )} = Z{H(jo)}. Z{X(jo)} where H(jo) is simply H(j)lowpassed to 6kHz. Substituting in (3) and reducing gives us: (5) H(ey) = Z{H(jto)} This is an interesting result in itself. It shows that when using the input invariance filter design method with a band limited input, only the bandwidth of the input affects the shape of the filter. Using this method results in a phase that is discontinuous, hence we still have a problem of ripple to deal with, which is the result of lowpass filtering H(jw) with an ideal filter. Method 2 - The all pole model. In order to achieve a better result in the discrete time and frequency domain, we can attempt to find a model that will give a close fit to the results of method 1, without the ripple, which is the byproduct of moving between continuous and discrete frequency. One possible model is the approximation of H(eje) by an all pole transfer function. Using an autoregressive method [Makhoul, 1975], we approximated H(ei) with a 40 pole transfer function/H(eJÂ~). Stability of H(e1Â~) is ensured by the using the autocorrelation method [Makhoul, 1975]. The results in time and frequency domain appear in Fig. 2. It is apparent that the ripple in the time domain has disappeared, and the discontinuity in the phase also. CL a Time [seconds] Frequency [Hz] Fig. 2: Comparison of the filter to exp. measurements ICMC Proceedings 1993 385 4P.05

Page  386 ï~~From The comparison of the results of the allpole model to the experimental results in Fig. 2, it is apparent that the phase ofH(e ) is close to the experimentally obtained phase, as are the magnitudes and the impulse responses. A point that should be mentioned with respect to the allpole model is that though it gives good results, it is a minimum phase model, meaning that the delay caused by the tube is missing - the peak of the impulse response is at the origin, as is the linear term from the phase of the frequency response. This delay can be recreated by adding the linear phase term to the phase of the frequency response predicted by this model. 4. The complete instrument Using the model described here for losses in a cylindrical tube, it is possible to incorporate it into a discrete model for an instrument with varying cross section. One lossless model for such, an instrument is the cylindrical segment model, which can be implemented as a filter called the "digital waveguide filter" [Smith, 1987]. Such a filter is composed of delays representing the cylindrical segments themselves, and of scattering junctions representing the junction between each two cylindrical segments. In order to account for losses, the delays in this filter can be replaced with the filters h(n) representing losses in each of the segments. An example of such a filter appears in Fig. 3, a filter we chose to call a "lossy digital waveguide filter". plasticine. The microphone is connected to a DT2801 sampling card with 12 bit resolution, in a PC. both the source tube and the tube under measurement must be long enough to ensure that the excitation pulse, the pulse reflected from the termination, and the further reflections from the source loudspeaker do not overlap over the microphone. The experiment was repeated 100 times in a few seconds, and then summed to average out noise. The response was deconvolved with the excitation, in the frequency domain, in order to obtain the impulse response. This deconvolution also removes any linear distortions introduced by the microphone and anti-aliasing filter, ensuring good results. 6. Summary and conclusions Our first point is that comparisons between theoretical results in continuous time and frequency must be transformed correctly before being compared to sampled resuits obtained from measurements. We show that for the case of an acoustic tube, this can be done by using the input invariance method of digital filter design. Though this method of filter design is dependent on the specific input used, we show that if it is band limited we only need to know its bandwidth. This result is important, since it enables to use the input invariance method despite the fact that we do not know the exact shape of the input in continuous time. We show that using an exact model for transforming from continuous to discrete time, an annoying ripple creeps into the results, due to windowing in the continuous frequency domain. Using an all-pole model for the filter in question, we show how this ripple can be removed, finally achieving a digital model for the tube that gives values close to the measured results. Applications of the attenuation model suggested here are numerous. An example given here is the "lossy digital waveguide filter", describing a discrete model of a lossy tube with varying cross section. References [Keefe, 1984] D. H. Keefe, Acoustical wave propagation in cylindrical ducts. JASA, 75(1): pp.58-62, 1984. [Makhoul, 1975] J. Makhoul, Linear prediction: a tutorial review. Proc. of the IEEE, 63(4): pp.561-580, 1975. [Smith, 1987] J. Smith, Waveguide filter tutorial, ICMC conference proceedings, pp.9-16, 1987. [Watkinson et al., 1982] P. S. Watkinson, R. Shepherd and J. M. Bowsher, Acoustic energy losses in brass instruments, Acustica, 51: pp.213-221. 1982. [Watson, 1989] A. P. Watson, Impulse measurements on tubular acoustic systems, Ph.D. Thesis, University of Surrey, 1989. I 1-r 1-r2 junction #1 junction #2 Fig. 3: A lossy digital waveguide filter This filter can be useful for synthesizing musical instruments based on physical modeling, or for other applications in researching the behavior of wind instruments. 5. The experimental apparatus The apparatus used here is very similar to that used by Watson[1989]. A dome tweeter is connected by an airtight adaptor to an aluminum tube. This tube is called the source tube. At its end it is connected to the tube being measured, having the same internal radius, with an adaptor which has a condenser microphone inserted flush with the tube wall. The second tube is of length 2.005 meters, and terminated with a lump of 4P.05 386 ICMC Proceedings 1993