Free Books

Non-Cylindrical Acoustic Tubes

In many situations, the wave impedance of the medium varies in a continuous manner rather than in discrete steps. This occurs, for example, in conical bores and flaring horns. In this section, based on [436], we will consider non-cylindrical acoustic tubes.

Horns as Waveguides

Waves in a horn can be analyzed as one-parameter waves, meaning that it is assumed that a constant-phase wavefront progresses uniformly along the horn. Each ``surface of constant phase'' composing the traveling wave has tangent planes normal to the horn axis and to the horn boundary. For cylindrical tubes, the surfaces of constant phase are planar, while for conical tubes, they are spherical [357,317,144]. The key property of a ``horn'' is that a traveling wave can propagate from one end to the other with negligible ``backscattering'' of the wave. Rather, it is smoothly ``guided'' from one end to the other. This is the meaning of saying that a horn is a ``waveguide''. The absence of backscattering means that the entire propagation path may be simulated using a pure delay line, which is very efficient computationally. Any losses, dispersion, or amplitude change due to horn radius variation (``spreading loss'') can be implemented where the wave exits the delay line to interact with other components.

Overview of Methods

We will first consider the elementary case of a conical acoustic tube. All smooth horns reduce to the conical case over sufficiently short distances, and the use of many conical sections, connected via scattering junctions, is often used to model an arbitrary bore shape [71]. The conical case is also important because it is essentially the most general shape in which there are exact traveling-wave solutions (spherical waves) [357].

Beyond conical bore shapes, however, one can use a Sturm-Liouville formulation to solve for one-parameter-wave scattering parameters [50]. In this formulation, the curvature of the bore's cross-section (more precisely, the curvature of the one-parameter wave's constant-phase surface area) is treated as a potential function that varies along the horn axis, and the solution is an eigenfunction of this potential. Sturm-Liouville analysis is well known in quantum mechanics for solving elastic scattering problems and for finding the wave functions (at various energy levels) for an electron in a nonuniform potential well.

When the one-parameter-wave assumption breaks down, and multiple acoustic modes are excited, the boundary element method (BEM) is an effective tool [190]. The BEM computes the acoustic field from velocity data along any enclosing surface. There also exist numerical methods for simulating wave propagation in varying cross-sections that include ``mode conversion'' [336,13,117].

This section will be henceforth concerned with non-cylindrical tubes in which backscatter and mode-conversion can be neglected, as treated in [317]. The only exact case is the cone, but smoothly varying horn shapes can be modeled approximately in this way.

Back to the Cone

Note that the cylindrical tube is a limiting case of a cone with its apex at infinity. Correspondingly, a plane wave is a limiting case of a spherical wave having infinite radius.

On a fundamental level, all pressure waves in 3D space are composed of spherical waves [357]. You may have learned about the Huygens-Fresnel principle in a physics class when it covered waves [295]. The Huygens-Fresnel principle states that the propagation of any wavefront can be modeled as the superposition of spherical waves emanating from all points along the wavefront [122, page 344]. This principle is especially valuable for intuitively understanding diffraction and related phenomena such as mode conversion (which happens, for example, when a plane wave in a horn hits a sharp bend or obstruction and breaks up into other kinds of waves in the horn).

Conical Acoustic Tubes

The conical acoustic tube is a one-dimensional waveguide which propagates circular sections of spherical pressure waves in place of the plane wave which traverses a cylindrical acoustic tube [22,349]. The wave equation in the spherically symmetric case is given by

$\displaystyle c^2 p''_x = \ddot{p}_x \protect$ (C.144)


c & \isdef & \mbox{sound speed} & \qqu...
...'_x & \isdef & \dfrac{\partial}{\partial x}p_x(t,x)

and $ p(t,x)$ is the pressure at time $ t$ and radial position $ x$ along the cone axis (or wall). In terms of $ p$ (rather than $ p_x=xp$), Eq.$ \,$(C.145) expands to Webster's horn equation [357]:

$\displaystyle \frac{1}{c^2} \ddot{p}
\eqsp \frac{1}{A}\left(A p'\right)'
\eqsp p'' + \frac{A'}{A}p'

where $ A=\alpha x^2$ is the area of the spherical wavefront, so that $ A'/A=2/x$ (as discussed further in §C.18.4 below).

Spherical coordinates are appropriate because simple closed-form solutions to the wave equation are only possible when the forced boundary conditions lie along coordinate planes. In the case of a cone, the boundary conditions lie along a conical section of a sphere. It can be seen that the wave equation in a cone is identical to the wave equation in a cylinder, except that $ p$ is replaced by $ xp$. Thus, the solution is a superposition of left- and right-going traveling wave components, scaled by $ 1/x$:

$\displaystyle p(t,x) = \frac{ f\left(t-\frac{x}{c}\right)}{x} + \frac{ g\left(t+\frac{x}{c} \right)}{x} \protect$ (C.145)

where $ f(\cdot)$ and $ g(\cdot)$ are arbitrary twice-differentiable continuous functions that are made specific by the boundary conditions. A function of $ (t-x/c)$ may be interpreted as a fixed waveshape traveling to the right, (i.e., in the positive $ x$ direction), with speed $ c$. Similarly, a function of $ (t+x/c)$ may be seen as a wave traveling to the left (negative $ x$ direction) at $ c$ meters per second. The point $ x=0$ corresponds to the tip of the cone (center of the sphere), and $ p(t,x)$ may be singular there.

In cylindrical tubes, the velocity wave is in phase with the pressure wave. This is not the case with conical or more general tubes. The velocity of a traveling may be computed from the corresponding traveling pressure wave by dividing by the wave impedance.

Digital Simulation

A discrete-time simulation of the above solution may be obtained by simply sampling the traveling-wave amplitude at intervals of $ T$ seconds, which implies a spatial sampling interval of $ X\isdeftext cT$ meters. Sampling is carried out mathematically by the change of variables

x& \to & x_m&=& x_0+ mX\\
t & \to & t_n&=& nT

and this substitution in Eq.$ \,$(C.146) gives

x_mp(t_n,x_m) &\,\mathrel{\mathop=}\,& f(t_n- x_m/c)+g(t_n+ x...
...hop=}\,& f\left[(n-m)T-x_0/c\right]+ g\left[(n+m)T+x_0/c\right].


$\displaystyle p^+(n) \isdef f(nT-x_0/c) \qquad\qquad p^-(n) \isdef g(nT+x_0/c)

where $ x_0$ is arbitrarily chosen as the position along the cone closest to the tip. (We avoid a sample at the tip itself, where a pressure singularity may exist.) Then a section of the ideal cone can be simulated as shown in Figure C.43 (where pressure outputs are shown for $ x=x_0$ and $ x=x_0+3X$). A particle-velocity output is formed by dividing the traveling pressure waves by the wave impedance of air. Since the wave impedance is now a function of frequency and propagation direction (as derived below), a digital filter will replace what was a real number for cylindrical tubes.

Figure C.43: Digital simulation of the ideal, lossless, conical waveguide with observation points at $ x=x_0$ and $ x=x_0+3X=x_0+3cT$. The symbol ``$ z^{-1}$'' denotes a one-sample delay.

A more compact simulation diagram which stands for either sampled or continuous simulation is shown in Figure C.44. The figure emphasizes that the ideal, lossless waveguide is simulated by a bidirectional delay line.

Figure C.44: Simplified picture of ideal waveguide simulation.

As in the case of uniform waveguides, the digital simulation of the traveling-wave solution to the lossless wave equation in spherical coordinates is exact at the sampling instants, to within numerical precision, provided that the traveling waveshapes are initially bandlimited to less than half the sampling frequency. Also as before, bandlimited interpolation can be used to provide time samples or position samples at points off the simulation grid. Extensions to include losses, such as air absorption and thermal conduction, or dispersion, can be carried out as described in §2.3 and §C.5 for plane-wave propagation (through a uniform wave impedance).

The simulation of Fig.C.44 suffices to simulate an isolated conical frustum, but what if we wish to interconnect two or more conical bores? Even more importantly, what driving-point impedance does a mouthpiece ``see'' when attached to the narrow end of a conical bore? The preceding only considered pressure-wave behavior. We must now also find the velocity wave, and form their ratio to obtain the driving-point impedance of a conical tube.

Momentum Conservation in Nonuniform Tubes

Newton's second law ``force equals mass times acceleration'' implies that the pressure gradient in a gas is proportional to the acceleration of a differential volume element in the gas. Let $ A(x)$ denote the area of the surface of constant phase at radial coordinate $ x$ in the tube. Then the total force acting on the surface due to pressure is $ f(t,x)=A(x)p(t,x)$, as shown in Fig.C.45.

Figure C.45: Differential volume element for the conical acoustic tube.

The net force $ df(t,x)$ to the right across the volume element between $ x$ and $ x+dx$ is then

$\displaystyle df(t,x) = f(t,x)-f(t,x+dx)$ $\displaystyle =$ $\displaystyle f(t,x)-[f(t,x) + dx \cdot f'(x)]$  
  $\displaystyle =$ $\displaystyle - dx \cdot f'(x)$  
  $\displaystyle =$ $\displaystyle - dx \cdot[A(x)p(t,x)]'$  
  $\displaystyle =$ $\displaystyle -dx \cdot[A' p + A'p],$  

where, when time and/or position arguments have been dropped, as in the last line above, they are all understood to be $ t$ and $ x$, respectively. To apply Newton's second law equating net force to mass times acceleration, we need the mass of the volume element

dM(t,x) &=& \int_x^{x+dx} \rho(t,\xi) A(\xi)\,d\xi \\ [5pt]
...\rho A' \right)\frac{(dx)^2}{2}\\ [5pt]
&\approx& \rho\, A\,dx,

where $ \rho(t,x)$ denotes air density.

The center-of-mass acceleration of the volume element can be written as $ {\dot u}(t,x)$ where $ u(t,x)$ is particle velocity.C.16 Applying Newton's second law $ df = dM\cdot {\dot u}$, we obtain

$\displaystyle -dx \cdot (A'p + Ap') \eqsp \rho\, A\,dx\, {\dot u} \protect$ (C.146)

or, dividing through by $ -A\,dx$,

$\displaystyle p' + p \, \frac{A'}{A} \eqsp - \rho\,{\dot u}. \protect$ (C.147)

In terms of the logarithmic derivative of $ A$, this can be written

$\displaystyle p' + p \,$   ln$\displaystyle 'A \eqsp - \rho\,{\dot u}. \protect$ (C.148)

Note that $ p$ denotes small-signal acoustic pressure, while $ \rho$ denotes the full gas density (not just an acoustic perturbation in the density). We may therefore treat $ \rho$ as a constant.

Cylindrical Tubes

In the case of cylindrical tubes, the logarithmic derivative of the area variation, ln$ 'A(x) = A'/A$, is zero, and Eq.$ \,$(C.148) reduces to the usual momentum conservation equation $ p' = -\rho {\dot u}$ encountered when deriving the wave equation for plane waves [318,349,47]. The present case reduces to the cylindrical case when

$\displaystyle \frac{A'}{A} \;\ll\; \frac{p'}{p}

i.e., when the relative change in cross-sectional area is much less than the relative change in pressure along the tube. In other words, the tube area variation must be slower than the spatial variation of the wave itself. This assumption is also necessary for the ``one-parameter-wave'' approximation to hold in the first place.

If we look at sinusoidal spatial waves, $ p=A_p e^{j k_p x}$ and $ A=A_A
e^{j k_A x}$, then $ A'/A = k_A$ and $ p'/p = k_p$, and the condition for cylindrical-wave behavior becomes $ k_A\ll k_p$, i.e., the spatial frequency of the wall variation must be much less than that of the wave. Another way to say this is that the wall must be approximately flat across a wavelength. This is true for smooth horns/bores at sufficiently high wave frequencies.

Wave Impedance in a Cone

From Eq.$ \,$(C.146) we have that the traveling-wave solution of the wave equation in spherical coordinates can be expressed as

$\displaystyle p(t,x) \eqsp \frac{f\left(t \mp \frac{x}{c}\right)}{x}

where the minus sign is associated with an expanding spherical wave, and the plus sign corresponds to a converging wave. The spatial derivative of this function is

$\displaystyle p' \eqsp \frac{f'}{x}-\frac{f}{x^2} \eqsp \mp \frac{{\dot f}}{c x}-\frac{f}{x^2} \eqsp \mp \frac{{\dot p}}{c}-\frac{p}{x} \protect$ (C.149)

i.e., it can be expressed in terms of its own time derivative. This is a general property of any traveling wave.

Figure C.46: Cone parameters.

Referring to Fig.C.46, the area function $ A(x)$ can be written for any cone in terms of the distance from its apex as

$\displaystyle A(x) \eqsp \alpha x^2

for some $ \alpha$. (We chose $ x=0$ at the tip of the cone in arriving at the basic traveling-wave solution to the wave equation.) This formula holds for both the planar cross-section indicated in Fig.C.46 and the spherical section having the same perimeter (i.e., spherical surface area is proportional to radius squared). The logarithmic derivative of $ A(x)=\alpha x^2$ is then

   ln$\displaystyle ' A \eqsp \frac{ A' }{A} \eqsp \frac{2}{x},

which is independent of the conical taper angle $ \theta =
\tan\left(\sqrt{\alpha/\pi}\right)$. That is, one conical section of a spherical wave is like any other, as it must be due to spherical symmetry.

Substituting the logarithmic derivative of $ A(x)$ and $ p'$ from Eq.$ \,$(C.150) into the momentum-conservation equation Eq. (C.148) yields

$\displaystyle \mp \frac{{\dot p}}{c}-\frac{p}{x} + p\frac{2}{x} \eqsp - \rho {\dot u}


$\displaystyle \mp \frac{{\dot p}}{c} + \frac{p}{x} \eqsp - \rho {\dot u},

where the minus sign is for an expanding spherical wave, and the plus sign is for a converging spherical wave. Taking the Laplace transform of the above expression gives

$\displaystyle \mp \frac{s P(s)}{c} + \frac{P(s)}{x} \eqsp - \rho s U(s)

in the case of zero initial conditions $ p(0,x)=u(0,x)=0$ for all $ x$.

We can now solve for the wave impedance in each direction, where the wave impedance may be defined (§7.1) as the Laplace transform of the traveling pressure divided by the Laplace transform of the corresponding traveling velocity wave:

$\displaystyle R^{+}(s)$ $\displaystyle \isdef$ $\displaystyle \frac{P^+(s)}{U^{+}(s)}$  
$\displaystyle R^{-}(s)$ $\displaystyle \isdef$ $\displaystyle \frac{P^-(s)}{U^{-}(s)}$  

We introduce the shorthand

$\displaystyle P^\pm(s)\eqsp \pm R^\pm (s)U^\pm(s)

where all of the upper signs or all of the lower signs are taken together, corresponding to expanding or converging spherical waves, respectively. In this notation, we may solve for the conical wave impedance as

$\displaystyle R^\pm (s)\isdefs \frac{P^\pm(s)}{\pm U^\pm(s)} \eqsp \frac{\rho c}{1 \mp \frac{c}{sx}}.

Along the frequency axis $ s=j\omega$ we get

$\displaystyle R^\pm (j\omega) \eqsp \frac{\rho c}{1 \pm j\frac{c}{\omega x}}
\eqsp \frac{\rho c}{1 \pm j\frac{1}{kx}}

where $ \omega = 2\pi f$ is radian temporal frequency, and $ k \isdef
\omega/c$ is the radian spatial frequency, or wavenumber.

Note that for a cylindrical tube, the wave impedance in both directions is $ R^\pm (s)= \rho c$, and there is no frequency dependence. A wavelength or more away from the conical tip, i.e., for $ x\gg 1/k = \lambda/2\pi$, where $ \lambda$ is the spatial wavelength, the wave impedance again approaches that of a cylindrical bore. However, in conical musical instruments, the fundamental wavelength is typically twice the bore length, so the complex nature of the wave impedance is important throughout the bore and approaches being purely imaginary near the mouthpiece. This is especially relevant to conical-bore double-reeds, such as the bassoon.

Writing the wave impedance as

$\displaystyle R^\pm (s)\isdefs \frac{P^\pm(s)}{\pm U^\pm(s)}
\eqsp \frac{1}{\frac{1}{\rho c} \mp \frac{1}{s\rho x}},

the reader familiar with circuit theory (§7.2.3) will recognize this as the parallel combination of the wave impedance $ \rho c$ of a cylindrical bore and the lumped impedance of a mass formally equal to $ \mp\rho x$. This equivalent circuit is explored in [37]. As $ x$ goes to 0 at the tip, the mass goes to zero and ``shorts out'' the traveling wave. Note that in the case of a converging cone, the mass is negative. This will give rise to an unstable one-pole filter in the waveguide model, and various solutions have been proposed for this problem [528,50].

Up to now, we have been defining wave impedance as pressure divided by particle velocity. In acoustic tubes, volume velocity is what is conserved at a junction between two different acoustic tube types. Therefore, in acoustic tubes, we define the wave impedance as the ratio of pressure to volume velocity

$\displaystyle R_A^\pm (s)\isdefs \frac{P^\pm(s)}{\pm U^\pm(s)A(x)} \eqsp \frac{R^\pm (s)}{A(x)} \eqsp
\frac{\rho c}{A(x)} \frac{1}{1 \pm \frac{c}{sx}}

or, defining $ t_x \isdef x/c$ as the time to propagate over the distance $ x$,

$\displaystyle \fbox{$\displaystyle R_A^\pm (s)\eqsp \frac{\rho c}{A(x)} \frac{s}{s \pm \frac{1}{t_x}}.$} \protect$ (C.150)

This is the wave impedance we use to compute the generalized reflection and transmission coefficients at a change in cross-sectional area and/or taper angle in a conical acoustic tube. Note that it has a zero at $ s=0$ and a pole at $ s=\pm 1/t_x$.

In this case, the equivalent mass is $ \pm\rho x/A(x) = \pm\rho/(\alpha
x)$. It would perhaps be more satisfying if the equivalent mass in the conical wave impedance were instead $ \pm\rho x A(x)$ which is the mass of air contained in a cylinder of radius $ A(x)$ projected back to the tip of the cone. However, the ``acoustic mass'' cannot be physically equivalent to mechanical mass. To see this, consider that the impedance of a mechanical mass is $ ms$ which is in physical units of mass per unit time, and by definition of mechanical impedance this equals force over velocity. The impedance in an acoustic tube, on the other hand, must be in units of pressure (force/area) divided by volume velocity (velocity $ \times$ area) and this reduces to

$\displaystyle \frac{\mbox{force/area}}{\mbox{velocity}\cdot\mbox{area}}
\eqsp \...
\eqsp \frac{\mbox{mass-volume-density}}{\mbox{time}\cdot\mbox{distance}}

which is what we found.

The real part of the wave impedance corresponds to transportation of wave energy, the imaginary part is a so-called ``reactance'' and does not correspond to power transfer. Instead, it corresponds to a ``standing wave'' which is created by equal and opposite power flow, or an ``evanescent wave'' (§C.8.2), which is a non-propagating, exponentially decaying, limiting form of a traveling wave in which the ``propagation constant'' is purely imaginary due to being at a frequency above or below a ``cut off'' frequency for the waveguide [295,122]. Driving an ideal mass at the end of a waveguide results in total reflection of all incident wave energy along with a quarter-cycle phase shift. Another interpretation is that the traveling wave becomes a standing wave at the tip of the cone. This is one way to see how the resonances of a cone can be the same as those of a cylinder the same length which is open on both ends. (One might first expect the cone to behave like a cylinder which is open on one end and closed on the other.) Because the impedance approaches a purely imaginary zero at the tip, it looks like a mass (with impedance $ \pm\rho x s$). The ``piston of air'' at the open end similarly looks like a mass [285].

More General One-Parameter Waves

The wave impedance derivation above made use of known properties of waves in cones to arrive at the wave impedances in the two directions of travel in cones. We now consider how this solution might be generalized to arbitrary bore shapes. The momentum conservation equation is already applicable to any wavefront area variation $ A(x)$:

$\displaystyle p' + p \,$   ln$\displaystyle 'A \eqsp -\rho{\dot u}.

As we did for vibrating stringsC.3.4), suppose the pressure is sinusoidally driven so that we have

$\displaystyle p(t,x) \eqsp P(x) e^{st}

where $ s=j\omega$, $ \omega = 2\pi f$, and $ f$ is the driving frequency. The partial derivatives become

$\displaystyle p' \eqsp ($ln$\displaystyle 'P) p, \qquad \qquad \dot{p}\eqsp s p.

Substituting into the momentum equation gives

$\displaystyle p($ln$\displaystyle 'P +$   ln$\displaystyle 'A) \eqsp -\rho {\dot u}.

Because the medium is linear and time-invariant, the velocity $ u(t,x)$ must be of the form $ U(x) e^{st}$, and we can define the spatially instantaneous wave impedance as

$\displaystyle R(x) \isdefs \frac{P(x)}{U(x)}.

The corresponding instantaneous wave admittance is then $ \Gamma(x)\isdef U(x)/P(x)$. Then $ \u=\Gamma(x) p$, and the momentum equation becomes

   ln$\displaystyle 'P +$   ln$\displaystyle 'A \eqsp -\rho s\Gamma

Solving for the wave impedance gives

$\displaystyle R(x) \eqsp -\frac{s\rho}{\mbox{ln}'P + \mbox{ln}'A}.

Expressing $ P(x)$ in exponential form as

$\displaystyle P(x) \isdefs e^{j\theta(x)}

where $ \theta(x)$ may be complex, we may define the instantaneous spatial frequency (wavenumber) as

$\displaystyle k(x) \isdefs \theta'(x)

and since ln$ 'P = j k(x)$, we have

$\displaystyle R(x) \eqsp -\frac{s\rho}{j k(x) + \mbox{ln}'A} \eqsp -\frac{\omega\rho}{k(x) + \mbox{ln}'A}

Defining the spatially instantaneous phase velocity as

$\displaystyle c(x)\isdefs \frac{\omega}{k(x)}

we have

$\displaystyle R(x) \eqsp -\frac{\rho c(x)}{1 + \frac{1}{j k(x)}{\mbox{ln}'A(x)}} \eqsp -\frac{\rho c(x)}{1 + \frac{\mbox{ln}'A(x)}{\mbox{ln}'P(x)}} \protect$ (C.151)

This reduces to the simple case of the uniform waveguide when the logarithmic derivative of cross-sectional area $ A(x)$ is small compared with the logarithmic derivative of the amplitude $ P(x)$ which is proportional to the instantaneous spatial frequency. A traveling wave solution interpretation makes sense when the instantaneous wavenumber $ k(x)$ is approximately real, and the phase velocity $ c(x)$ is approximately constant over a number of wavelengths $ \lambda(x) = 2\pi/k(x)$.

Generalized Wave Impedance

Figure C.47 depicts a section of a conical bore which widens to the right connected to a section which narrows to the right. In addition, the cross-sectional areas are not matched at the junction.

The horizontal $ x$ axis (taken along the boundary of the cone) is chosen so that $ x=0$ corresponds to the apex of the cone. Let $ A(x)=\alpha x^2$ denote the cross-sectional area of the bore.

Figure C.47: a) Physical picture. b) Waveguide implementation.

Since a piecewise-cylindrical approximation to a general acoustic tube can be regarded as a ``zeroth-order hold'' approximation. A piecewise conical approximation then uses first-order (linear) segments. One might expect that quadratic, cubic, etc., would give better and better approximations. However, such a power series expansion has a problem: In zero-order sections (cylinders), plane waves propagate as traveling waves. In first-order sections (conical sections), spherical waves propagate as traveling waves. However, there are no traveling wave types for higher-order waveguide flare (e.g., quadratic or higher) [357].

Since the digital waveguide model for a conic section is no more expensive to implement than that for a cylindrical section, (both are simply bidirectional delay lines), it would seem that modeling accuracy can be greatly improved for non-cylindrical bores (or parts of bores such as the bell) essentially for free. However, while the conic section itself costs nothing extra to implement, the scattering junctions between adjoining cone segments are more expensive computationally than those connecting cylindrical segments. However, the extra expense can be small. Instead of a single, real, reflection coefficient occurring at the interface between two cylinders of differing diameter, we obtain a first-order reflection filter at the interface between two cone sections of differing taper angle, as seen in the next section.

Generalized Scattering Coefficients

Generalizing the scattering coefficients at a multi-tube intersection (§C.12) by replacing the usual real tube wave impedance $ R_i=\rho c/A_i$ by the complex generalized wave impedance

$\displaystyle R(x)=-\frac{\rho c(x)}{1 + \frac{\mbox{ln}'A(x)}{\mbox{ln}'P(x)}}

from Eq.$ \,$(C.152), or, as a special case, the conical-section wave impedance $ R_A^\pm (s)=[\rho c/A(x)]/[s/(s \pm 1/t_x)]$ from Eq.$ \,$(C.151), we obtain the junction-pressure phasor [436]

$\displaystyle P_J = \left(G_J + \sum_{i=1}^N G_i^-\right)^{-1} \sum_{i=1}^N
\left(G_i^+ + G_i^- \right)P_i^+

where $ G_i^+
$ is the complex, frequency-dependent, incoming, acoustic admittance of the $ i$th branch at the junction, $ G_i^-$ is the corresponding outgoing acoustic admittance, $ P_i^+$ is the incoming traveling pressure-wave phasor in branch $ i$, $ P_i^- = P_J -
P_i^+$ is the outgoing wave, and $ G_J$ is the admittance of a load at the junction, such as a coupling to another simulation. For generality, the formula is given as it appears in the multivariable case.

Cylinder with Conical Cap

Consider a cylindrical acoustic tube adjoined to a converging conical cap, as depicted in Figure C.48a. We may consider the cylinder to be either open or closed on the left side, but everywhere else it is closed. Since such a physical system is obviously passive, an interesting test of acoustic theory is to check whether theory predicts passivity in this case.

Figure C.48: Cylindrical acoustic tube with a conical cap. a) Physical outline. b) Continuous-time waveguide formulation.

It is well known that a growing exponential appears at the junction of two conical waveguides when the waves in one conical taper angle reflect from a section with a smaller (or more negative) taper angle [7,300,8,160,9]. The most natural way to model a growing exponential in discrete time is to use an unstable one-pole filter [506]. Since unstable filters do not normally correspond to passive systems, we might at first expect passivity to not be predicted. However, it turns out that all unstable poles are ultimately canceled, and the model is stable after all, as we will see. Unfortunately, as is well known in the field of automatic control, it is not practical to attempt to cancel an unstable pole in a real system, even when it is digital. This is because round-off errors will grow exponentially in the unstable feedback loop and eventually dominate the output.

The need for an unstable filter to model reflection and transmission at a converging conical junction has precluded the use of a straightforward recursive filter model [406]. Using special ``truncated infinite impulse response'' (TIIR) digital filters [540], an unstable recursive filter model can in fact be used in practice [528]. All that is then required is that the infinite-precision system be passive, and this is what we will show in the special case of Fig.C.48.

Scattering Filters at the Cylinder-Cone Junction

As derived in §C.18.4, the wave impedance (for volume velocity) at frequency $ \omega $ rad/sec in a converging cone is given by

$\displaystyle Z_\xi(j\omega) = \frac{\rho c}{S(\xi)} \frac{j\omega}{j\omega-c/\xi}$ (C.152)

where $ \xi$ is the distance to the apex of the cone, $ S(\xi)$ is the cross-sectional area, and $ \rho c$ is the wave impedance in open air. A cylindrical tube is the special case $ \xi=\infty$, giving $ Z_\infty(j\omega) = \rho c/S$, independent of position in the tube. Under normal assumptions such as pressure continuity and flow conservation at the cylinder-cone junction (see, e.g., [300]), the junction reflection transfer function (reflectance) seen from the cylinder looking into the cone is derived to be

$\displaystyle R(s) = -\frac{c/\xi}{c/\xi - 2s}$ (C.153)

(where $ s$ is the Laplace transform variable which generalizes $ s=j\omega$) while the junction transmission transfer function (transmittance) to the right is given by

$\displaystyle T(s) = 1 + R(s) = -\frac{2s}{c/\xi - 2s}$ (C.154)

The reflectance and transmittance from the right of the junction are the same when there is no wavefront area discontinuity at the junction [300]. Both $ R(s)$ and $ T(s)$ are first-order transfer functions: They each have a single real pole at $ s=c/(2\xi)$. Since this pole is in the right-half plane, it corresponds to an unstable one-pole filter.

Reflectance of the Conical Cap

Let $ {t_{\xi}}\isdef \xi/c$ denote the time to propagate across the length of the cone in one direction. As is well known [22], the reflectance at the tip of an infinite cone is $ -1$ for pressure waves. I.e., it reflects like an open-ended cylinder. We ignore any absorption losses propagating in the cone, so that the transfer function from the entrance of the cone to the tip is $ e^{-s{t_{\xi}}}$. Similarly, the transfer function from the conical tip back to the entrance is also $ e^{-s{t_{\xi}}}$. The complete reflection transfer function from the entrance to the tip and back is then

$\displaystyle R_{t_{\xi}}(s) = - e^{-2s{t_{\xi}}}$ (C.155)

Note that this is the reflectance a distance $ \xi=c{t_{\xi}}$ from a conical tip inside the cone.

We now want to interface the conical cap reflectance $ R_{t_{\xi}}(s)$ to the cylinder. Since this entails a change in taper angle, there will be reflection and transmission filtering at the cylinder-cone junction given by Eq.$ \,$(C.154) and Eq.$ \,$(C.155).

From inside the cylinder, immediately next to the cylinder-cone junction shown in Fig.C.48, the reflectance of the conical cap is readily derived from Fig.C.48b and Equations (C.154) and (C.155) to be

$\displaystyle R_J(s)$ $\displaystyle =$ $\displaystyle \frac{R(s) + 2 R(s) R_{t_{\xi}}(s) + R_{t_{\xi}}(s)}{1 - R(s)R_{t...
= \frac{1 + (1+2s{t_{\xi}})R_{t_{\xi}}(s)}{2s{t_{\xi}}-1-R_{t_{\xi}}(s)}$  
  $\displaystyle =$ $\displaystyle \frac{1 - (1+2s{t_{\xi}})e^{-2s{t_{\xi}}}}{2s{t_{\xi}}-1+e^{-2s{t_{\xi}}}}
\isdef \frac{N(s)}{D(s)}$ (C.156)


$\displaystyle N(s) \isdef 1 - e^{-2s{t_{\xi}}} - 2s{t_{\xi}}e^{-2s{t_{\xi}}}$ (C.157)

is the numerator of the conical cap reflectance, and

$\displaystyle D(s) \isdef 2 s{t_{\xi}}- 1 + e^{-2s{t_{\xi}}}$ (C.158)

is the denominator. Note that for very large $ {t_{\xi}}$, the conical cap reflectance approaches $ R_J = -e^{-2s{t_{\xi}}}$ which coincides with the impedance of a length $ \xi=c{t_{\xi}}$ open-end cylinder, as expected.

Stability Proof

A transfer function $ R_J(s) = N(s)/D(s)$ is stable if there are no poles in the right-half $ s$ plane. That is, for each zero $ s_i$ of $ D(s)$, we must have re$ \left\{s_i\right\}\leq 0$. If this can be shown, along with $ \left\vert R_J(j\omega)\right\vert\leq 1$, then the reflectance $ R_J$ is shown to be passive. We must also study the system zeros (roots of $ N(s)$) in order to determine if there are any pole-zero cancellations (common factors in $ D(s)$ and $ N(s)$).

Since re$ \left\{s{t_{\xi}}\right\}\geq 0$ if and only if re$ \left\{s\right\}\geq 0$, for $ {t_{\xi}}>0$, we may set $ {t_{\xi}}=1$ without loss of generality. Thus, we need only study the roots of

N(s) &=& 1 - e^{-2s} - 2s e^{-2s} \\
D(s) &=& 2s - 1 + e^{-2s}

If this system is stable, we have stability also for all $ {t_{\xi}}>0$. Since $ e^{-2s}$ is not a rational function of $ s$, the reflectance $ R_J(s)$ may have infinitely many poles and zeros.

Let's first consider the roots of the denominator

$\displaystyle D(s) = 2s - 1 + e^{-2s}.$ (C.159)

At any solution $ s$ of $ D(s)=0$, we must have

$\displaystyle s = \frac{1-e^{-2s}}{2}$ (C.160)

To obtain separate equations for the real and imaginary parts, set $ s=\sigma+j\omega$, where $ \sigma$ and $ \omega $ are real, and take the real and imaginary parts of Eq.$ \,$(C.161) to get

\mbox{re}\left\{D(s)\right\} &=& (2\sigma - 1) + e^{-2\sigma}\...
...}\left\{D(s)\right\} &=& 2\omega - e^{-2\sigma}\sin(2\omega) = 0

Both of these equations must hold at any pole of the reflectance. For stability, we further require $ \sigma\leq 0$. Defining $ \tau = 2\sigma$ and $ \nu=2\omega$, we obtain the somewhat simpler conditions

$\displaystyle e^\tau ( 1 - \tau)$ $\displaystyle =$ $\displaystyle \cos(\nu)$ (C.161)
$\displaystyle e^\tau$ $\displaystyle =$ $\displaystyle \frac{\sin(\nu)}{\nu}$ (C.162)

For any poles of $ R_J(s)$ on the $ j\omega $ axis, we have $ \tau=0$, and Eq.$ \,$(C.163) reduces to

$\displaystyle \nu = \sin(\nu)$ (C.163)

It is well known that the ``sinc function'' $ \sin(\nu)/\nu$ is less than $ 1$ in magnitude at all $ \nu$ except $ \nu=0$. Therefore, Eq.$ \,$(C.164) can hold only at $ \omega=\nu=0$.

We have so far proved that any poles on the $ j\omega $ axis must be at $ \omega=0$.

The same argument can be extended to the entire right-half plane as follows. Going back to the more general case of Eq.$ \,$(C.163), we have

$\displaystyle \frac{\sin(\nu)}{\nu} = e^\tau$ (C.164)

Since $ \left\vert\sin(\nu)/\nu\right\vert\leq 1$ for all real $ \nu$, and since $ \left\vert e^\tau\right\vert>1$ for $ \tau>0$, this equation clearly has no solutions in the right-half plane. Therefore, the only possible poles in the right-half plane, including the $ j\omega $ axis, are at $ s=0$.

In the left-half plane, there are many potential poles. Equation (C.162) has infinitely many solutions for each $ \tau<0$ since the elementary inequality $ 1-\tau \leq e^{-\tau}$ implies $ e^\tau ( 1
- \tau) < e^\tau e^{-\tau} = 1$. Also, Eq.$ \,$(C.163) has an increasing number of solutions as $ \tau $ grows more and more negative; in the limit of $ \tau=-\infty$, the number of solutions is infinite and given by the roots of $ \sin(\nu)$ ( $ \nu = k\pi$ for any integer $ k$). However, note that at $ \tau=-\infty$, the solutions of Eq.$ \,$(C.162) converge to the roots of $ \cos(\nu)$ ( $ \nu = (2k+1)(\pi/2)$ for any integer $ k$). The only issue is that the solutions of Eq.$ \,$(C.162) and Eq.$ \,$(C.163) must occur together.

Figure: Locus of solutions to Eq.$ \,$(C.162) (dashed) (and Eq.$ \,$(C.163) (solid). The dotted lines are lines of construction for an analytic proof that there are infinitely many roots in the left-half $ s$-plane.

Figure C.49 plots the locus of real-part zeros (solutions to Eq.$ \,$(C.162)) and imaginary-part zeros (Eq.$ \,$(C.163)) in a portion the left-half plane. The roots at $ s=0$ can be seen at the middle-right. Also, the asymptotic interlacing of these loci can be seen along the left edge of the plot. It is clear that the two loci must intersect at infinitely many points in the left-half plane near the intersections indicated in the graph. As $ \left\vert\nu\right\vert$ becomes large, the intersections evidently converge to the peaks of the imaginary-part root locus (a log-sinc function rotated 90 degrees). At all frequencies $ \nu$, the roots occur near the log-sinc peaks, getting closer to the peaks as $ \left\vert\nu\right\vert$ increases. The log-sinc peaks thus provide a reasonable estimate number and distribution in the left-half $ s$-plane. An outline of an analytic proof is as follows:

  • Rotate the loci in Fig.C.49 counterclockwise by 90 degrees.

  • Prove that the two root loci are continuous, single-valued functions of $ \nu$ (as the figure suggests).

  • Prove that for $ \left\vert\nu\right\vert>3$, there are infinitely many extrema of the log-sinc function (imaginary-part root-locus) which have negative curvature and which lie below $ \tau<-1$ (as the figure suggests). The $ \tau=-1$ and $ \nu=\pm 3$ lines are shown in the figure as dotted lines.

  • Prove that the other root locus (for the real part) has infinitely many similar extrema which occur for $ \tau>-1$ (again as the figure suggests).

  • Prove that the two root-loci interlace at $ \tau=-\infty$ (already done above).

  • Then topologically, the continuous functions must cross at infinitely many points in order to achieve interlacing at $ \tau=-\infty$.

The peaks of the log-sinc function not only indicate approximately where the left-half-plane roots occur

Reflectance Magnitude

We have shown that the conical cap reflectance has no poles in the strict right-half plane. For passivity, we also need to show that its magnitude is bounded by unity for all $ s$ on the $ j\omega $ axis.

We have

$\displaystyle R_J(j\omega) = \frac{1 - e^{-2j\omega} - 2j\omega e^{-2j\omega}}{...
...mega} - 1 + 2j\omega}
= e^{-2j\omega} \frac{\overline{D(j\omega)}}{D(j\omega)}

so that $ \left\vert R_J(j\omega)\right\vert = 1$, which is exactly lossless, as expected.

Poles at $ s=0$

We know from the above that the denominator of the cone reflectance has at least one root at $ s=0$. In this subsection we investigate this ``dc behavior'' of the cone more thoroughly.

A hasty analysis based on the reflection and transmission filters in Equations (C.154) and (C.155) might conclude that the reflectance of the conical cap converges to $ -1$ at dc, since $ R(0)=-1$ and $ T(0)=0$. However, this would be incorrect. Instead, it is necessary to take the limit as $ \omega\to0$ of the complete conical cap reflectance $ R_J(s)$:

$\displaystyle R_J(s) = \frac{1 - e^{-2s} - 2s e^{-2s}}{2s - 1 + e^{-2s}}$ (C.165)

We already discovered a root at $ s=0$ in the denominator in the context of the preceding stability proof. However, note that the numerator also goes to zero at $ s=0$. This indicates a pole-zero cancellation at dc. To find the reflectance at dc, we may use L'Hospital's rule to obtain

$\displaystyle R_J(0) = \lim_{s\to0} \frac{N^\prime(s)}{D^\prime(s)} = \lim_{s\to 0}\frac{4s e^{-2s}}{2-2e^{-2s}}$ (C.166)

and once again the limit is an indeterminate $ 0/0$ form. We therefore apply L'Hospital's rule again to obtain

$\displaystyle R_J(0) = \lim_{s\to0} \frac{N^{\prime\prime}(s)}{D^{\prime\prime}(s)} = \lim_{s\to 0}\frac{(4-8s) e^{-2s}}{4e^{-2s}} = +1$ (C.167)

Thus, two poles and zeros cancel at dc, and the dc reflectance is $ +1$, not $ -1$ as an analysis based only on the scattering filters would indicate. From a physical point of view, it makes more sense that the cone should ``look like'' a simple rigid termination of the cylinder at dc, since its length becomes small compared with the wavelength in the limit.

Another method of showing this result is to form a Taylor series expansion of the numerator and denominator:

$\displaystyle N(s)$ $\displaystyle =$ $\displaystyle 2 {s^2} - {{8 {s^3}}\over 3} + 2 {s^4} + \cdots$ (C.168)
$\displaystyle D(s)$ $\displaystyle =$ $\displaystyle 2 {s^2} - {{4 {s^3}}\over 3} + {{2 {s^4}}\over 3} + \cdots$ (C.169)

Both series begin with the term $ 2s^2$ which means both the numerator and denominator have two roots at $ s=0$. Hence, again the conclusion is two pole-zero cancellations at dc. The series for the conical cap reflectance can be shown to be

$\displaystyle R_J(s) = 1 - {{2 s}\over 3} + {{2 {s^2}}\over 9} - {{4 {s^3}}\over {135}} - {{2 {s^4}}\over {405}} + \cdots$ (C.170)

which approaches $ +1$ as $ s\to0$.

An alternative analysis of this issue is given by Benade in [37].

Next Section:
Finite-Difference Schemes
Previous Section:
The Digital Waveguide Oscillator