Free Books Physical Audio Signal Processing

Digital Waveguide Theory

In this appendix, the basic principles of digital waveguide acoustic modeling are derived from a mathematical point of view. For this, the reader is expected to have some background in linear systems and elementary physics. In particular, facility with Laplace transforms [284], Newtonian mechanics [180], and basic differential equations is assumed.

We begin with the partial differential equation (PDE) describing the ideal vibrating string, which we first digitize by converting partial derivatives to finite differences. This yields a discrete-time recursion which approximately simulates the ideal string. Next, we go back and solve the original PDE, obtaining continuous traveling waves as the solution. These traveling waves are then digitized by ordinary sampling, resulting in the digital waveguide model for the ideal string. The digital waveguide simulation is then shown to be equivalent to a particular finite-difference recursion. (This only happens for the lossless ideal vibrating string with a particular choice of sampling intervals, so it is an interesting case.) Next digital waveguides simulating lossy and dispersive vibrating strings are derived, and alternative choices of wave variables (displacement, velocity, slope, force, power, etc.) are derived. Finally, an introduction to scattering theory for digital waveguides is presented.

The Ideal Vibrating String

**Figure C.1:** The ideal vibrating string.
$\includegraphics[width=\twidth]{eps/Fphysicalstring}$

The wave equation for the ideal (lossless, linear, flexible) vibrating string, depicted in Fig.C.1, is given by

$\displaystyle Ky''= \epsilon {\ddot y}$

(C.1)

where

$\begin{displaymath}\begin{array}{rclrcl} K& \isdef & \mbox{string tension} & \qq... ...isdef & \frac{\partial}{\partial x}y(t,x) \nonumber \end{array}\end{displaymath}$

where `` $\isdef$ '' means ``is defined as.'' The wave equation is derived in §B.6. It can be interpreted as a statement of Newton's second law, ``force = mass $\times$ acceleration,'' on a microscopic scale. Since we are concerned with transverse vibrations on the string, the relevant restoring force (per unit length) is given by the string tension times the curvature of the string (

); the restoring force is balanced at all times by the inertial force per unit length of the string which is equal to mass density times transverse acceleration ( $\epsilon {\ddot y}$ ).

The Finite Difference Approximation

In the musical acoustics literature, the normal method for creating a computational model from a differential equation is to apply the so-called finite difference approximation (FDA) in which differentiation is replaced by a finite difference (see Appendix D) [481,311]. For example

$\displaystyle {\dot y}(t,x)\approx \frac{y(t,x)-y(t-T,x)}{T} \protect$

(C.2)

and

$\displaystyle y'(t,x)\approx \frac{y(t,x)-y(t,x-X)}{X} \protect$

(C.3)

where

is the time sampling interval to be used in the simulation, and

is a spatial sampling interval. These approximations can be seen as arising directly from the definitions of the partial derivatives with respect to

and

. The approximations become exact in the limit as

and

approach zero. To avoid a delay error, the second-order finite-differences are defined with a compensating time shift:

$\displaystyle {\ddot y}(t,x) \approx \frac{y(t+T,x) - 2 y(t,x) + y(t-T,x) }{T^2} \protect$

(C.4)

$\displaystyle y''(t,x) \approx \frac{y(t,x+X) - 2 y(t,x) + y(t,x-X) }{X^2} \protect$

(C.5)

The odd-order derivative approximations suffer a half-sample delay error while all even order cases can be compensated as above.

FDA of the Ideal String

Substituting the FDA into the wave equation gives

$\displaystyle K\frac{y(t,x+X) - 2 y(t,x) + y(t,x-X)}{X^2} = \epsilon \frac{y(t+T,x) - 2 y(t,x) + y(t-T,x)}{T^2}$

which can be solved to yield the following recursion for the string displacement:

$\displaystyle y(t+T,x)$	$\displaystyle =$	$\displaystyle \frac{KT^2}{\epsilon X^2} \left[ y(t,x+X) - 2 y(t,x) + y(t,x-X)\right]$
		$\displaystyle \qquad\qquad\qquad\qquad + 2 y(t,x) - y(t-T,x) \protect$

In a practical implementation, it is common to set $T=1,\, X=(\sqrt{K/\epsilon })T$ , and evaluate on the integers

and

to obtain the difference equation

$\displaystyle y(n+1,m) = y(n,m+1) + y(n,m-1) - y(n-1,m). \protect$

(C.6)

Thus, to update the sampled string displacement, past values are needed for each point along the string at time instants

and

. Then the above recursion can be carried out for time

by iterating over all

along the string.

Perhaps surprisingly, it is shown in Appendix E that the above recursion is exact at the sample points in spite of the apparent crudeness of the finite difference approximation [442]. The FDA approach to numerical simulation was used by Pierre Ruiz in his work on vibrating strings [392], and it is still in use today [74,75].

When more terms are added to the wave equation, corresponding to complex losses and dispersion characteristics, more terms of the form appear in (C.6). These higher-order terms correspond to frequency-dependent losses and/or dispersion characteristics in the FDA. All linear differential equations with constant coefficients give rise to some linear, time-invariant discrete-time system via the FDA. A general subclass of the linear, time-invariant case giving rise to ``filtered waveguides'' is

$\displaystyle \sum_{k=0}^\infty \alpha_k \frac{\partial^k y(t,x)}{\partial t^k} = \sum_{l=0}^\infty \beta_l \frac{\partial^l y(t,x)}{\partial x^l},$

(C.7)

while the fully general linear, time-invariant 2D case is

$\displaystyle \sum_{k=0}^\infty \sum_{l=0}^\infty \alpha_{k,l} \frac{\partial^k... ...nfty \beta_{m,n} \frac{\partial^m\partial^n y(t,x)}{\partial x^m \partial x^n}.$

(C.8)

A nonlinear example is

$\displaystyle \frac{\partial y(t,x)}{\partial t} = \left(\frac{\partial y(t,x)}{\partial x}\right)^2,$

(C.9)

and a time-varying example can be given by

$\displaystyle \frac{\partial y(t,x)}{\partial t} = t^2\frac{\partial y(t,x)}{\partial x}.$

(C.10)

Traveling-Wave Solution

It is easily shown that the lossless 1D wave equation $Ky''=\epsilon {\ddot y}$ is solved by any string shape which travels to the left or right with speed $c \isdeftext \sqrt{K/\epsilon }$ . Denote right-going traveling waves in general by and left-going traveling waves by , where and are assumed twice-differentiable.^C.1Then a general class of solutions to the lossless, one-dimensional, second-order wave equation can be expressed as

$\displaystyle y(t,x) = y_r\left(t-\frac{x}{c}\right) + y_l\left(t+\frac{x}{c}\right). \protect$

(C.11)

The next section derives the result that ${\ddot y}_r= c^2y''_r$ and ${\ddot y}_l= c^2y''_l$ , establishing that the wave equation is satisfied for all traveling wave shapes

and

. However, remember that the derivation of the wave equation in §B.6 assumes the string slope is much less than

at all times and positions. Finally, we show in §C.3.6 that the traveling-wave picture is general; that is, any physical state of the string can be converted to a pair of equivalent traveling force- or velocity-wave components.

An important point to note about the traveling-wave solution of the 1D wave equation is that a function of two variables has been replaced by two functions of a single variable in time units. This leads to great reductions in computational complexity.

The traveling-wave solution of the wave equation was first published by d'Alembert in 1747 [100]. See Appendix A for more on the history of the wave equation and related topics.

Traveling-Wave Partial Derivatives

Because we have defined our traveling-wave components and as having arguments in units of time, the partial derivatives with respect to time are identical to simple derivatives of these functions. Let ${\dot y}_r$ and ${\dot y}_l$ denote the (partial) derivatives with respect to time of and , respectively. In contrast, the partial derivatives with respect to are

$\begin{eqnarray*} \frac{\partial}{\partial x} y_r\left(t-\frac{x}{c}\right) &=&... ...c}\right) &=& \frac{1}{c}{\dot y}_l\left(t+ \frac{x}{c}\right). \end{eqnarray*}$

Denoting the spatial partial derivatives by and , respectively, we can write more succinctly

$\begin{eqnarray*} y'_r&=& -\frac{1}{c}{\dot y}_r\\ [5pt] y'_l&=& \frac{1}{c}{\dot y}_l, \end{eqnarray*}$

where this argument-free notation assumes the same and for all terms in each equation, and the subscript or determines whether the omitted argument is or .

Now we can see that the second partial derivatives in are

$\begin{eqnarray*} y''_r&=& \left(-\frac{1}{c}\right)^2 {\ddot y}_r= \frac{1}{c^2... ...eft(\frac{1}{c}\right)^2 {\ddot y}_l= \frac{1}{c^2} {\ddot y}_l. \end{eqnarray*}$

These relations, together with the fact that partial differention is a linear operator, establish that

$\displaystyle y(t,x) = y_r\left(t-\frac{x}{c}\right) + y_l\left(t+\frac{x}{c}\right). \protect$

obeys the ideal wave equation ${\ddot y}= c^2y''$ for all twice-differentiable functions

and

Use of the Chain Rule

These traveling-wave partial-derivative relations may be derived a bit more formally by means of the chain rule from calculus, which states that, for the composition of functions and , i.e.,

$\displaystyle y(x) = f(g(x)),$

the derivative of the composition with respect to

can be expressed according to the chain rule as

$\displaystyle y'(x) = f^\prime(g(x))g^\prime(x),$

where $f^\prime(x)$ denotes the derivative of

with respect to

To apply the chain rule to the spatial differentiation of traveling waves, define

$\begin{eqnarray*} g_r(t,x) &=& t - \frac{x}{c}\\ [10pt] g_l(t,x) &=& t + \frac{x}{c}. \end{eqnarray*}$

Then the traveling-wave components can be written as and , and their partial derivatives with respect to become

$\begin{eqnarray*} y'_r\;\isdef \; \frac{\partial}{\partial x} y_r\left[g_r(t,x)\... ...t \left(-\frac{1}{c}\right) \;\isdef \; -\frac{1}{c}{\dot y}_r, \end{eqnarray*}$

and similarly for .

String Slope from Velocity Waves

Let's use the above result to derive the slope of the ideal vibrating string From Eq. $\,$ (C.11), we have the string displacement given by

$\displaystyle y(t,x) = y_r(t-x/c) + y_l(t+x/c). \protect$

By linearity of differentiation, the string slope is given by

$\displaystyle s(t,x) \isdef \frac{\partial}{\partial x} y(t,x) = \frac{\partial}{\partial x}y_r(t-x/c) + \frac{\partial}{\partial x}y_l(t+x/c).$

Consider only the right-going component, and define

$\displaystyle {\dot y}_r(\tau) \isdef \frac{d}{d\tau} y_r(\tau)$

with $\tau\isdef t-x/c$ . By the chain rule,

$\displaystyle \frac{\partial}{\partial x}y_r(t-x/c) = \frac{dy_r}{d\tau} \cdot \frac{d\tau}{dx} = {\dot y}_r(t-x/c)\cdot\left(-\frac{1}{c}\right).$

The left-going component is similar, but with

. Thus, the string slope in terms of traveling velocity-wave components can be written as

$\displaystyle s(t,x) = \frac{1}{c}{\dot y}_l(t+x/c) - \frac{1}{c}{\dot y}_r(t-x/c).$

Wave Velocity

Because $e^{st}$ is an eigenfunction under differentiation (i.e., the exponential function is its own derivative), it is often profitable to replace it with a generalized exponential function, with maximum degrees of freedom in its parametrization, to see if parameters can be found to fulfill the constraints imposed by differential equations.

In the case of the one-dimensional ideal wave equation (Eq. $\,$ (C.1)), with no boundary conditions, an appropriate choice of eigensolution is

$\displaystyle y(t,x) = e^{st+vx}$

(C.12)

Substituting into the wave equation yields

$\begin{displaymath} \begin{array}{rclcrcl} {\dot y}& \,\mathrel{\mathop=}\,& sy... ...\quad & y''& \,\mathrel{\mathop=}\,& v^2y \nonumber \end{array}\end{displaymath}$

Defining the wave velocity (or phase velocity^C.2) as $c \isdeftext {s/v}$ , the wave equation becomes

$\displaystyle Kv^2y$	$\displaystyle =$	$\displaystyle \epsilon s^2y$	(C.13)
$\displaystyle \,\,\Rightarrow\,\,\frac{K}{\epsilon }$	$\displaystyle =$	$\displaystyle \frac{s^2}{v^2} \isdef c^2$
$\displaystyle \,\,\Rightarrow\,\,v$	$\displaystyle =$	$\displaystyle \pm \frac{s}{c}.$

Thus

$\displaystyle y(t,x) = e^{s(t\pm x/c)}$

is a solution for all

. By superposition,

$\displaystyle y(t,x) = \sum\limits_i^{} A^{+}(s_i) e^{s_i(t-x/c)}+ A^{-}(s_i) e^{s_i(t+x/c)}$

is also a solution, where $A^{+}(s_i)$ and $A^{-}(s_i)$ are arbitrary complex-valued functions of arbitrary points

in the complex plane.

D'Alembert Derived

Setting $s\isdeftext j \omega$ , and extending the summation to an integral, we have, by Fourier's theorem,

$\displaystyle y(t,x) = y_r\left(t-\frac{x}{c}\right) + y_l\left(t+\frac{x}{c}\right)$

(C.14)

for arbitrary continuous functions $y_r(\cdot)$ and $y_l(\cdot)$ . This is again the traveling-wave solution of the wave equation attributed to d'Alembert, but now derived from the eigen-property of sinusoids and Fourier theory rather than ``guessed''.

An example of the appearance of the traveling wave components shortly after plucking an infinitely long string at three points is shown in Fig.C.2.

**Figure C.2:** An infinitely long string, ``plucked'' simultaneously at *three* points, labeled ``p'' in the figure, so as to produce an initial triangular displacement. The initial displacement is modeled as the sum of two identical triangular pulses which are exactly on top of each other at time 0. At time shortly after time 0, the traveling waves centers are separated by meters, and their sum gives the trapezoidal physical string displacement at time which is also shown. Note that only three short string segments are in motion at that time: the flat top segment which is heading to zero where it will halt forever, and two short pieces on the left and right which are the leading edges of the left- and right-going traveling waves. The string is not moving where the traveling waves overlap at the same slope. When the traveling waves fully separate, the string will be at rest everywhere but for two half-amplitude triangular pulses heading off to plus and minus infinity at speed .
$\includegraphics[width=\twidth]{eps/f_t_waves_no_term}$

Converting Any String State to Traveling Slope-Wave Components

We verified in §C.3.1 above that traveling-wave components and in Eq. $\,$ (C.14) satisfy the ideal string wave equation ${\ddot y}= c^2y''$ . By definition, the physical string displacement is given by the sum of the traveling-wave components, or

$\displaystyle y(t,x) \eqsp y_r\left(t-\frac{x}{c}\right) + y_l\left(t+\frac{x}{c}\right). \protect$

(C.15)

Thus, given any pair of traveling waves

and

, we can compute a corresponding string displacement

. This leads to the question whether any initial string state can be converted to a pair of equivalent traveling-wave components. If so, then d'Alembert's traveling-wave solution is complete, and all solutions to the ideal string wave equation can be expressed in terms of traveling waves.

The state of an ideal string at time is classically specified by its displacement and velocity

$\displaystyle v(t,x)\isdefs {\dot y}(t,x)\isdefs \frac{\partial}{\partial t} y(t,x)$

for all

[317]. Equation (C.15) gives us

as a simple sum of the traveling-wave components, and now we need a formula for

in terms of them as well. It will be derived in §C.7.3 (see Equations (C.44-C.46)) that we can write

$\displaystyle v(t,x) \eqsp -cy_r^\prime\left(t-\frac{x}{c}\right) + cy_l^\prime\left(t+\frac{x}{c}\right).$

where

denotes the partial derivative with respect to

as usual. We have

$\displaystyle \left[\begin{array}{c} y(t,x) \\ [2pt] v(t,x) \end{array}\right] ... ...ght] \left[\begin{array}{c} y_r(t-x/c) \\ [2pt] y_l(t+x/c) \end{array}\right].$

Inverting the two-by-two differential operator matrix yields left- and right-going slope waves as a function of an arbitrary initial slope and velocity:

$\displaystyle \left[\begin{array}{c} y'^{+} \\ [2pt] y'^{-} \end{array}\right] ... ...eft[\begin{array}{c} y'-\frac{v}{c} \\ [2pt] y'+\frac{v}{c} \end{array}\right]$

Integrating both sides with respect to

, and choosing the constant of integration to give the correct constant component of

, we obtain the displacement-wave components

$\displaystyle \left[\begin{array}{c} y^{+} \\ [2pt] y^{-} \end{array}\right] \eqsp \frac{1}{2}\left[\begin{array}{c} y-w \\ [2pt] y+w \end{array}\right]$

where

$\displaystyle w(t,x) \isdefs \frac{1}{c}\int_{-\infty}^x v(t,\xi)\,d\xi.$

Notice that if the initial velocity is zero, each of the initial traveling displacement waves is simply half the initial displacement, as expected. On the other hand, if the initial displacement is zero and there is a uniform initial velocity (the whole string is moving), the initial displacement-wave components are unbounded as the string length goes to infinity. Related discussion appears in Appendix E.

It will be seen in §C.7.4 that state conversion between physical variables and traveling-wave components is simpler when force and velocity are chosen as the physical state variables (as opposed to displacement and velocity used here).

Sampled Traveling Waves

To carry the traveling-wave solution into the ``digital domain,'' it is necessary to sample the traveling-wave amplitudes at intervals of seconds, corresponding to a sampling rate $f_s \isdeftext 1/T$ samples per second. For CD-quality audio, we have kHz. The natural choice of spatial sampling interval is the distance sound propagates in one temporal sampling interval , or $X\isdeftext cT$ meters. In a lossless traveling-wave simulation, the whole wave moves left or right one spatial sample each time sample; hence, lossless simulation requires only digital delay lines. By lumping losses parsimoniously in a real acoustic model, most of the traveling-wave simulation can in fact be lossless even in a practical application.

Formally, sampling is carried out by the change of variables

$\begin{displaymath} \begin{array}{rclcl} x &\to& x_m&=& mX, \nonumber \\ t &\to& t_n&=& nT. \nonumber \end{array}\end{displaymath}$

Substituting into the traveling-wave solution of the wave equation gives

$\displaystyle y(t_n,x_m)$	$\displaystyle \,\mathrel{\mathop=}\,$	$\displaystyle y_r(t_n- x_m/c) + y_l(t_n+ x_m/c) \protect$
	$\displaystyle \,\mathrel{\mathop=}\,$	$\displaystyle y_r(nT- mX/c) + y_l(nT+ mX/c) \nonumber$
	$\displaystyle \,\mathrel{\mathop=}\,$	$\displaystyle y_r\left[(n-m)T\right]+ y_l\left[(n+m)T\right]. \protect$

Since

multiplies all arguments, let's suppress it by defining

$\displaystyle y^{+}(n) \isdefs y_r(nT), \qquad\qquad y^{-}(n) \isdefs y_l(nT). \protect$

(C.16)

This new notation also introduces a ``

'' superscript to denote a traveling-wave component propagating to the right, and a ``

'' superscript to denote propagation to the left. This notation is similar to that used for acoustic tubes [297].

Digital Waveguide Model

**Figure C.3:** Digital simulation of the ideal, lossless waveguide with observation points at and . (The symbol `` $z^{-1}$ '' denotes a one-sample delay.)
$\includegraphics[scale=0.9]{eps/fideal}$

In this section, we interpret the sampled d'Alembert traveling-wave solution of the ideal wave equation as a digital filtering framework. This is an example of what are generally known as digital waveguide models [430,431,433,437,442].

The term $y_r\left[(n-m)T\right]\isdef y^{+}(n-m)$ in Eq. $\,$ (C.16) can be thought of as the output of an -sample delay line whose input is $y^{+}(n)$ . In general, subtracting a positive number from a time argument corresponds to delaying the waveform by samples. Since $y^{+}$ is the right-going component, we draw its delay line with input $y^{+}(n)$ on the left and its output $y^{+}(n-m)$ on the right. This can be seen as the upper ``rail'' in Fig.C.3

Similarly, the term $y_l\left[(n+m)T\right]\isdeftext y^{-}(n+m)$ can be thought of as the input to an -sample delay line whose output is $y^{-}(n)$ . Adding to the time argument produces an -sample waveform advance. Since $y^{-}$ is the left-going component, it makes sense to draw the delay line with its input $y^{-}(n+m)$ on the right and its output $y^{-}(n)$ on the left. This can be seen as the lower ``rail'' in Fig.C.3.

Note that the position along the string, meters, is laid out from left to right in the diagram, giving a physical interpretation to the horizontal direction in the diagram. Finally, the left- and right-going traveling waves must be summed to produce a physical output according to the formula

$\displaystyle y(t_n,x_m) = y^{+}(n-m) + y^{-}(n+m). \protect$

(C.17)

We may compute the physical string displacement at any spatial sampling point

by simply adding the upper and lower rails together at position

along the delay-line pair. In Fig.C.3, ``transverse displacement outputs'' have been arbitrarily placed at

and

. The diagram is similar to that of well known ladder and lattice digital filter structures (§C.9.4,[297]), except for the delays along the upper rail, the absence of scattering junctions, and the direct physical interpretation. (A scattering junction implements partial reflection and partial transmission in the waveguide.) We could proceed to ladder and lattice filters by (1) introducing a perfectly reflecting (rigid or free) termination at the far right, and (2) commuting the delays rightward from the upper rail down to the lower rail [432,434]. The absence of scattering junctions is due to the fact that the string has a uniform wave impedance. In acoustic tube simulations, such as for voice [87,297] or wind instruments [195], lossless scattering junctions are used at changes in cross-sectional tube area and lossy scattering junctions are used to implement tone holes. In waveguide bowed-string synthesis (discussed in a later section), the bow itself creates an active, time-varying, and nonlinear scattering junction on the string at the bowing point.

Any ideal, one-dimensional waveguide can be simulated in this way. It is important to note that the simulation is exact at the sampling instants, to within the numerical precision of the samples themselves. To avoid aliasing associated with sampling, we require all waveshapes traveling along the string to be initially bandlimited to less than half the sampling frequency. In other words, the highest frequencies present in the signals and may not exceed half the temporal sampling frequency $f_s \isdeftext 1/T$ ; equivalently, the highest spatial frequencies in the shapes and may not exceed half the spatial sampling frequency $\nu_s \isdeftext 1/X$ .

A C program implementing a plucked/struck string model in the form of Fig.C.3 is available at http://ccrma.stanford.edu/~jos/pmudw/.

Digital Waveguide Interpolation

A more compact simulation diagram which stands for either sampled or continuous waveguide simulation is shown in Fig.C.4. The figure emphasizes that the ideal, lossless waveguide is simulated by a bidirectional delay line, and that bandlimited spatial interpolation may be used to construct a displacement output for an arbitrary not a multiple of , as suggested by the output drawn in Fig.C.4. Similarly, bandlimited interpolation across time serves to evaluate the waveform at an arbitrary time not an integer multiple of (§4.4).

**Figure C.4:** Simplified picture of ideal waveguide simulation.
$\includegraphics[scale=0.9]{eps/fcompact}$

Ideally, bandlimited interpolation is carried out by convolving a continuous ``sinc function'' sinc $(x)\isdeftext \sin(\pi x)/\pi x$ with the signal samples. Specifically, convolving a sampled signal with sinc ``evaluates'' the signal at an arbitrary continuous time . The sinc function is the impulse response of the ideal lowpass filter which cuts off at half the sampling rate.

In practice, the interpolating sinc function must be windowed to a finite duration. This means the associated lowpass filter must be granted a ``transition band'' in which its frequency response is allowed to ``roll off'' to zero at half the sampling rate. The interpolation quality in the ``pass band'' can always be made perfect to within the resolution of human hearing by choosing a sufficiently large product of window-length times transition-bandwidth. Given ``audibly perfect'' quality in the pass band, increasing the transition bandwidth reduces the computational expense of the interpolation. In fact, they are approximately inversely proportional. This is one reason why oversampling at rates higher than twice the highest audio frequency is helpful. For example, at a kHz sampling rate, the transition bandwidth above the nominal audio upper limit of kHz is only kHz, while at a kHz sampling rate (used in DAT machines) the guard band is kHz wide--nearly double. Since the required window length (impulse response duration) varies inversely with the provided transition bandwidth, we see that increasing the sampling rate by less than ten percent reduces the filter expense by almost fifty percent. Windowed-sinc interpolation is described further in §4.4. Many more techniques for digital resampling and delay-line interpolation are reviewed in [267].

Relation to the Finite Difference Recursion

In this section we will show that the digital waveguide simulation technique is equivalent to the recursion produced by the finite difference approximation (FDA) applied to the wave equation [442, pp. 430-431]. A more detailed derivation, with examples and exploration of implications, appears in Appendix E. Recall from (C.6) that the time update recursion for the ideal string digitized via the FDA is given by

$\displaystyle y(n+1,m) = y(n,m+1) + y(n,m-1) - y(n-1,m). \protect$

(C.18)

To compare this with the waveguide description, we substitute the traveling-wave decomposition $y(n,m) = y^{+}(n-m) + y^{-}(n+m)$ (which is exact in the ideal case at the sampling instants) into the right-hand side of the FDA recursion above and see how good is the approximation to the left-hand side $y(n+1,m) = y^{+}(n+1-m) + y^{-}(n+1+m)$ . Doing this gives

$\displaystyle y(n+1,m)$	$\displaystyle =$	$\displaystyle y(n,m+1) + y(n,m-1) - y(n-1,m)$	(C.19)
	$\displaystyle =$	$\displaystyle y^{+}(n-m-1) + y^{-}(n+m+1)$
		$\displaystyle + y^{+}(n-m+1) + y^{-}(n+m-1)$
		$\displaystyle - y^{+}(n-m-1) - y^{-}(n+m-1)$
	$\displaystyle =$	$\displaystyle y^{-}(n+m+1) + y^{+}(n-m+1)$
	$\displaystyle =$	$\displaystyle y^{+}[(n+1)-m] + y^{-}[(n+1)+m]$
	$\displaystyle \isdef$	$\displaystyle y(n+1,m).$

Thus, we obtain the result that the FDA recursion is also exact in the lossless case, because it is equivalent to the digital waveguide method which we know is exact at the sampling points. This is surprising since the FDA introduces artificial damping when applied to lumped, mass-spring systems, as discussed earlier.

The last identity above can be rewritten as

$\displaystyle y(n+1,m)$	$\displaystyle \isdef$	$\displaystyle y^{+}[(n+1)-m] + y^{-}[(n+1)+m]$	(C.20)
	$\displaystyle =$	$\displaystyle y^{+}[n-(m-1)] + y^{-}[n+(m+1)]$

which says the displacement at time

, position

, is the superposition of the right-going and left-going traveling wave components at positions

and

, respectively, from time

. In other words, the physical wave variable can be computed for the next time step as the sum of incoming traveling wave components from the left and right. This picture also underscores the lossless nature of the computation.

This results extends readily to the digital waveguide mesh (§C.14), which is essentially a lattice-work of digital waveguides for simulating membranes and volumes. The equivalence is important in higher dimensions because the finite-difference model requires less computations per node than the digital waveguide approach.

Even in one dimension, the digital waveguide and finite-difference methods have unique advantages in particular situations, and as a result they are often combined together to form a hybrid traveling-wave/physical-variable simulation [351,352,222,124,123,224,263,223]. In this hybrid simulations, the traveling-wave variables are called ``W variables'' (where `W' stands for ``Wave''), while the physical variables are caled ``K variables'' (where `K' stands for ``Kirchoff''). Each K variable, such as displacement on a vibrating string, can be regarded as the sum of two traveling-wave components, or W variables:

$\displaystyle y(nT,mX) = y_r(nT-mX/c) + y_l(nT+mX/c)$

Conversion between K variables and W variables can be non-trivial due to the non-local dependence of one set of state variables on the other, in general. A detailed examination of this issue is given in Appendix E.

A Lossy 1D Wave Equation

In any real vibrating string, there are energy losses due to yielding terminations, drag by the surrounding air, and internal friction within the string. While losses in solids generally vary in a complicated way with frequency, they can usually be well approximated by a small number of odd-order terms added to the wave equation. In the simplest case, force is directly proportional to transverse string velocity, independent of frequency. If this proportionality constant is $\mu$ , we obtain the modified wave equation

$\displaystyle Ky''= \epsilon {\ddot y}+ \mu{\dot y}.$

(C.21)

Thus, the wave equation has been extended by a ``first-order'' term, i.e., a term proportional to the first derivative of

with respect to time. More realistic loss approximations would append terms proportional to ${\dddot y}$ , ${\partial^5 y/\partial t^5}$ , and so on, giving frequency-dependent losses.

Setting $y(t,x) = e^{st+vx}$ in the wave equation to find the relationship between temporal and spatial frequencies in the eigensolution, the wave equation becomes

$\displaystyle K\left(v^2 y\right)$	$\displaystyle =$	$\displaystyle \epsilon \left(s^2 y\right)+ \mu\left(s y\right) \,\,\Rightarrow\,\,Kv^2 = \epsilon s^2 + \mu s$
$\displaystyle \,\,\Rightarrow\,\,v^2$	$\displaystyle =$	$\displaystyle \frac{\epsilon }{K} s^2 + \frac{\mu}{K} s = \frac{\epsilon }{K} s... ...on s}} \right) \isdef \frac{s^2}{c^2}\left(1 + {\frac{\mu}{\epsilon s}} \right)$
$\displaystyle \,\,\Rightarrow\,\,v$	$\displaystyle =$	$\displaystyle \pm \frac{s}{c} \sqrt{1 + {\frac{\mu}{\epsilon s}}}$

where $c \isdef \sqrt{K/\epsilon }$ is the wave velocity in the lossless case. At high frequencies (large $\vert s\vert$ ), or when the friction coefficient $\mu$ is small relative to the mass density $\epsilon$ at the lowest frequency of interest, we have the approximation

$\displaystyle \left(1 + {\frac{\mu}{\epsilon s}}\right)^\frac{1}{2} \approx 1 + \frac{1}{2}{\frac{\mu}{\epsilon s}}$

(C.22)

by the binomial theorem. For this small-loss approximation, we obtain the following relationship between temporal and spatial frequency:

$\displaystyle v \approx \pm \frac{1}{c}\left({s + \frac{\mu}{2\epsilon }} \right)$

(C.23)

The eigensolution is then

$\displaystyle e^{st+vx} = \exp{\left[st\pm \left({s + \frac{\mu}{2\epsilon }}\r... ...c{x}{c}\right)\right]} \exp{\left(\pm\frac{\mu}{2\epsilon }\frac{x}{c}\right)}.$

(C.24)

The right-going part of the eigensolution is

$\displaystyle e^{-{\left(\mu/2\epsilon \right)}{x/c}} e^{s \left(t - {x/c}\right)}$

(C.25)

which decays exponentially in the direction of propagation, and the left-going solution is

$\displaystyle e^{{\left(\mu/2\epsilon \right)}{x/c}} e^{s \left(t + {x/c}\right)}$

(C.26)

which also decays exponentially in the direction of travel.

Setting $s=j\omega$ and using superposition to build up arbitrary traveling wave shapes, we obtain the general class of solutions

$\displaystyle y(t,x) = e^{-{\left(\mu/2\epsilon \right)}{x/c}} y_r\left(t-{x/c}\right) + e^{{\left(\mu/2\epsilon \right)}{x/c}} y_l\left(t+{x/c}\right).$

(C.27)

Sampling these exponentially decaying traveling waves at intervals of seconds (or meters) gives

$\begin{eqnarray*} y(t_n,x_m) &=& e^{-{\left(\mu/2\epsilon \right)}{x_m/c}} y_r\l... ... y^{-}(n+m) \\ &\isdef & g^{m} y^{+}(n-m) + g^{-m} y^{-}(n+m). \end{eqnarray*}$

The simulation diagram for the lossy digital waveguide is shown in Fig.C.5.

**Figure C.5:** Discrete simulation of the ideal, lossy waveguide.
$\includegraphics[scale=0.9]{eps/floss}$

Again the discrete-time simulation of the decaying traveling-wave solution is an exact implementation of the continuous-time solution at the sampling positions and instants, even though losses are admitted in the wave equation. Note also that the losses which are distributed in the continuous solution have been consolidated, or lumped, at discrete intervals of meters in the simulation. The loss factor $g = e^{-{\mu T/2\epsilon }}$ summarizes the distributed loss incurred in one sampling interval. The lumping of distributed losses does not introduce an approximation error at the sampling points. Furthermore, bandlimited interpolation can yield arbitrarily accurate reconstruction between samples. The only restriction is again that all initial conditions and excitations be bandlimited to below half the sampling rate.

Loss Consolidation

In many applications, it is possible to realize vast computational savings in digital waveguide models by commuting losses out of unobserved and undriven sections of the medium and consolidating them at a minimum number of points. Because the digital simulation is linear and time invariant (given constant medium parameters $K,\epsilon ,\mu$ ), and because linear, time-invariant elements commute, the diagram in Fig.C.6 is exactly equivalent (to within numerical precision) to the previous diagram in Fig.C.5.

**Figure C.6:** Discrete simulation of the ideal, lossy waveguide. Each per-sample loss factor may be ``pushed through'' delay elements and combined with other loss factors until an input or output is encountered which inhibits further migration. If further consolidation is possible on the other side of a branching node, a loss factor can be pushed *through* the node by pushing a copy into each departing branch. If there are other *inputs* to the node, the *inverse* of the loss factor must appear on each of them. Similar remarks apply to pushing backwards through a node.
$\includegraphics[scale=0.9]{eps/flloss}$

Frequency-Dependent Losses

In nearly all natural wave phenomena, losses increase with frequency. Distributed losses due to air drag and internal bulk losses in the string tend to increase monotonically with frequency. Similarly, air absorption increases with frequency, adding loss for sound waves in acoustic tubes or open air [318].

Perhaps the apparently simplest modification to Eq. $\,$ (C.21) yielding frequency-dependent damping is to add a third-order time-derivative term [392]:

$\displaystyle Ky''= \epsilon {\ddot y}+ \mu{\dot y}+ \mu_3{\dddot y} \protect$

(C.28)

While this model has been successful in practice [77], it turns out to go unstable at very high sampling rates. The technical term for this problem is that the PDE is ill posed [45].

A well posed replacement for Eq. $\,$ (C.28) is given by

$\displaystyle Ky''= \epsilon {\ddot y}+ \mu{\dot y}+ \mu_2{\dot y''} \protect$

(C.29)

in which the third-order partial derivative with respect to time, ${\dddot y}$ , has been replaced by a third-order mixed partial derivative--twice with respect to

and once with respect to

The solution of a lossy wave equation containing higher odd-order derivatives with respect to time yields traveling waves which propagate with frequency-dependent attenuation. Instead of scalar factors distributed throughout the diagram as in Fig.C.5, each factor becomes a lowpass filter having some frequency-response per sample denoted by $G(\omega)$ . Because propagation is passive, we will always have $\left\vert G(\omega)\right\vert\leq 1$ .

More specically, As shown in [392], odd-order partial derivatives with respect to time in the wave equation of the form

$\displaystyle \frac{\partial^n}{\partial t^n} y(t,x), \quad n=1,3,5,\ldots,$

correspond to attenuation of traveling waves on the string. (The even-order time derivatives can be associated with variations in dispersion as a function of frequency, which is considered in §C.6 below.) For

, the losses are frequency dependent, and the per-sample amplitude-response ``rolls off'' proportional to

$\displaystyle G(\omega) \propto \frac{1}{\omega^{n-1}},$ $\displaystyle \mbox{($n$\ odd)}$ $\displaystyle .$

In particular, if the wave equation (C.21) is modified by adding terms proportional to $\mu_3{\dddot y}$ and $\mu_5{\partial^5 y/\partial t^5}$ , for instance, then the per-sample propagation gain $G(\omega)$ has the form

$\displaystyle G(\omega) = g_0 + g_2 \omega^2 + g_4 \omega^4$

where the

are constants depending on the constants $\mu_3$ and $\mu_5$ in the wave equation. Since these per-sample loss filters are linear and time-invariant [449], they may also be consolidated at a minimum number of points in the waveguide without introducing any approximation error, just like the constant gains

in Fig.C.5. This result does not extend precisely to the waveguide mesh (§C.14).

In view of the above, we see that we can add odd-order time derivatives to the wave equation to approximate experimentally observed frequency-dependent damping characteristics in vibrating strings [73]. However, we then have the problem that such wave equations are ill posed, leading to possible stability failures at high sampling rates. As a result, it is generally preferable to use mixed derivatives, as in Eq. $\,$ (C.29), and try to achieve realistic damping using higher order spatial derivatives instead.

Well Posed PDEs for Modeling Damped Strings

A large class of well posed PDEs is given by [45]

$\displaystyle {\ddot y} + 2\sum_{k=0}^M q_k \frac{\partial^{2k+1}y}{\partial x^{2k}\partial t} + \sum_{k=0}^N r_k\frac{\partial^{2k}y}{\partial x^{2k}}. \protect$

(C.30)

Thus, to the ideal string wave equation Eq. $\,$ (C.1) we add any number of even-order partial derivatives in

, plus any number of mixed odd-order partial derivatives in

and

, where differentiation with respect to

occurs only once. Because every member of this class of PDEs is only second-order in time, it is guaranteed to be well posed, as shown in §D.2.2.

Digital Filter Models of Damped Strings

In an efficient digital simulation, lumped loss factors of the form $G^k(\omega)$ are approximated by a rational frequency response ${\hat G}_k(e^{j\omega T})$ . In general, the coefficients of the optimal rational loss filter are obtained by minimizing $\vert\vert\,G^k(\omega) - {\hat G}_k(e^{j\omega T})\,\vert\vert$ with respect to the filter coefficients or the poles and zeros of the filter. To avoid introducing frequency-dependent delay, the loss filter should be a zero-phase, finite-impulse-response (FIR) filter [362]. Restriction to zero phase requires the impulse response ${\hat g}_k(n)$ to be finite in length (i.e., an FIR filter) and it must be symmetric about time zero, i.e., ${\hat g}_k(-n)={\hat g}_k(n)$ . In most implementations, the zero-phase FIR filter can be converted into a causal, linear phase filter by reducing an adjacent delay line by half of the impulse-response duration.

Lossy Finite Difference Recursion

We will now derive a finite-difference model in terms of string displacement samples which correspond to the lossy digital waveguide model of Fig.C.5. This derivation generalizes the lossless case considered in §C.4.3.

Figure C.7 depicts a digital waveguide section once again in ``physical canonical form,'' as shown earlier in Fig.C.5, and introduces a doubly indexed notation for greater clarity in the derivation below [442,222,124,123].

**Figure C.7:** Lossy digital waveguide--frequency-independent loss-factors .
$\includegraphics{eps/wglossy}$

Referring to Fig.C.7, we have the following time-update relations:

$\begin{eqnarray*} y^{+}_{n+1,m}&=& gy^{+}_{n,m-1}\;=\; g\cdot(y_{n,m-1}- y^{-}_{... ..._{n+1,m}&=& gy^{-}_{n,m+1}\;=\; g\cdot(y_{n,m+1}- y^{+}_{n,m+1}) \end{eqnarray*}$

Adding these equations gives

$\displaystyle y_{n+1,m}$	$\displaystyle =$	$\displaystyle g\cdot(y_{n,m-1}+y_{n,m+1}) - g\cdot(\underbrace{y^{-}_{n,m-1}}_{gy^{-}_{n-1,m}} + \underbrace{y^{+}_{n,m+1}}_{gy^{+}_{n-1,m}})$
	$\displaystyle =$	$\displaystyle g\cdot(y_{n,m-1}+y_{n,m+1}) - g^2 y_{n-1,m} \protect$	(C.31)

This is now in the form of the finite-difference time-domain (FDTD) scheme analyzed in [222]:

$\displaystyle y_{n+1,m}= g^{+}_my_{n,m-1}+ g^{-}_my_{n,m+1}+ a_my_{n-1,m},$

with $g^{+}_m= g^{-}_m= g$ , and

. In [124], it was shown by von Neumann analysis (§D.4) that these parameter choices give rise to a stable finite-difference scheme (§D.2.3), provided $\vert g\vert\leq 1$ . In the present context, we expect stability to follow naturally from starting with a passive digital waveguide model.

Frequency-Dependent Losses

The preceding derivation generalizes immediately to frequency-dependent losses. First imagine each in Fig.C.7 to be replaced by , where for passivity we require

$\displaystyle \left\vert G(e^{j\omega T})\right\vert\leq 1.$

In the time domain, we interpret

as the impulse response corresponding to

. We may now derive the frequency-dependent counterpart of Eq. $\,$ (C.31) as follows:

$\begin{eqnarray*} y^{+}_{n+1,m}&=& g\ast y^{+}_{n,m-1}\;=\; g\ast (y_{n,m-1}- y^... ... &=& g\ast \left[(y_{n,m-1}+y_{n,m+1}) - g\ast y_{n-1,m}\right] \end{eqnarray*}$

where $\ast$ denotes convolution (in the time dimension only). Define filtered node variables by

$\begin{eqnarray*} y^f_{n,m}&=& g\ast y_{n,m}\\ y^{ff}_{n,m}&=& g\ast y^f_{n,m}. \end{eqnarray*}$

Then the frequency-dependent FDTD scheme is simply

$\displaystyle y_{n+1,m}= y^f_{n,m-1}+ y^f_{n,m+1}- y^{ff}_{n-1,m}.$

We see that generalizing the FDTD scheme to frequency-dependent losses requires a simple filtering of each node variable $y_{n,m}$ by the per-sample propagation filter

. For computational efficiency, two spatial lines should be stored in memory at time

: $y^f_{n,m}$ and $y^{ff}_{n-1,m}$ , for all

. These state variables enable computation of $y_{n+1,m}$ , after which each sample of $y^f_{n,m}$ ( $\forall m$ ) is filtered by

to produce $y^{ff}_{n-1,m}$ for the next iteration, and $y_{n+1,m}$ is filtered by

to produce $y^f_{n,m}$ for the next iteration.

The frequency-dependent generalization of the FDTD scheme described in this section extends readily to the digital waveguide mesh. See §C.14.5 for the outline of the derivation.

The Dispersive 1D Wave Equation

In the ideal vibrating string, the only restoring force for transverse displacement comes from the string tension (§C.1 above); specifically, the transverse restoring force is equal the net transverse component of the axial string tension. Consider in place of the ideal string a bundle of ideal strings, such as a stranded cable. When the cable is bent, there is now a new restoring force arising from some of the fibers being compressed and others being stretched by the bending. This force sums with that due to string tension. Thus, stiffness in a vibrating string introduces a new restoring force proportional to bending angle. It is important to note that string stiffness is a linear phenomenon resulting from the finite diameter of the string.

In typical treatments,^C.3bending stiffness adds a new term to the wave equation that is proportional to the fourth spatial derivative of string displacement:

$\displaystyle \epsilon {\ddot y}= Ky''- \kappa y'''' \protect$

(C.32)

where the moment constant $\kappa = YI$ is the product of Young's modulus

(the ``relative-displacement spring constant per unit cross-sectional area,'' discussed in §B.5.1) and the area moment of inertia

(§B.4.8); as derived in §B.4.9, a cylindrical string of radius

has area moment of inertia equal to $\pi a^2 \cdot (a/2)^2 = \pi a^4/4$ . This wave equation works well enough for small amounts of bending stiffness, but it is clearly missing some terms because it predicts that deforming the string into a parabolic shape will incur no restoring force due to stiffness. See §6.9 for further discussion of wave equations for stiff strings.

To solve the stiff wave equation Eq. $\,$ (C.32), we may set $y(t,x) = e^{st+vx}$ to get

$\displaystyle \epsilon s^2 = Kv^2 - \kappa v^4.$

At very low frequencies, or when stiffness is negligible in comparison with

, we obtain again the non-stiff string: $\epsilon s^2\approx Kv^2 \,\,\Rightarrow\,\,v=\pm s/c$ .

At very high frequencies, or when the tension is negligible relative to $\kappa v^2$ , we obtain the ideal bar (or rod) approximation:

$\displaystyle \epsilon s^2 \approx -\kappa v^4 \,\,\Rightarrow\,\,v \approx \pm e^{\pm j\frac{\pi}{4}} \left(\frac{\epsilon }{\kappa} \right)^{1/4}\sqrt{s}$

In an ideal bar, the only restoring force is due to bending stiffness. Setting $s=j\omega$ gives solutions $v=\pm j\left( \epsilon /\kappa\right)^{1/4}\sqrt{\omega}$ and $v=\pm\left( \epsilon /\kappa\right)^{1/4}\sqrt{\omega}$ . In the first case, the wave velocity becomes proportional to $\sqrt{\omega}$ . That is, waves travel faster along the ideal bar as oscillation frequency increases, going up as the square root of frequency. The second solution corresponds to a change in the wave shape which prevents sharp corners from forming due to stiffness [95,118].

At intermediate frequencies, between the ideal string and the ideal bar, the stiffness contribution can be treated as a correction term [95]. This is the region of most practical interest because it is the principal operating region for strings, such as piano strings, whose stiffness has audible consequences (an inharmonic, stretched overtone series). Assuming $\kappa_0 \isdeftext \kappa/K\ll 1$ ,

$\displaystyle s^2$	$\displaystyle =$	$\displaystyle \frac{K}{\epsilon } v^2 - \frac{\kappa}{\epsilon } v^4 = \frac{K}... ...(1 - \frac{\kappa}{K} v^2\right) \isdef c_0^2 v^2 \left(1 - \kappa_0 v^2\right)$
$\displaystyle \,\,\Rightarrow\,\,v^2$	$\displaystyle \approx$	$\displaystyle \frac{s^2}{c_0^2} \left(1+\kappa_0 v^2\right) \approx \frac{s^2}{c_0^2} \left(1+\kappa_0 \frac{s^2}{c_0^2}\right)$
$\displaystyle \,\,\Rightarrow\,\,v$	$\displaystyle \approx$	$\displaystyle \pm \frac{s}{c_0} \sqrt{1+\kappa_0 \frac{s^2}{c_0^2}} \approx \pm \frac{s}{c_0} \left(1+\frac{1}{2}\kappa_0 \frac{s^2}{c_0^2} \right).$

Substituting for

in terms of

in $e^{st+vx}$ gives the general eigensolution

$\displaystyle e^{st+vx} = \exp{\left\{{s\left[t\pm \frac{x}{c_0}\left( 1+\frac{1}{2}\kappa_0 \frac{s^2}{c_0^2} \right)\right]}\right\}}.$

Setting $s=j\omega$ as before, corresponding to driving the medium sinusoidally over time at frequency $\omega$ , the medium response is

$\displaystyle e^{st+vx} = e^{j\omega\left[t\pm {x/ c(\omega)}\right]}$

where

$\displaystyle c(\omega) \isdef c_0\left(1 + \frac{\kappa\omega^2}{2Kc_0^2}\right).$

Because the effective wave velocity depends on $\omega$ , we cannot use Fourier's theorem to construct arbitrary traveling shapes by superposition. At

, we can construct any function of time, but the waveshape disperses as it propagates away from

. The higher-frequency Fourier components travel faster than the lower-frequency components.

Since the temporal and spatial sampling intervals are related by , this must generalize to $X= c(\omega)T\,\,\Rightarrow\,\, T(\omega)=X/c(\omega)=c_0T_0/c(\omega)$ , where is the size of a unit delay in the absence of stiffness. Thus, a unit delay $z^{-1}$ may be replaced by

$\displaystyle z^{-1}\to z^{-c_0/c(\omega)}$ (for frequency-dependent wave velocity).

That is, each delay element becomes an allpass filter which approximates the required delay versus frequency. A diagram appears in Fig.C.8, where

denotes the allpass filter which provides a rational approximation to $z^{-c_0/c(\omega)}$ .

**Figure C.8:** Section of a stiff string where allpass filters play the role of unit delay elements.
$\includegraphics[scale=0.9]{eps/fstiffstring}$

The general, order , allpass filter is given by [449]

$\displaystyle H_a(z) \isdef z^{-L} {A(z^{-1})/A(z)}$

where

$\displaystyle A(z) \isdef 1 + a_1 z^{-1}+ a_2 z^{-2} + \cdots + a_L z^{-L}$

and the roots of

must all have modulus less than

. That is, the numerator polynomial is just the reverse of the denominator polynomial. This implies each pole

is gain-compensated by a zero at

For computability of the string simulation in the presence of scattering junctions, there must be at least one sample of pure delay along each uniform section of string. This means for at least one allpass filter in Fig.C.8, we must have $H_a(\infty)=0$ which implies can be factored as $z^{-1}H_a'(z)$ . In a systolic VLSI implementation, it is desirable to have at least one real delay from the input to the output of every allpass filter, in order to be able to pipeline the computation of all of the allpass filters in parallel. Computability can be arranged in practice by deciding on a minimum delay, (e.g., corresponding to the wave velocity at a maximum frequency), and using an allpass filter to provide excess delay beyond the minimum.

Because allpass filters are linear and time invariant, they commute like gain factors with other linear, time-invariant components. Fig.C.9 shows a diagram equivalent to Fig.C.8 in which the allpass filters have been commuted and consolidated at two points. For computability in all possible contexts (e.g., when looped on itself), a single sample of delay is pulled out along each rail. The remaining transfer function, in the example of Fig.C.9, can be approximated using any allpass filter design technique [1,2,267,272,551]. Alternatively, both gain and dispersion for a stretch of waveguide can be provided by a single filter which can be designed using any general-purpose filter design method which is sensitive to frequency-response phase as well as magnitude; examples include equation error methods (such as used in the matlab invfreqz function (§8.6.4), and Hankel norm methods [177,428,36].

**Figure C.9:** Section of a stiff string where the allpass delay elements are consolidated at two points, and a sample of pure delay is extracted from each allpass chain.
$\includegraphics[scale=0.9]{eps/flstiffstring}$

In the case of a lossless, stiff string, if denotes the consolidated allpass transfer function, it can be argued that the filter design technique used should minimize the phase-delay error, where phase delay is defined by [362]

$\displaystyle P_c(\omega) \isdefs - \frac{\angle H_c\left(e^{j\omega T}\right)}{\omega}.$ (Phase Delay)

Minimizing the Chebyshev norm of the phase-delay error,

$\displaystyle \vert\vert\,P_c(\omega)-c_0/c(\omega)\,\vert\vert _\infty$

approximates minimization of the error in mode tuning for the freely vibrating string [428, pp. 182-184]. Since the stretching of the overtone series is typically what we hear most in a stiff, vibrating string, the worst-case phase-delay error is a good choice in such a case.

Alternatively, a lumped allpass filter can be designed by minimizing group delay,

$\displaystyle D_c(\omega) \isdefs - \frac{ d\angle H_c\left(e^{j\omega T}\right)}{d\omega}.$ (Group Delay)

The group delay of a filter gives the delay experienced by the amplitude envelope of a narrow frequency band centered at $\omega$ , while the phase delay applies to the ``carrier'' at $\omega$ , or a sinusoidal component at frequency $\omega$ [342]. As a result, for proper tuning of overtones, phase delay is what matters, while for precisely estimating (or controlling) the decay time in a lossy waveguide, group delay gives the effective filter delay ``seen'' by the exponential decay envelope.

See §9.4.1 for designing allpass filters with a prescribed delay versus frequency. To model stiff strings, the allpass filter must supply a phase delay which decreases as frequency increases. A good approximation may require a fairly high-order filter, adding significantly to the cost of simulation. (For low-pitched piano strings, order 8 allpass filters work well perceptually [1].) To a large extent, the allpass order required for a given error tolerance increases as the number of lumped frequency-dependent delays is increased. Therefore, increased dispersion consolidation is accompanied by larger required allpass filters, unlike the case of resistive losses.

The function piano_dispersion_filter in the Faust distribution (in effect.lib) designs and implements an allpass filter modeling the dispersion due to stiffness in a piano string [154,170,368].

Higher Order Terms

The complete, linear, time-invariant generalization of the lossy, stiff string is described by the differential equation

$\displaystyle \sum_{k=0}^\infty \alpha_k \frac{\partial^k y(t,x)}{\partial t^k} = \sum_{l=0}^\infty \beta_l \frac{\partial^l y(t,x)}{\partial x^l}. \protect$

(C.33)

which, on setting $y(t,x) = e^{st+vx}$ , (or taking the 2D Laplace transform with zero initial conditions), yields the algebraic equation,

$\displaystyle \sum_{k=0}^\infty \alpha_k s^k = \sum_{l=0}^\infty \beta_l v^l$

(C.34)

Solving for

in terms of

is, of course, nontrivial in general. However, in specific cases, we can determine the appropriate attenuation per sample $G(\omega)$ and wave propagation speed $c(\omega)$ by numerical means. For example, starting at

, we normally also have

(corresponding to the absence of static deformation in the medium). Stepping

forward by a small differential $j{{\Delta}}\omega$ , the left-hand side can be approximated by $\alpha_0+\alpha_1j{{\Delta}}\omega$ . Requiring the generalized wave velocity

to be continuous, a physically reasonable assumption, the right-hand side can be approximated by $\beta_0+\beta_1 \Delta v$ , and the solution is easy. As

steps forward, higher order terms become important one by one on both sides of the equation. Each new term in

spawns a new solution for

in terms of

, since the order of the polynomial in

is incremented. It appears possible that homotopy continuation methods [316] can be used to keep track of the branching solutions of

as a function of

. For each solution

, let ${v_r}(\omega)$ denote the real part of $v(j\omega)$ and let ${v_i}(\omega)$ denote the imaginary part. Then the eigensolution family can be seen in the form $\exp{\left\{j\omega t\pm v(j\omega)x\right\}}=\exp{\left\{\pm{v_r}(\omega)x\right\}}\cdot\exp{\left\{j\omega\left(t \pm {v_i}(\omega)x/\omega\right)\right\}}$ . Defining $c(\omega)\isdeftext \omega/{v_i}(\omega)$ , and sampling according to $x\to x_m\isdeftext mX$ and $t\to t_n\isdeftext nT(\omega)$ , with $X\isdeftext c(\omega)T(\omega)$ as before, (the spatial sampling period is taken to be frequency invariant, while the temporal sampling interval is modulated versus frequency using allpass filters), the left- and right-going sampled eigensolutions become

$\displaystyle e^{j\omega t_n\pm v(j\omega)x_m}$	$\displaystyle =$	$\displaystyle e^{\pm{v_r}(\omega)x_m}\cdot e^{ j\omega\left(t_n\pm x_m/c(\omega)\right)}$	(C.35)
	$\displaystyle =$	$\displaystyle G^m(\omega)\cdot e^{ j\omega\left(n \pm m\right)T(\omega)}$

where $G(\omega)\isdef e^{\pm{v_r}(\omega)X}$ . Thus, a general map of

versus

, corresponding to a partial differential equation of any order in the form (C.33), can be translated, in principle, into an accurate, local, linear, time-invariant, discrete-time simulation. The boundary conditions and initial state determine the initial mixture of the various solution branches as usual.

We see that a large class of wave equations with constant coefficients, of any order, admits a decaying, dispersive, traveling-wave type solution. Even-order time derivatives give rise to frequency-dependent dispersion and odd-order time derivatives correspond to frequency-dependent losses. The corresponding digital simulation of an arbitrarily long (undriven and unobserved) section of medium can be simplified via commutativity to at most two pure delays and at most two linear, time-invariant filters.

Every linear, time-invariant filter can be expressed as a zero-phase filter in series with an allpass filter. The zero-phase part can be interpreted as implementing a frequency-dependent gain (damping in a digital waveguide), and the allpass part can be seen as frequency-dependent delay (dispersion in a digital waveguide).

Alternative Wave Variables

We have thus far considered discrete-time simulation of transverse displacement in the ideal string. It is equally valid to choose velocity $v\isdeftext {\dot y}$ , acceleration $a\isdeftext {\ddot y}$ , slope , or perhaps some other derivative or integral of displacement with respect to time or position. Conversion between various time derivatives can be carried out by means integrators and differentiators, as depicted in Fig.C.10. Since integration and differentiation are linear operators, and since the traveling wave arguments are in units of time, the conversion formulas relating , , and hold also for the traveling wave components $y^\pm , v^\pm , a^\pm$ .

**Figure C.10:** Conversions between various time derivatives of displacement: displacement, $v = {\dot y}=$ velocity, $a = {\ddot y}=$ acceleration, where ${\dot y}$ denotes and ${\ddot y}$ denotes .
$\includegraphics[scale=0.9]{eps/fwaveconversions}$

Differentiation and integration have a simple form in the frequency domain. Denoting the Laplace Transform of by

$\displaystyle Y(s,x) \isdefs {\cal L}_s\{y(\cdot,x)\} \isdefs \int_0^\infty y(t,x) e^{-st} dt,$

(C.36)

where `` $\cdot$ '' in the time argument means ``for all time,'' we have, according to the differentiation theorem for Laplace transforms [284],

$\displaystyle {\cal L}_s\{{\dot y}(\cdot,x)\} \eqsp s Y(s,x) - y(0,x).$

(C.37)

Similarly, ${\cal L}_s\{\dot y^{+}\} = s Y^{+}(s) - y^{+}(0)$ , and so on. Thus, in the frequency domain, the conversions between displacement, velocity, and acceleration appear as shown in Fig.C.11.

**Figure C.11:** Conversions between various time derivatives of displacement in the frequency domain.
$\includegraphics[scale=0.9]{eps/ffdwaveconversions}$

In discrete time, integration and differentiation can be accomplished using digital filters [362]. Commonly used first-order approximations are shown in Fig.C.12.

**Figure C.12:** Simple approximate conversions between time derivatives in the discrete-time case: a) The first-order difference ${\hat v}(n) = y(n) - y(n-1)$ . b) The first-order ``leaky'' integrator ${\hat y}(n) = v(n) + g {\hat y}(n-1)$ with loss factor (slightly less than ) used to avoid infinite dc build-up.
$\includegraphics[width=\twidth]{eps/fdigitaldiffint}$

If discrete-time acceleration is defined as the sampled version of continuous-time acceleration, i.e., $a_d(n) \isdeftext a(nT,x)={\ddot y}(nT,x)$ , (for some fixed continuous position which we suppress for simplicity of notation), then the frequency-domain form is given by the transform [485]:

$\displaystyle A_d(z) \isdefs \sum_{n=0}^\infty a_d(n) z^{-n}$

(C.38)

In the frequency domain for discrete-time systems, the first-order approximate conversions appear as shown in Fig.C.13.

**Figure C.13:** Frequency-domain description of the approximate conversions between time derivatives in the discrete-time case. The subscript ``'' denotes the ``digital'' case. A ``hat'' over a variable indicates it is an approximation to the variable without the hat.
$\includegraphics[scale=0.6]{eps/ffddigitaldiffint}$

The transform plays the role of the Laplace transform for discrete-time systems. Setting $z=e^{sT}$ , it can be seen as a sampled Laplace transform (divided by ), where the sampling is carried out by halting the limit of the rectangle width at in the definition of a Reimann integral for the Laplace transform. An important difference between the two is that the frequency axis in the Laplace transform is the imaginary axis (the `` $j\omega$ axis''), while the frequency axis in the plane is on the unit circle $z = e^{j\omega T}$ . As one would expect, the frequency axis for discrete-time systems has unique information only between frequencies $-\pi/T$ and $\pi/T$ while the continuous-time frequency axis extends to plus and minus infinity.

These first-order approximations are accurate (though scaled by ) at low frequencies relative to half the sampling rate, but they are not ``best'' approximations in any sense other than being most like the definitions of integration and differentiation in continuous time. Much better approximations can be obtained by approaching the problem from a digital filter design viewpoint, as discussed in §8.6.

Spatial Derivatives

In addition to time derivatives, we may apply any number of spatial derivatives to obtain yet more wave variables to choose from. The first spatial derivative of string displacement yields slope waves

$\displaystyle y'(t,x)$	$\displaystyle \isdef$	$\displaystyle \frac{\partial}{\partial x}y(t,x)$
	$\displaystyle =$	$\displaystyle y'_r(t-x/c) + y'_l(t+x/c)$
	$\displaystyle =$	$\displaystyle -\frac{1}{c} {\dot y}_r(t-x/c) + \frac{1}{c}{\dot y}_l(t+x/c), \protect$	(C.39)

or, in discrete time,

$\displaystyle y'(t_n,x_m)$	$\displaystyle \isdef$	$\displaystyle y'(nT,mX)$
	$\displaystyle =$	$\displaystyle y'_r\left[(n-m)T\right]+ y'_l\left[(n+m)T\right]$
	$\displaystyle \isdef$	$\displaystyle y'^{+}(n-m) + y'^{-}(n+m)$
	$\displaystyle =$	$\displaystyle -\frac{1}{c} \dot y^{+}(n-m) + \frac{1}{c}\dot y^{-}(n+m)$
	$\displaystyle \isdef$	$\displaystyle -\frac{1}{c} v^{+}(n-m) + \frac{1}{c}v^{-}(n+m)$
	$\displaystyle =$	$\displaystyle \frac{1}{c} \left[v^{-}(n+m) - v^{+}(n-m) \right]. \protect$	(C.40)

From this we may conclude that $v^{-}= cy'^{-}$ and $v^{+}= -cy'^{+}$ . That is, traveling slope waves can be computed from traveling velocity waves by dividing by

and negating in the right-going case. Physical string slope can thus be computed from a velocity-wave simulation in a digital waveguide by subtracting the upper rail from the lower rail and dividing by

By the wave equation, curvature waves, $y''= {\ddot y}/c^2$ , are simply a scaling of acceleration waves, in the case of ideal strings.

In the field of acoustics, the state of a vibrating string at any instant of time is often specified by the displacement and velocity ${\dot y}(t_0,x)$ for all [317]. Since velocity is the sum of the traveling velocity waves and displacement is determined by the difference of the traveling velocity waves, viz., from Eq. $\,$ (C.39),

$\displaystyle y(t,x) \eqsp \int_0^{x} y'(t,\xi)d\xi \eqsp -\frac{1}{c}\int_0^{x} \left[v_r(t-\xi/c) - v_l(t+\xi/c)\right]d\xi,$

one state description can be converted to the other.

In summary, all traveling-wave variables can be computed from any one, as long as both the left- and right-going component waves are available. Alternatively, any two linearly independent physical variables, such as displacement and velocity, can be used to compute all other wave variables. Wave variable conversions requiring differentiation or integration are relatively expensive since a large-order digital filter is necessary to do it right (§8.6.1). Slope and velocity waves can be computed from each other by simple scaling, and curvature waves are identical to acceleration waves to within a scale factor.

In the absence of factors dictating a specific choice, velocity waves are a good overall choice because (1) it is numerically easier to perform digital integration to get displacement than it is to differentiate displacement to get velocity, (2) slope waves are immediately computable from velocity waves. Slope waves are important because they are a simple scaling of force waves.

Force Waves

**Figure C.14:** Transverse force propagation in the ideal string.
$\includegraphics[width=\twidth]{eps/fStringForce}$

Referring to Fig.C.14, at an arbitrary point along the string, the vertical force applied at time to the portion of string to the left of position by the portion of string to the right of position is given by

$\displaystyle f_{rl}(t,x) = K\sin(\theta) \approx K\tan(\theta) = Ky'(t,x)$

(C.41)

assuming $\left\vert y'(t,x)\right\vert \ll 1$ , as is assumed in the derivation of the wave equation. Similarly, the force applied by the portion to the left of position

to the portion to the right is given by

$\displaystyle f_{lr}(t,x) = - K\sin(\theta) \approx - Ky'(t,x).$

(C.42)

These forces must cancel since a nonzero net force on a massless point would produce infinite acceleration. I.e., we must have $f_{rl} + f_{lr} \equiv 0$ at all times

and positions

Vertical force waves propagate along the string like any other transverse wave variable (since they are just slope waves multiplied by tension ). We may choose either $f_{rl}$ or $f_{lr}$ as the string force wave variable, one being the negative of the other. It turns out that to make the description for vibrating strings look the same as that for air columns, we have to pick $f_{lr}$ , the one that acts to the right. This makes sense intuitively when one considers longitudinal pressure waves in an acoustic tube: a compression wave traveling to the right in the tube pushes the air in front of it and thus acts to the right. We therefore define the force wave variable to be

$\displaystyle f(t,x) \isdefs f_{lr}(t,x) = -Ky'(t,x).$

(C.43)

Note that a negative slope pulls up on the segment to the right. At this point, we have not yet considered a traveling-wave decomposition.

Wave Impedance

Using the above identities, we have that the force distribution along the string is given in terms of velocity waves by

$\displaystyle f(t,x) = \frac{K}{c} \left[{\dot y}_r(t-x/c) - {\dot y}_l(t+x/c) \right], \protect$

(C.44)

where $K/c \isdef K/\sqrt{K/\epsilon } = \sqrt{K\epsilon }$ . This is a fundamental quantity known as the wave impedance of the string (also called the characteristic impedance), denoted as

$\displaystyle R\isdefs \sqrt{K\epsilon } \eqsp \frac{K}{c} \eqsp \epsilon c.$

(C.45)

The wave impedance can be seen as the geometric mean of the two resistances to displacement: tension (spring force) and mass (inertial force).

The digitized traveling force-wave components become

$\begin{displaymath}\begin{array}{rcrl} f^{{+}}(n)&=&&R\,v^{+}(n) \\ f^{{-}}(n)&=&-&R\,v^{-}(n) \end{array} \protect\end{displaymath}$

(C.46)

which gives us that the right-going force wave equals the wave impedance times the right-going velocity wave, and the left-going force wave equals minus the wave impedance times the left-going velocity wave.^C.4Thus, in a traveling wave, force is always in phase with velocity (considering the minus sign in the left-going case to be associated with the direction of travel rather than a

degrees phase shift between force and velocity). Note also that if the left-going force wave were defined as the string force acting to the left, the minus sign would disappear. The fundamental relation $f^{{+}}= Rv^{+}$ is sometimes referred to as the mechanical counterpart of Ohm's law for traveling waves, and

in c.g.s. units can be called acoustical ohms [261].

In the case of the acoustic tube [317,297], we have the analogous relations

$\begin{displaymath}\begin{array}{rcrl} p^+(n) &=& &R_{\hbox{\tiny T}}\, u^{+}(n)... ...p^-(n) &=& -&R_{\hbox{\tiny T}}\, u^{-}(n) \end{array} \protect\end{displaymath}$

(C.47)

where

is the right-going traveling longitudinal pressure wave component,

is the left-going pressure wave, and $u^\pm (n)$ are the left and right-going volume velocity waves. In the acoustic tube context, the wave impedance is given by

$\displaystyle R_{\hbox{\tiny T}}= \frac{\rho c}{A}$ (Acoustic Tubes)

(C.48)

where $\rho$ is the mass per unit volume of air,

is sound speed in air, and

is the cross-sectional area of the tube. Note that if we had chosen particle velocity rather than volume velocity, the wave impedance would be $R_0=\rho c$ instead, the wave impedance in open air. Particle velocity is appropriate in open air, while volume velocity is the conserved quantity in acoustic tubes or ``ducts'' of varying cross-sectional area.

State Conversions

In §C.3.6, an arbitrary string state was converted to traveling displacement-wave components to show that the traveling-wave representation is complete, i.e., that any physical string state can be expressed as a pair of traveling-wave components. In this section, we revisit this topic using force and velocity waves.

By definition of the traveling-wave decomposition, we have

$\begin{eqnarray*} f&=&f^{{+}}+f^{{-}}\\ v&=&v^{+}+v^{-}. \end{eqnarray*}$

Using Eq. $\,$ (C.46), we can eliminate $v^{+}=f^{{+}}/R$ and $v^{-}=-f^{{+}}/R$ , giving, in matrix form,

$\displaystyle \left[\begin{array}{c} f \\ [2pt] v \end{array}\right] = \left[\b... ...ay}\right] \left[\begin{array}{c} f^{{+}} \\ [2pt] f^{{-}} \end{array}\right].$

Thus, the string state (in terms of force and velocity) is expressed as a linear transformation of the traveling force-wave components. Using the Ohm's law relations to eliminate instead $f^{{+}}= Rv^{+}$ and $f^{{-}}=-Rv^{-}$ , we obtain

$\displaystyle \left[\begin{array}{c} f \\ [2pt] v \end{array}\right] = \left[\b... ...d{array}\right]\left[\begin{array}{c} v^{+} \\ [2pt] v^{-} \end{array}\right].$

To convert an arbitrary initial string state

to either a traveling force-wave or velocity-wave simulation, we simply must be able to invert the appropriate two-by-two matrix above. That is, the matrix must be nonsingular. Requiring both determinants to be nonzero yields the condition

$\displaystyle 0 < R < \infty.$

That is, the wave impedance must be a positive, finite number. This restriction makes good physical sense because one cannot propagate a finite-energy wave in either a zero or infinite wave impedance.

Carrying out the inversion to obtain force waves $(f^{{+}},f^{{-}})$ from yields

$\displaystyle \left[\begin{array}{c} f^{{+}} \\ [2pt] f^{{-}} \end{array}\right... ...ft[\begin{array}{c} \frac{f+Rv}{2} \\ [2pt] \frac{f-Rv}{2} \end{array}\right].$

Similarly, velocity waves $(v^{+},v^{-})$ are prepared from

according to

$\displaystyle \left[\begin{array}{c} v^{+} \\ [2pt] v^{-} \end{array}\right] = ... ...[\begin{array}{c} \frac{v+f/R}{2} \\ [2pt] \frac{v-f/R}{2} \end{array}\right].$

Power Waves

Basic courses in physics teach us that power is work per unit time, and work is a measure of energy which is typically defined as force times distance. Therefore, power is in physical units of force times distance per unit time, or force times velocity. It therefore should come as no surprise that traveling power waves are defined for strings as follows:

$\begin{displaymath}\begin{array}{rcrl} {\cal P}^{+}(n) &\isdef &&f^{{+}}(n)v^{+}(n) \\ {\cal P}^{-}(n) &\isdef &-&f^{{-}}(n)v^{-}(n) \end{array}\end{displaymath}$

From the Ohm's-law relations $f^{{+}}= Rv^{+}$ and $f^{{-}}=-Rv^{-}$ , we also have

$\begin{eqnarray*} {\cal P}^{+}(n)&=&R[v^{+}(n)]^2 \eqsp [f^{{+}}(n)]^2/R, \\ {\cal P}^{-}(n)&=&R[v^{-}(n)]^2 \eqsp [f^{{-}}(n)]^2/R. \end{eqnarray*}$

Thus, both the left- and right-going components are nonnegative. The sum of the traveling powers at a point gives the total power at that point in the waveguide:

$\displaystyle {\cal P}(t_n,x_m) \isdefs {\cal P}^{+}(n-m) + {\cal P}^{-}(n+m)$

(C.49)

If we had left out the minus sign in the definition of left-going power waves, the sum would instead be a net power flow.

Power waves are important because they correspond to the actual ability of the wave to do work on the outside world, such as on a violin bridge at the end of a string. Because energy is conserved in closed systems, power waves sometimes give a simpler, more fundamental view of wave phenomena, such as in conical acoustic tubes. Also, implementing nonlinear operations such as rounding and saturation in such a way that signal power is not increased, gives suppression of limit cycles and overflow oscillations [432], as discussed in the section on signal scattering.

For example, consider a waveguide having a wave impedance which increases smoothly to the right. A converging cone provides a practical example in the acoustic tube realm. Then since the energy in a traveling wave must be in the wave unless it has been transduced elsewhere, we expect ${\cal P}^{+}$ to propagate unchanged along the waveguide. However, since the wave impedance is increasing, it must be true that force is increasing and velocity is decreasing according to ${\cal P}^{+}= R(v^{+})^2 = (f^{{+}})^2/R$ . Looking only at force or velocity might give us the mistaken impression that the wave is getting stronger (looking at force) or weaker (looking at velocity), when really it was simply sailing along as a fixed amount of energy. This is an example of a transformer action in which force is converted into velocity or vice versa. It is well known that a conical tube acts as if it's open on both ends even though we can plainly see that it is closed on one end. A tempting explanation is that the cone acts as a transformer which exchanges pressure and velocity between the endpoints of the tube, so that a closed end on one side is equivalent to an open end on the other. However, this view is oversimplified because, while spherical pressure waves travel nondispersively in cones, velocity propagation is dispersive [22,50].

Energy Density Waves

The vibrational energy per unit length along the string, or wave energy density [317] is given by the sum of potential and kinetic energy densities:

$\displaystyle W(t,x) \isdefs \underbrace{\frac{1}{2} Ky'^2(t,x)}_{\mbox{potential}} + \underbrace{\frac{1}{2} \epsilon {\dot y}^2(t,x)}_{\mbox{kinetic}}$

(C.50)

Sampling across time and space, and substituting traveling wave components, one can show in a few lines of algebra that the sampled wave energy density is given by

$\displaystyle W(t_n,x_m) \isdefs W^{+}(n-m) + W^{-}(n+m),$

(C.51)

where

$\begin{eqnarray*} W^{+}(n) &=& \frac{{\cal P}^{+}(n)}{c} \,\mathrel{\mathop=}\,\... ...ht]^2 \,\mathrel{\mathop=}\,\frac{\left[f^{{-}}(n)\right]^2}{K}. \end{eqnarray*}$

Thus, traveling power waves (energy per unit time) can be converted to energy density waves (energy per unit length) by simply dividing by , the speed of propagation. Quite naturally, the total wave energy in the string is given by the integral along the string of the energy density:

$\displaystyle {\cal E}(t) \,\mathrel{\mathop=}\,\int_{x=-\infty}^\infty W(t,x)dx \approx \sum_{m = -\infty}^\infty W(t,x_m)X$

(C.52)

In practice, of course, the string length is finite, and the limits of integration are from the

coordinate of the left endpoint to that of the right endpoint, e.g., 0 to

Root-Power Waves

It is sometimes helpful to normalize the wave variables so that signal power is uniformly distributed numerically. This can be especially helpful in fixed-point implementations.

From (C.49), it is clear that power normalization is given by

$\begin{displaymath}\begin{array}{rclrcl} \tilde{f}^{+}&\isdef & f^{{+}}/\sqrt{R}... ...rt{R}\qquad & \tilde{v}^{-}& \isdef & v^{-}\sqrt{R} \end{array}\end{displaymath}$

(C.53)

where we have dropped the common time argument `

' for simplicity. As a result, we obtain

$\begin{displaymath}\begin{array}{rcccl} {\cal P}^{+}& = & f^{{+}}v^{+}&=& \tilde... ...\ &=&(f^{{+}})^2 / R&=& (\tilde{f}^{+})^2 \nonumber \end{array}\end{displaymath}$

and

$\begin{displaymath}\begin{array}{rcccl} {\cal P}^{-}& = & -f^{{-}}v^{-}&=& -\til... ... &=&(f^{{-}})^2 / R&=& (\tilde{f}^{-})^2. \nonumber \end{array}\end{displaymath}$

The normalized wave variables $\tilde{f}^\pm$ and $\tilde{v}^\pm$ behave physically like force and velocity waves, respectively, but they are scaled such that either can be squared to obtain instantaneous signal power. Waveguide networks built using normalized waves have many desirable properties [174,172,432]. One is the obvious numerical advantage of uniformly distributing signal power across available dynamic range in fixed-point implementations. Another is that only in the normalized case can the wave impedances be made time varying without modulating signal power. In other words, use of normalized waves eliminates ``parametric amplification'' effects; signal power is decoupled from parameter changes.

Total Energy in a Rigidly Terminated String

The total energy ${\cal E}$ in a length , rigidly terminated, freely vibrating string can be computed as

$\displaystyle {\cal E}(t)$	$\displaystyle \isdef$	$\displaystyle \int_{0}^L W(t,x)dx = \int_{t_0}^{t_0+2L/c}{\cal P}(\tau,x) d\tau$	(C.54)
	$\displaystyle \approx$	$\displaystyle \sum_{m=0}^{N/2-1} W(t_n,x_m)X = \sum_{n=1}^{N}{\cal P}(t_0 + t_n,x_m) T$	(C.55)

for any $x\in[0,L]$ . Since the energy never decays,

and

are also arbitrary. Thus, because free vibrations of a doubly terminated string must be periodic in time, the total energy equals the integral of power over any period at any point along the string.

Scattering at Impedance Changes

When a traveling wave encounters a change in wave impedance, scattering occurs, i.e., a traveling wave impinging on an impedance discontinuity will partially reflect and partially transmit at the junction in such a way that energy is conserved. This is a classical topic in transmission line theory [295], and it is well covered for acoustic tubes in a variety of references [297,363]. However, for completeness, we derive the basic scattering relations below for plane waves in air, and for longitudinal stress waves in rods.

Plane-Wave Scattering

**Figure C.15:** Plane wave propagation in a medium changing from wave impedance to .
$\includegraphics{eps/planewavescat}$

Consider a plane wave with peak pressure amplitude propagating from wave impedance into a new wave impedance , as shown in Fig.C.15. (Assume and are real and positive.) The physical constraints on the wave are that

pressure must be continuous everywhere, and
velocity in must equal velocity out (the junction has no state).

Since power is pressure times velocity, these constraints imply that signal power is conserved at the junction.^C.5Expressed mathematically, the physical constraints at the junction can be written as follows:

$\begin{eqnarray*} p^+_1+p^-_1 &=& p^+_2\quad\mbox{(pressure continuous across ju... ...+}_1+v^{-}_1 &=& v^{+}_2\quad\mbox{(velocity in = velocity out)} \end{eqnarray*}$

As derived in §C.7.3, we also have the Ohm's law relations:

$\begin{displaymath} \begin{array}{rcrl} p^+_i&=&&R_iv^{+}_i\\ p^-_i&=&-&R_iv^{-}_i \end{array}\end{displaymath}$

These equations determine what happens at the junction.

To obey the physical constraints at the impedance discontinuity, the incident plane-wave must split into a reflected plane wave and a transmitted plane-wave such that pressure is continuous and signal power is conserved. The physical pressure on the left of the junction is , and the physical pressure on the right of the junction is , since according to our set-up.

Scattering Solution

Define the junction pressure and junction velocity by

$\begin{eqnarray*} p_j &\isdef & p^+_1+p^-_1 = p^+_2\quad\mbox{(pressure at junct... ...f & v^{+}_1+v^{-}_1 = v^{+}_2\quad\mbox{(velocity at junction).} \end{eqnarray*}$

Then we can write

$\begin{eqnarray*} p^+_1+p^-_1 &=& p^+_2\;=\;p_j\\ [10pt] \,\,\Rightarrow\,\,R_1v... ...\\ [10pt] \,\,\Rightarrow\,\,2\,R_1v^{+}_1 - R_1 v_j &=& R_2 v_j \end{eqnarray*}$

$\displaystyle \,\,\Rightarrow\,\,\zbox {v_j = \frac{2\,R_1}{R_1 + R_2}v^{+}_1.}$

Note that $v_j=v^{+}_2$ , so we have found the velocity of the transmitted wave. Since $v_j = v^{+}_1+v^{-}_1$ , the velocity of the reflected wave is simply

$\displaystyle v^{-}_1 = v_j - v^{+}_1 = \left[\frac{2\,R_1}{R_1+R_2} - 1\right]v^{+}_1 = \frac{R_1-R_2}{R_1+R_2} v^{+}_1.$

We have solved for the transmitted and reflected velocity waves given the incident wave and the two impedances.

Using the Ohm's law relations, the pressure waves follow easily:

$\begin{eqnarray*} p^+_2 &=& R_2v^{+}_2 = R_2 v_j = \frac{2\,R_2}{R_1+R_2}p^+_1\\ [10pt] p^-_1 &=& -R_1v^{-}_1 = \frac{R_2-R_1}{R_1+R_2} p^+_1 \end{eqnarray*}$

Reflection Coefficient

Define the reflection coefficient of the scattering junction as

$\displaystyle \zbox {\rho = \frac{R_2-R_1}{R_1+R_2} = \frac{\mbox{Impedance Step}}{\mbox{Impedance Sum}}.}$

Then we get the following scattering relations in terms of $\rho$ for pressure waves:

$\begin{eqnarray*} p^+_2 &=& (1+\rho)p^+_1\\ [3pt] p^-_1 &=& \rho\,p^+_1 \end{eqnarray*}$

Signal flow graphs for pressure and velocity are given in Fig.C.16.

**Figure C.16:** Signal flow graph for the pressure and velocity components of a plane wave scattering at an impedance discontinuity :.
$\includegraphics{eps/planewavescatdm2}$

It is a simple exercise to verify that signal power is conserved by checking that $p^+_1v^{+}_1 = p^+_2v^{+}_2 + ( - p^-_1v^{-}_1)$ . (Left-going power is negated to account for its opposite direction-of-travel.)

So far we have only considered a plane wave incident on the left of the junction. Consider now a plane wave incident from the right. For that wave, the impedance steps from to , so the reflection coefficient it ``sees'' is $-\rho$ . By superposition, the signal flow graph for plane waves incident from either side is given by Fig.C.17. Note that the transmission coefficient is one plus the reflection coefficient in either direction. This signal flow graph is often called the ``Kelly-Lochbaum'' scattering junction [297].

**Figure C.17:** Signal flow graph for plane waves incident on either the left or right of an impedance discontinuity. Also shown are delay lines corresponding to sampled traveling plane-wave components propagating on either side of the scattering junction.
$\includegraphics{eps/planewavescatdmfull}$

There are some simple special cases:

$R_2=\infty\,\,\Rightarrow\,\,\rho = 1\quad$ (e.g., rigid wall reflection)
$R_2=0\,\,\Rightarrow\,\,\rho = -1\quad$ (e.g., open-ended tube)
$R_2=R_1\,\,\Rightarrow\,\,\rho = 0\quad$ (no reflection)

Plane-Wave Scattering at an Angle

Figure C.18 shows the more general situation (as compared to Fig.C.15) of a sinusoidal traveling plane wave encountering an impedance discontinuity at some arbitrary angle of incidence, as indicated by the vector wavenumber $\underline{k}_1^+$ . The mathematical details of general sinusoidal plane waves in air and vector wavenumber are reviewed in §B.8.1.

**Figure C.18:** Sinusoidal plane wave scattering at an impedance discontinuity--oblique angle of incidence $\theta _1^+$ .
$\includegraphics{eps/planewavescatangle}$

At the boundary between impedance and , we have, by continuity of pressure,

$\begin{eqnarray*} k_1\sin(\theta_1^+) &=&k_1\sin(\theta_1^-) \;=\;k_2\sin(\theta_2^+) \end{eqnarray*}$

as we will now derive.

Let the impedance change be in the $\underline{x}=(0,y,z)$ plane. Thus, the impedance is for $x\le0$ and for . There are three plane waves to consider:

The incident plane wave with wave vector $\underline{k}_1^+$
The reflected plane wave with wave vector $\underline{k}_1^-$
The transmitted plane wave with wave vector $\underline{k}_2^+$

By continuity, the waves must agree on boundary plane:

$\displaystyle \left<\underline{k}_1^+,\underline{r}\right> = \left<\underline{k}_1^-,\underline{r}\right> = \left<\underline{k}_2^+,\underline{r}\right>$

where $\underline{r}=(0,y,z)$ denotes any vector in the boundary plane. Thus, at

we have

$\displaystyle k_{1y}^+\,y + k_{1z}^+\,z = k_{1y}^-\,y + k_{1z}^-\,z = k_{2y}^+\,y + k_{2z}^+\,z.$

If the incident wave is constant along

, then $k_{1z}^+=0$ , requiring $k_{1z}^- = k_{2z}^+ = 0$ , leaving

$\displaystyle k_{1y}^+\,y = k_{1y}^-\,y =k_{2y}^+\,y$

$\displaystyle \zbox {k_1\sin(\theta_1^+) =k_1\sin(\theta_1^-) =k_2\sin(\theta_2^+)} \protect$

(C.56)

where $\theta$ is defined as zero when traveling in the direction of positive

for the incident ( $\underline{k}_1^+$ ) and transmitted ( $\underline{k}_2^+$ ) wave vector, and along negative

for the reflected ( $\underline{k}_1^-$ ) wave vector (see Fig.C.18).

Reflection and Refraction

The first equality in Eq. $\,$ (C.56) implies that the angle of incidence equals angle of reflection:

$\displaystyle \zbox {\theta_1^+=\theta_1^-} % \isdef\theta_1}$

We now wish to find the wavenumber in medium 2. Let denote the phase velocity in wave impedance :

$\displaystyle c_i = \frac{\omega}{k_i}, \quad i=1,2$

In impedance

, we have in particular

$\displaystyle \omega^2 \eqsp c_2^2 k_2^2 \eqsp c_2^2 \left[(k^+_{2x})^2 + (k^+_{2y})^2\right].$

Solving for $k^+_{2x}$ gives

$\displaystyle k^+_{2x} \eqsp \sqrt{\frac{\omega^2}{c_2^2} - (k^+_{2y})^2} \eqsp \sqrt{\frac{\omega^2}{c_2^2} - k_2^2\sin^2(\theta_2^+)}.$

Since $k_1\sin(\theta_1^+)=k_2\sin(\theta_2^+)$ from above,

$\displaystyle k^+_{2x} \eqsp \sqrt{\frac{\omega^2}{c_2^2} - k_1^2\sin^2(\theta... ...\eqsp \sqrt{\frac{\omega^2}{c_2^2}-\frac{\omega^2}{c_1^2}\sin^2(\theta_1^+)}.$

We have thus derived

$\displaystyle \zbox {k^+_{2x} \eqsp \frac{\omega}{c_2}\sqrt{1 - \frac{c_2^2}{c_1^2}\sin^2(\theta_1^+)}.}$

We earlier established $k^+_{2y} = k^+_{1y}$ (for agreement along the boundary, by pressure continuity). This describes the refraction of the plane wave as it passes through the impedance-change boundary. The refraction angle depends on ratio of phase velocities

. This ratio is often called the index of refraction:

$\displaystyle n \isdef \frac{c_2}{c_1}$

and the relation $k_1\sin(\theta_1^+)=k_2\sin(\theta_2^+)$ is called Snell's Law (of refraction).

Evanescent Wave due to Total Internal Reflection

Note that if $c_1 < c_2 \vert\sin(\theta_1^+)\vert$ , the horizontal component of the wavenumber in medium 2 becomes imaginary. In this case, the wave in medium 2 is said to be evanescent, and the wave in medium 1 undergoes total internal reflection (no power travels from medium 1 to medium 2). The evanescent-wave amplitude decays exponentially to the right and oscillates ``in place'' (like a standing wave). ``Tunneling'' is possible given a medium 3 beyond medium 2 in which wave propagation resumes.

To show explicitly the exponential decay and in-place oscillation in an evanescent wave, express the imaginary wavenumber as $k_x\isdef -j\kappa_x$ . Then we have

$\begin{eqnarray*} p(t,\underline{x}) &=& \cos\left(\omega t - \underline{k}^T\... ...-k_x x}\right\}}}\\ [5pt] &=& e^{-k_x x} \cos(\omega t - k_y y). \end{eqnarray*}$

Thus, an imaginary wavenumber corresponds to an exponentially decaying evanescent wave. Note that the time dependence (cosine term) applies to all points to the right of the boundary. Since evanescent waves do not really ``propagate,'' it is perhaps better to speak of an ``evanescent acoustic field'' or ``evanescent standing wave'' instead of ``evanescent waves''.

For more on the physics of evanescent waves and tunneling, see [295].

Longitudinal Waves in Rods

In this section, elementary scattering relations will be derived for the case of longitudinal force and velocity waves in an ideal string or rod. In solids, force-density waves are referred to as stress waves [169,261]. Longitudinal stress waves in strings and rods have units of (compressive) force per unit area and are analogous to longitudinal pressure waves in acoustic tubes.

**Figure:** A waveguide section between two partial sections. a) Physical picture indicating traveling waves in a continuous medium whose wave impedance changes from to to . b) Digital simulation diagram for the same situation. The section propagation delay is denoted as $z^{-T}$ . The behavior at an impedance discontinuity is characterized by a lossless splitting of an incoming wave into transmitted and reflected components.
$\includegraphics[width=\twidth]{eps/Fwgfs}$

A single waveguide section between two partial sections is shown in Fig.C.19. The sections are numbered 0 through from left to right, and their wave impedances are , , and , respectively. Such a rod might be constructed, for example, using three different materials having three different densities. In the th section, there are two stress traveling waves: $f^{{+}}_i$ traveling to the right at speed , and $f^{{-}}_i$ traveling to the left at speed . To minimize the numerical dynamic range, velocity waves may be chosen instead when .

As in the case of transverse waves (see the derivation of (C.46)), the traveling longitudinal plane waves in each section satisfy [169,261]

$\begin{displaymath}\begin{array}{rcrl} f^{{+}}_i(t)&=&&R_iv^{+}_i(t) \\ f^{{-}}_i(t)&=&-&R_iv^{-}_i(t) \end{array}\end{displaymath}$

(C.57)

where the wave impedance is now $R_i=\sqrt{E\rho}$ , with $\rho$ being the mass density, and

being the Young's modulus of the medium (defined as the stress over the strain, where strain means relative displacement--see §B.5.1) [169,261]. As before, velocity $v_i=v^{+}_i+v^{-}_i$ is defined as positive to the right, and $f^{{+}}_i$ is the right-going traveling-wave component of the stress, and it is positive when the rod is locally compressed.

If the wave impedance is constant, the shape of a traveling wave is not altered as it propagates from one end of a rod-section to the other. In this case we need only consider $f^{{+}}_i$ and $f^{{-}}_i$ at one end of each section as a function of time. As shown in Fig.C.19, we define $f^\pm _i(t)$ as the force-wave component at the extreme left of section . Therefore, at the extreme right of section , we have the traveling waves $f^{{+}}_i(t-T)$ and $f^{{-}}_i(t+T)$ , where is the travel time from one end of a section to the other.

For generality, we may allow the wave impedances to vary with time. A number of possibilities exist which satisfy (C.57) in the time-varying case. For the moment, we will assume the traveling waves at the extreme right of section are still given by $f^{{+}}_i(t-T)$ and $f^{{-}}_i(t+T)$ . This definition, however, implies the velocity varies inversely with the wave impedance. As a result, signal energy, being the product of force times velocity, is ``pumped'' into or out of the waveguide by a changing wave impedance. Use of normalized waves $\tilde{f}^\pm _i$ avoids this. However, normalization increases the required number of multiplications, as we will see in §C.8.6 below.

As before, the physical force density (stress) and velocity at the left end of section are obtained by summing the left- and right-going traveling wave components:

$\displaystyle f_i$	$\displaystyle =$	$\displaystyle f^{{+}}_i+ f^{{-}}_i$	(C.58)
$\displaystyle v_i$	$\displaystyle =$	$\displaystyle v^{+}_i+ v^{-}_i$

Let

denote the force at position

and time

in section

, where $x_i\in[0,cT]$ is measured from the extreme left of section

along its axis. Then we have, for example, $f_i(t,0)\isdef f^{{+}}_i(t)+f^{{-}}_i(t)$ and $f_i(t,cT)\isdef f^{{+}}_i(t-T)+f^{{-}}_i(t+T)$ at the boundaries of section

. More generally, within section

, the physical stress may be expressed in terms of its traveling-wave components by

$\displaystyle f_i(t,x_i) \eqsp f^{{+}}_i\left(t-\frac{x_i}{c}\right)+f^{{-}}_i\left(t+\frac{x_i}{c}\right), \quad 0\leq x_i\leq cT.$

Kelly-Lochbaum Scattering Junctions

Conservation of energy and mass dictate that, at the impedance discontinuity, force and velocity variables must be continuous

$\displaystyle f_{i-1}(t,cT)$	$\displaystyle =$	$\displaystyle f_i(t,0)$	(C.59)
$\displaystyle v_{i-1}(t,cT)$	$\displaystyle =$	$\displaystyle v_i(t,0)$

where velocity is defined as positive to the right on both sides of the junction. Force (or stress or pressure) is a scalar while velocity is a vector with both a magnitude and direction (in this case only left or right). Equations (C.57), (C.58), and (C.59) imply the following scattering equations (a derivation is given in the next section for the more general case of

waveguides meeting at a junction):

$\displaystyle f^{{+}}_i(t)$	$\displaystyle =$	$\displaystyle \left[1+k_i(t) \right]f^{{+}}_{i-1}(t-T) - k_i(t) f^{{-}}_i(t)$
$\displaystyle f^{{-}}_{i-1}(t+T)$	$\displaystyle =$	$\displaystyle k_i(t)f^{{+}}_{i-1}(t-T) + \left[1-k_i(t)\right]f^{{-}}_i(t)$	(C.60)

where

$\displaystyle k_i(t) \isdef \frac{ R_i(t)-R_{i-1}(t) }{R_i(t)+R_{i-1}(t) }$

(C.61)

is called the

th reflection coefficient. Since $R_i(t)\geq 0$ , we have $k_i(t)\in[-1,1]$ . It can be shown that if $\vert k_i\vert>1$ , then either

or $R_{i-1}$ is negative, and this implies an active (as opposed to passive) medium. Correspondingly, lattice and ladder recursive digital filters are stable if and only if all reflection coefficients are bounded by

in magnitude [297].

**Figure C.20:** The Kelly-Lochbaum scattering junction.
$\includegraphics[scale=0.9]{eps/Fkl}$

The scattering equations are illustrated in Figs. C.19b and C.20. In linear predictive coding of speech [482], this structure is called the Kelly-Lochbaum scattering junction, and it is one of several types of scattering junction used to implement lattice and ladder digital filter structures (§C.9.4,[297]).

One-Multiply Scattering Junctions

By factoring out in each equation of (C.60), we can write

$\displaystyle f^{{+}}_i(t)$	$\displaystyle =$	$\displaystyle f^{{+}}_{i-1}(t-T) + f_{{\Delta}}(t)$
$\displaystyle f^{{-}}_{i-1}(t+T)$	$\displaystyle =$	$\displaystyle f^{{-}}_i(t) + f_{{\Delta}}(t)$	(C.62)

where

$\displaystyle f_{{\Delta}}(t) \isdef k_i(t)\left[f^{{+}}_{i-1}(t-T) - f^{{-}}_i(t) \right]$

(C.63)

Thus, only one multiplication is actually necessary to compute the transmitted and reflected waves from the incoming waves in the Kelly-Lochbaum junction. This computation is shown in Fig.C.21, and it is known as the one-multiply scattering junction [297].

**Figure C.21:** The one-multiply scattering junction.
$\includegraphics[scale=0.9]{eps/Fom}$

Another one-multiply form is obtained by organizing (C.60) as

$\displaystyle f^{{+}}_i(t)$	$\displaystyle =$	$\displaystyle f^{{-}}_i(t) + \alpha_i(t)\tilde{f_d}(t)$
$\displaystyle f^{{-}}_{i-1}(t+T)$	$\displaystyle =$	$\displaystyle f^{{+}}_i(t) - \tilde{f_d}(t)$	(C.64)

where

$\displaystyle \alpha_i(t)$	$\displaystyle \isdef$	$\displaystyle 1+k_i(t)$
$\displaystyle \tilde{f_d}(t)$	$\displaystyle \isdef$	$\displaystyle f^{{+}}_{i-1}(t-T) - f^{{-}}_i(t).$	(C.65)

As in the previous case, only one multiplication and three additions are required per junction. This one-multiply form generalizes more readily to junctions of more than two waveguides, as we'll see in a later section.

A scattering junction well known in the LPC speech literature but not described here is the so-called two-multiply junction [297] (requiring also two additions). This omission is because the two-multiply junction is not valid as a general, local, physical modeling building block. Its derivation is tied to the reflectively terminated, cascade waveguide chain. In cases where it applies, however, it can be the implementation of choice; for example, in DSP chips having a fast multiply-add instruction, it may be possible to implement the inner loop of the two-multiply, two-add scattering junction using only two instructions.

Normalized Scattering Junctions

**Figure C.22:** The normalized scattering junction.
$\includegraphics[scale=0.9]{eps/scatnlf}$

Using (C.53) to convert to normalized waves $\tilde{f}^\pm$ , the Kelly-Lochbaum junction (C.60) becomes

$\displaystyle \tilde{f}^{+}_i(t)$	$\displaystyle =$	$\displaystyle \sqrt{1-k_i^2(t)}\, \tilde{f}^{+}_{i-1}(t-T) - k_i(t)\, \tilde{f}^{-}_i(t)$
$\displaystyle \tilde{f}^{-}_{i-1}(t+T)$	$\displaystyle =$	$\displaystyle k_i(t)\,\tilde{f}^{+}_{i-1}(t-T) + \sqrt{1-k_i^2(t)}\,\tilde{f}^{-}_i(t)$	(C.66)

as diagrammed in Fig.C.22. This is called the normalized scattering junction [297], although a more precise term would be the ``normalized-wave scattering junction.''

It is interesting to define $\theta_i \isdef \sin^{-1}(k_i)$ , always possible for passive junctions since $-1\leq k_i\leq 1$ , and note that the normalized scattering junction is equivalent to a 2D rotation:

$\displaystyle \tilde{f}^{+}_i(t)$	$\displaystyle =$	$\displaystyle \cos(\theta_i) \, \tilde{f}^{+}_{i-1}(t-T) - \sin(\theta_i) \, \tilde{f}^{-}_i(t)$
$\displaystyle \tilde{f}^{-}_{i-1}(t+T)$	$\displaystyle =$	$\displaystyle \sin(\theta_i)\, \tilde{f}^{+}_{i-1}(t-T) + \cos(\theta_i)\, \tilde{f}^{-}_i(t)$	(C.67)

where, for conciseness of notation, the time-invariant case is written.

While it appears that scattering of normalized waves at a two-port junction requires four multiplies and two additions, it is possible to convert this to three multiplies and three additions using a two-multiply ``transformer'' to power-normalize an ordinary one-multiply junction [432].

The transformer is a lossless two-port defined by [136]

$\displaystyle f^{{+}}_i(t)$	$\displaystyle =$	$\displaystyle g_i\, f^{{+}}_{i-1}(t-T)$
$\displaystyle f^{{-}}_{i-1}(t+T)$	$\displaystyle =$	$\displaystyle \frac{1}{g_i}\,f^{{-}}_i(t).$	(C.68)

The transformer can be thought of as a device which steps the wave impedance to a new value without scattering; instead, the traveling signal power is redistributed among the force and velocity wave variables to satisfy the fundamental relations $f^\pm =\pm Rv^\pm$ (C.57) at the new impedance. An impedance change from $R_{i-1}$ on the left to

on the right is accomplished using

$\displaystyle g_i \isdefs \sqrt\frac{R_i}{R_{i-1}} \eqsp \sqrt\frac{1+k_i(t)}{1-k_i(t)} \protect$

(C.69)

as can be quickly derived by requiring $(f^{{+}}_{i-1})^2/R_{i-1}= (f^{{+}}_i)^2/R_i$ . The parameter

can be interpreted as the ``turns ratio'' since it is the factor by which force is stepped (and the inverse of the velocity step factor).

**Figure C.23:** Three-multiply normalized-wave scattering junction.
$\includegraphics{eps/scatThreeMulNLFx}$

Figure C.23 illustrates a three-multiply normalized-wave scattering junction [432]. The impedance of all waveguides (bidirectional delay lines) may be taken to be . Scattering junctions may then be implemented as a denormalizing transformer $g=\sqrt{R_{i-1}}$ , a one-multiply scattering junction , and a renormalizing transformer $g=1/\sqrt{R_i}$ . Either transformer may be commuted with the junction and combined with the other transformer to give a three-multiply normalized-wave scattering junction. (The transformers are combined on the left in Fig.C.23).

In slightly more detail, a transformer $g=\sqrt{R_{i-1}}$ steps the wave impedance (left-to-right) from to $R=R_{i-1}$ . Equivalently, the normalized force-wave $\tilde{f}^{+}_{i-1}(t)$ is converted unnormalized form $f^{{+}}_{i-1}(t)$ . Next there is a physical scattering from impedance $R_{i-1}$ to (reflection coefficient $k_i=(R_i-R_{i-1})/(R_i+R_{i-1})$ ). The outgoing wave to the right is then normalized by transformer $g=1/\sqrt{R_i}$ to return the wave impedance back to for wave propagation within a normalized-wave delay line to the right. Finally, the right transformer is commuted left and combined with the left transformer to reduce total computational complexity to one multiply and three adds.

It is important to notice that transformer-normalized junctions may have a large dynamic range in practice. For example, if $k_i\in [-1+\epsilon,1-\epsilon]$ , then Eq. $\,$ (C.69) shows that the transformer coefficients may become as large as $\sqrt{2/\epsilon - 1}$ . If $\epsilon$ is the ``machine epsilon,'' i.e., $\epsilon = 2^{-(n-1)}$ for typical -bit two's complement arithmetic normalized to lie in , then the dynamic range of the transformer coefficients is bounded by $\sqrt{2^n-1}\approx 2^{n/2}$ . Thus, while transformer-normalized junctions trade a multiply for an add, they require up to % more bits of dynamic range within the junction adders. On the other hand, it is very nice to have normalized waves (unit wave impedance) throughout the digital waveguide network, thereby limiting the required dynamic range to root physical power in all propagation paths.

Junction Passivity

In fixed-point implementations, the round-off error and other nonlinear operations should be confined when possible to physically meaningful wave variables. When this is done, it is easy to ensure that signal power is not increased by the nonlinear operations. In other words, nonlinear operations such as rounding can be made passive. Since signal power is proportional to the square of the wave variables, all we need to do is make sure amplitude is never increased by the nonlinearity. In the case of rounding, magnitude truncation, sometimes called ``rounding toward zero,'' is one way to achieve passive rounding. However, magnitude truncation can attenuate the signal excessively in low-precision implementations and in scattering-intensive applications such as the digital waveguide mesh [518]. Another option is error power feedback in which case the cumulative round-off error power averages to zero over time.

A valuable byproduct of passive arithmetic is the suppression of limit cycles and overflow oscillations [432]. Formally, the signal power of a conceptually infinite-precision implementation can be taken as a Lyapunov function bounding the squared amplitude of the finite-precision implementation.

The Kelly-Lochbaum and one-multiply scattering junctions are structurally lossless [500,460] (see also §C.17) because they have only one parameter (or $\alpha_i$ ), and all quantizations of the parameter within the allowed interval (or ) correspond to lossless scattering.^C.6

In the Kelly-Lochbaum and one-multiply scattering junctions, because they are structurally lossless, we need only double the number of bits at the output of each multiplier, and add one bit of extended dynamic range at the output of each two-input adder. The final outgoing waves are thereby exactly computed before they are finally rounded to the working precision and/or clipped to the maximum representable magnitude.

For the Kelly-Lochbaum scattering junction, given -bit signal samples and -bit reflection coefficients, the reflection and transmission multipliers produce and bits, respectively, and each of the two additions adds one more bit. Thus, the intermediate word length required is bits, and this must be rounded without amplification down to bits for the final outgoing samples. A similar analysis gives also that the one-multiply scattering junction needs bits for the extended precision intermediate results before final rounding and/or clipping.

To formally show that magnitude truncation is sufficient to suppress overflow oscillations and limit cycles in waveguide networks built using structurally lossless scattering junctions, we can look at the signal power entering and leaving the junction. A junction is passive if the power flowing away from it does not exceed the power flowing into it. The total power flowing away from the th junction is bounded by the incoming power if

$\displaystyle \underbrace{\frac{[f^{{+}}_i(t)]^2}{R_i(t)} + \frac{[f^{{-}}_{i-1... ...t-T)]^2}{R_{i-1}(t)} + \frac{[f^{{-}}_i(t)]^2}{R_i(t)}}_{\mbox{incoming power}}$

(C.70)

Let ${\hat f}$ denote the finite-precision version of

. Then a sufficient condition for junction passivity is

$\displaystyle \left\vert{\hat f}^{+}_i(t)\right\vert$	$\displaystyle \leq$	$\displaystyle \left\vert f^{{+}}_i(t)\right\vert$	(C.71)
$\displaystyle \left\vert{\hat f}^{-}_{i-1}(t+T)\right\vert$	$\displaystyle \leq$	$\displaystyle \left\vert f^{{-}}_{i-1}(t+T)\right\vert$	(C.72)

Thus, if the junction computations do not increase either of the output force amplitudes, no signal power is created. An analogous conclusion is reached for velocity scattering junctions.

Unlike the structurally lossless cases, the (four-multiply) normalized scattering junction has two parameters, $s_i\isdeftext k_i$ and $c_i\isdeftext \sqrt{1-k_i^2}$ , and these can ``get out of synch'' in the presence of quantization. Specifically, let ${\hat s}_i \isdeftext s_i - \epsilon_s$ denote the quantized value of , and let ${\hat c}_i\isdeftext c_i -\epsilon_c$ denote the quantized value of . Then it is no longer the case in general that ${\hat s}^2_i + {\hat c}^2_i= 1$ . As a result, the normalized scattering junction is not structurally lossless in the presence of coefficient quantization. A few lines of algebra shows that a passive rounding rule for the normalized junction must depend on the sign of the wave variable being computed, the sign of the coefficient quantization error, and the sign of at least one of the two incoming traveling waves. We can assume one of the coefficients is exact for passivity purposes, so assume $\epsilon_s=0$ and define ${\hat c}_i=\left\lfloor \sqrt{1-s^2_i}\right\rfloor$ , where $\left\lfloor x\right\rfloor$ denotes largest quantized value less than or equal to . In this case we have $\epsilon_c\geq 0$ . Therefore,

$\displaystyle {\hat f}^{+}_i= {\hat c}_if^{{+}}_{i-1}- s_i f^{{-}}_i= f^{{+}}_i- \epsilon_c f^{{+}}_{i-1}$

and a passive rounding rule which guarantees $\vert{\hat f}^{+}_i\vert \leq \vert f^{{+}}_i\vert$ need only look at the sign bits of ${\hat f}^{+}_i$ and $f^{{+}}_{i-1}$ .

The three-multiply normalized scattering junction is easier to ``passify.'' While the transformer is not structurally lossless, its simplicity allows it to be made passive simply by using non-amplifying rounding on both of its coefficients as well as on its output wave variables. (The transformer is passive when the product of its coefficients has magnitude less than or equal to .) Since there are no additions following the transformer multiplies, double-precision adders are not needed. However, precision and a half is needed in the junction adders to accommodate the worst-case increased dynamic range. Since the one-multiply junction is structurally lossless, the overall junction is passive if non-amplifying rounding is applied to , , and the outgoing wave variables from the transformer and from the one-multiply junction.

In summary, a general means of obtaining passive waveguide junctions is to compute exact results internally, and apply saturation (clipping on overflow) and magnitude truncation (truncation toward zero) to the final outgoing wave variables. Because the Kelly-Lochbaum and one-multiply junctions are structurally lossless, exact intermediate results are obtainable using extended internal precision. For the (four-multiply) normalized scattering junction, a passive rounding rule can be developed based on two sign bits. For the three-multiply normalized scattering junction, it is sufficient to apply magnitude truncation to the transformer coefficients and all outgoing wave variables.

Digital Waveguide Filters

A Digital Waveguide Filter (DWF) can be defined as any digital filter built by stringing together digital waveguides.^C.7 In the case of a cascade combination of unit-length waveguide sections, they are essentially the same as microwave filters [372] and unit element filters [136] (Appendix F). This section shows how DWFs constructed as a cascade chain of waveguide sections, and reflectively terminated, can be transformed via elementary manipulations to well known ladder and lattice digital filters, such as those used in speech modeling [297,363].

DWFs may be derived conceptually by sampling the unidirectional traveling waves in a system of ideal, lossless waveguides. Sampling is across time and space. Thus, variables in a DWF structure are equal exactly (at the sampling times and positions, to within numerical precision) to variables propagating in the physical medium in an interconnection of uniform 1D waveguides.

Ladder Waveguide Filters

An important class of DWFs can be constructed as a linear cascade chain of digital waveguide sections, as shown in Fig.C.24 (a slightly abstracted version of Fig.C.19). Each block labeled stands for a state-less Kelly-Lochbaum scattering junction (§C.8.4) interfacing between wave impedances $R_{i-1}$ and . We call this a ladder waveguide filter structure. It is an exact bandlimited discrete-time model of a sequence of wave impedances , as used in the Kelly-Lochbaum piecewise-cylindrical acoustic-tube model for the vocal tract in the context of voice synthesis (Fig.6.2). The delays between scattering junctions along both the top and bottom signal paths make it possible to compute all of the scattering junctions in parallel--a property that is becoming increasingly valuable in signal-processing architectures.

figure[htbp] $\includegraphics{eps/DWF}$

To transform the DWF of Fig.C.24 to a conventional ladder digital filter structure, as used in speech modeling [297,363], we need to (1) terminate on the right with a pure reflection and (2) eliminate the delays along the top signal path. We will do this in stages so as to point out a valuable intermediate case.

Reflectively Terminated Waveguide Filters

A reflecting termination occurs when the wave impedance of the next section is either zero or infinity (§C.8.1). Since Fig.C.24 is labeled for force waves, an infinite terminating impedance on the right results in

$\displaystyle f^{{-}}_{i+1}(t+T) \eqsp f^{{+}}_{i+1}(t-T).$

Similarly, a zero terminating impedance on the right gives

$\displaystyle f^{{-}}_{i+1}(t+T) \eqsp -f^{{+}}_{i+1}(t-T).$

For velocity waves, the signs are interchanged. Thus, a reflectively terminated ladder digital waveguide corresponds to a final feedback with gain $\pm 1$ at the end of the ladder. Such a termination will be used below to derive conventional ladder/lattice filters below.

Half-Rate Ladder Waveguide Filters

The delays preceding the two inputs to each scattering junction can be ``pushed'' into the junction and then ``pulled'' out to the outputs and combine with the delays there. (This is easy to show using the Kelly-Lochbaum scattering junction derived in §C.8.4.) By performing this operation on every other section in the DWF chain, the half-rate ladder waveguide filter shown in Fig.C.25 is obtained [432].

figure[htbp] $\includegraphics{eps/DWFHR}$

This ``half-rate'' ladder DWF is so-called because its sampling rate can be cut in half due to each delay being two-samples long. It has advantages worth considering, which will become clear after we have derived conventional ladder digital filters below. Note that now pairs of scattering junctions can be computed in parallel, and an exact physical interpretation remains. That is, it can still be used as a general purpose physical modeling building block in this more efficient (half-rate) form.

Note that when the sampling rate is halved, the physical wave variables (computed by summing two traveling-wave components) are at the same time only every other spatial sample. In particular, the physical transverse forces on the right side of scattering junctions and in Fig.C.25 are

$\begin{eqnarray*} f_i(t+T)&=&f^{{+}}_i(t+T)+f^{{-}}_i(t+T)\\ [5pt] f_{i+1}(t)&=&f^{{+}}_{i+1}(t)+f^{{-}}_{i+1}(t) \end{eqnarray*}$

respectively. In the half-rate case, adjacent spatial samples are separated in time by half a temporal sample . If physical variables are needed only for even-numbered (or odd-numbered) spatial samples, then there is no relative time skew, but more generally, things like collision detection, such as for slap-bass string-models (§9.1.6), can be affected. In summary, the half-rate ladder waveguide filter has an alternating half-sample time skew from section to section when used as a physical modeling building block.

Conventional Ladder Filters

Given a reflecting termination on the right, the half-rate DWF chain of Fig.C.25 can be reduced further to the conventional ladder/lattice filter structure shown in Fig.C.26.

figure[htbp] $\includegraphics{eps/DWFHRT}$

To make a standard ladder/lattice filter, the sampling rate is cut in half (i.e., replace $z^{-2}$ by $z^{-1}$ ), and the scattering junctions are typically implemented in one-multiply form (§C.8.5) or normalized form (§C.8.6), etc. Conventionally, if the graph of the scattering junction is nonplanar, as it is for the one-multiply junction, the filter is called a lattice filter; it is called a ladder filter when the graph is planar, as it is for normalized and Kelly-Lochbaum scattering junctions. For all-pole transfer functions $H(z)=F^{-}_0(z)/F^{+}_0(z)$ , the Durbin recursion can be used to compute the reflection coefficients from the desired transfer-function denominator polynomial coefficients [449]. To implement arbitrary transfer-function zeros, a linear combination of delay-element outputs is formed using weights that are called ``tap parameters'' [173,297].

To create Fig.C.26 from Fig.C.24, all delays along the top rail are pushed to the right until they have all been worked around to the bottom rail. In the end, each bottom-rail delay becomes seconds instead of seconds. Such an operation is possible because of the termination at the right by an infinite (or zero) wave impedance. Note that there is a progressive one-sample time advance from section to section. The time skews for the right-going (or left-going) traveling waves can be determined simply by considering how many missing (or extra) delays there are between that signal and the unshifted signals at the far left.

Due to the reflecting termination, conventional lattice filters cannot be extended to the right in any physically meaningful way. Also, creating network topologies more complex than a simple linear cascade (or acyclic tree) of waveguide sections is not immediately possible because of the delay-free path along the top rail. In particular, the output $f^{{+}}_M$ cannot be fed back to the input $f^{{+}}_0$ . Nevertheless, as we have derived, there is an exact physical interpretation (with time skew) for the conventional ladder/lattice digital filter.

Power-Normalized Waveguide Filters

Above, we adopted the convention that the time variation of the wave impedance did not alter the traveling force waves $f^\pm _i$ . In this case, the power represented by a traveling force wave is modulated by the changing wave impedance as it propagates. The signal power becomes inversely proportional to wave impedance:

$\displaystyle {\cal I}_i(t,x) = {\cal I}^{+}_i(t,x)+{\cal I}^{-}_i(t,x) = \frac{[f^{{+}}_i(t,x)]^2-[f^{{-}}_i(t,x)]^2}{R_i(t)}$

In some applications (such as time-varying waveguide reverberation [430]), it may be preferable to compensate for the power modulation so that changes in the wave impedances of the waveguides do not affect the power of the signals propagating within.

In [432,433], three methods are discussed for making signal power invariant with respect to time-varying branch impedances:

The normalized waveguide scheme compensates for power modulation by scaling the signals leaving the delays so as to give them the same power coming out as they had going in. It requires two additional scaling multipliers per waveguide junction.
The normalized wave approach [297] propagates rms-normalized waves in the waveguide. In this case, each delay-line contains $\tilde{f}^{+}_i(t,x) = f^{{+}}_i(t,x)/\sqrt{R_i(t)}$ and $\tilde{f}^{-}_i(t,x) = f^{{-}}_i(t,x)/\sqrt{R_i(t)}$ . In this case, the power stored in the delays does not change when the wave impedance changes. This is the basis of the normalized ladder filter (NLF) [174,297]. Unfortunately, four multiplications are obtained at each scattering junction.
The transformer-normalized waveguide approach changes the wave impedance at the output of the delay back to what it was at the time it entered the delay using a ``transformer'' (defined in §C.16).

The transformer-normalized DWF junction is shown in Fig.C.27 [432]. As derived in §C.16, the transformer ``turns ratio'' is given by

$\displaystyle g_i \eqsp \sqrt{\frac{1-k_i}{1+k_i}}.$

We can now modulate a single scattering junction, even in arbitrary network topologies, by inserting a transformer immediately to the left or right of the junction. Conceptually, the wave impedance is not changed over the delay-line portion of the waveguide section; instead, it is changed to the new time-varying value just before (or after) it meets the junction. When velocity is the wave variable, the coefficients

and $g_i^{-1}$ in Fig.C.27 are swapped (or inverted).

**Figure:** Transformer-normalized one-multiply scattering junction, equivalent to the normalized scattering junction shown in Fig.C.22
$\includegraphics{eps/tnscat}$

So, as in the normalized waveguide case, for the price of two extra multiplies per section, we can implement time-varying digital filters which do not modulate stored signal energy. Moreover, transformers enable the scattering junctions to be varied independently, without having to propagate time-varying impedance ratios throughout the waveguide network.

It can be shown [433] that cascade waveguide chains built using transformer-normalized waveguides are equivalent to those using normalized-wave junctions. Thus, the transformer-normalized DWF in Fig.C.27 and the wave-normalized DWF in Fig.C.22 are equivalent. One simple proof is to start with a transformer (§C.16) and a Kelly-Lochbaum junction (§C.8.4), move the transformer scale factors inside the junction, combine terms, and arrive at Fig.C.22. One practical benefit of this equivalence is that the normalized ladder filter (NLF) can be implemented using only three multiplies and three additions instead of the usual four multiplies and two additions.

The transformer-normalized scattering junction is also the basis of the digital waveguide oscillator (§C.17).

``Traveling Waves'' in Lumped Systems

One of the topics in classical network theory is the reflection and transmission, or scattering formulation for lumped networks [35]. Lumped scattering theory also serves as the starting point for deriving wave digital filters (the subject of Appendix F). In this formulation, forces (voltages) and velocities (currents) are replaced by so-called wave variables

$\begin{eqnarray*} f^{{+}}(t) &\isdef & \frac{f(t) + R_0v(t)}{2} \\ f^{{-}}(t) &\isdef & \frac{f(t) - R_0v(t)}{2} \end{eqnarray*}$

where is an arbitrary reference impedance. Since the above wave variables have dimensions of force, they are specifically force waves. The corresponding velocity waves are

$\begin{eqnarray*} v^{+}(t) &\isdef & \frac{1}{2}[v(t) + f(t)/R_0], \\ v^{-}(t) &\isdef & \frac{1}{2}[v(t) - f(t)/R_0]. \end{eqnarray*}$

Dropping the time argument, since it is always `(t)', we see that

$\displaystyle f^{{+}}$	$\displaystyle =$	$\displaystyle R_0v^{+}$
$\displaystyle f^{{-}}$	$\displaystyle =$	$\displaystyle -R_0v^{-}\protect$	(C.73)

and

$\displaystyle f$	$\displaystyle =$	$\displaystyle f^{{+}}+ f^{{-}}$
$\displaystyle v$	$\displaystyle =$	$\displaystyle v^{+}+ v^{-}\protect$	(C.74)

These are the basic relations for traveling waves in an ideal medium such as an ideal vibrating string or acoustic tube. Using voltage and current gives elementary transmission line theory.

Reflectance of an Impedance

Let denote the driving-point impedance of an arbitrary continuous-time LTI system. Then, by definition, where and denote the Laplace transforms of the applied force and resulting velocity, respectively. The wave variable decomposition in Eq. $\,$ (C.74) gives

$\displaystyle F(s)$	$\displaystyle =$	$\displaystyle R(s) V(s)$
$\displaystyle \,\,\Rightarrow\,\,F^{+}(s) + F^{-}(s)$	$\displaystyle =$	$\displaystyle R(s) \left[V^{+}(s) + V^{-}(s)\right]$
	$\displaystyle =$	$\displaystyle R(s) \left[\frac{F^{+}(s) - F^{-}(s)}{R_0}\right]$
$\displaystyle \,\,\Rightarrow\,\,F^{-}(s) \left[\frac{R(s)}{R_0}+1\right]$	$\displaystyle =$	$\displaystyle F^{+}(s) \left[\frac{R(s)}{R_0}-1\right]$
$\displaystyle \,\,\Rightarrow\,\,F^{-}(s)$	$\displaystyle =$	$\displaystyle F^{+}(s) \left[\frac{R(s)-R_0}{R(s)+R_0}\right]$
	$\displaystyle \isdef$	$\displaystyle F^{+}(s)\, \hat{\rho}_f(s) \protect$	(C.75)

We may call $\hat{\rho}_f(s)$ the reflectance of impedance

relative to

. For example, if a transmission line with characteristic impedance

were terminated in a lumped impedance

, the reflection transfer function at the termination, as seen from the end of the transmission line, would be $\hat{\rho}_f(s)$ .

We are working with reflectance for force waves. Using the elementary relations Eq. $\,$ (C.73), i.e., $F^{+}(s) = R_0V^{+}(s)$ and $F^{-}(s) = -R_0V^{-}(s)$ , we immediately obtain the corresponding velocity-wave reflectance:

$\displaystyle \hat{\rho}_v(s) \isdefs \frac{V^{-}(s)}{V^{+}(s)} \eqsp \frac{-F^... .../R_0}{F^{+}(s)/R_0} \eqsp - \frac{F^{-}(s)}{F^{+}(s)} \eqsp - \hat{\rho}_f(s)$

Properties of Passive Impedances

It is well known that a real impedance (in Ohms, for example) is passive so long as $R\geq 0$ . A passive impedance cannot create energy. On the other hand, if , the impedance is active and has some energy source. The concept of passivity can be extended to complex frequency-dependent impedances $R(j\omega)$ as well: A complex impedance $R(j\omega)$ is passive if is positive real, where is the Laplace-transform variable. The positive-real property is discussed in §C.11.2 below.

This section explores some implications of the positive real condition for passive impedances. Specifically, §C.11.1 considers the nature of waves reflecting from a passive impedance in general, looking at the reflection transfer function, or reflectance, of a passive impedance. To provide further details, Section C.11.2 derives some mathematical properties of positive real functions, particularly for the discrete-time case. Application examples appear in §9.2.1 and §9.2.1.

Passive Reflectances

From Eq. $\,$ (C.75), we have that the reflectance seen at a continuous-time impedance is given for force waves by

$\displaystyle \hat{\rho}_f(s) \isdefs \frac{F^{-}(s)}{F^{+}(s)} \eqsp \frac{R(s)-R_0}{R(s)+R_0} \protect$

(C.76)

where

is the wave impedance connected to the impedance

, and the corresponding velocity reflectance is $\hat{\rho}_v(s)= -\hat{\rho}_f(s)$ . As mentioned above, all passive impedances are positive real. As shown in §C.11.2,

is positive real if and only if $\hat{\rho}_f(s)$ is stable and has magnitude less than or equal to

on the $j\omega$ axis (and hence over the entire left-half plane, by the maximum modulus theorem), i.e.,

$\displaystyle \left\vert\hat{\rho}_f(s)\right\vert \leq 1,$ re $\displaystyle \left\{s\right\} \leq 0. \protect$

(C.77)

In particular, $\left\vert\hat{\rho}_f(j\omega)\right\vert \leq 1$ for all radian frequencies $\omega\in(-\infty,\infty)$ . Any stable $\hat{\rho}_f(s)$ satisfying Eq. $\,$ (C.77) may be called a passive reflectance.

If the impedance goes to infinity (becomes rigid), then $\hat{\rho}_f(s)$ approaches , a result which agrees with an analysis of rigid string terminations (p. ). Similarly, when the impedance goes to zero, $\hat{\rho}_f(s)$ becomes , which agrees with the physics of a string with a free end. In acoustic stringed instruments, bridges are typically quite rigid, so that $\hat{\rho}_f(j\omega)\approx 1$ for all $\omega$ . If a body resonance is strongly coupled through the bridge, $\vert\hat{\rho}_f(j\omega_c)\vert$ can be significantly smaller than 1 at the resonant frequency $\omega_c$ .

Solving for in Eq. $\,$ (C.77), we can characterize every impedance in terms of its reflectance:

$\displaystyle R(s) = R_0\frac{1+\hat{\rho}_f(s)}{1-\hat{\rho}_f(s)}$

Rewriting Eq. $\,$ (C.76) in the form

$\displaystyle \hat{\rho}_f(s) \eqsp \frac{\dfrac{R(s)}{R_0}-1}{\dfrac{R(s)}{R_0}+1},$

we see that the reflectance is determined by the ratio of the ``new impedance''

to the ``old'' impedance

in which the incoming waves travel. In other words, the incoming waves see the wave impedance ``step'' from

, which results in a ``scattering'' of the incident wave into reflected and transmitted components, as discussed in §C.8. The reflection and transmission coefficients depend on frequency when $R(j\omega)$ is not constant with respect to $\omega$ .

In the discrete-time case, which may be related to the continuous-time case by the bilinear transform (§7.3.2), we have the same basic relations, but in the plane:

$\displaystyle \hat{\rho}_f(z)$	$\displaystyle \isdef$	$\displaystyle \frac{F^{-}(z)}{F^{+}(z)} \eqsp \frac{R(z)-R_0}{R(z)+R_0}$
$\displaystyle R(z)$	$\displaystyle =$	$\displaystyle R_0\frac{1+\hat{\rho}_f(z)}{1-\hat{\rho}_f(z)}$
$\displaystyle \Gamma(z)$	$\displaystyle =$	$\displaystyle \Gamma _0\frac{1-\hat{\rho}_f(z)}{1+\hat{\rho}_f(z)} \protect$	(C.78)

where $\Gamma\isdef 1/R$ denotes admittance, with

$\displaystyle \left\vert\hat{\rho}_f(z)\right\vert \leq 1, \quad \left\vert z\right\vert \leq 1. \protect$

(C.79)

Mathematically, any stable transfer function having these properties may be called a Schur function. Thus, the discrete-time reflectance $\hat{\rho}_f(z)$ of an impedance

is a Schur function if and only if the impedance is passive (positive real).

Note that Eq. $\,$ (C.79) may be obtained from the general formula for scattering at a loaded waveguide junction for the case of a single waveguide () terminated by a lumped load (§C.12).

In the limit as damping goes to zero (all poles of converge to the unit circle), the reflectance $\hat{\rho}_f(z)$ becomes a digital allpass filter. Similarly, $\hat{\rho}_f(s)$ becomes a continuous-time allpass filter as the poles of approach the $j\omega$ axis.

Recalling that a lossless impedance is called a reactance (§7.1), we can say that every reactance gives rise to an allpass reflectance. Thus, for example, waves reflecting off a mass at the end of a vibrating string will be allpass filtered, because the driving-point impedance of a mass () is a pure reactance. In particular, the force-wave reflectance of a mass terminating an ideal string having wave impedance is $\hat{\rho}_f(s)= (ms-R_0)/(ms+R_0)$ , which is a continuous-time allpass filter having a pole at and a zero at .

It is intuitively reasonable that a passive reflection gain cannot exceed at any frequency (i.e., the reflectance is a Schur filter, as defined in Eq. $\,$ (C.79)). It is also reasonable that lossless reflection would have a gain of 1 (i.e., it is allpass).

Note that reflection filters always have an equal number of poles and zeros, as can be seen from Eq. $\,$ (C.76) above. This property is preserved by the bilinear transform, so it holds in both the continuous- and discrete-time cases.

Reflectance and Transmittance of a Yielding String Termination

**Figure C.28:** Ideal vibrating string terminated on the right by a passive impedance .
$\includegraphics{eps/yieldingterm}$

Consider the special case of a reflection and transmission at a yielding termination, or ``bridge'', of an ideal vibrating string on its right end, as shown in Fig.C.28. Denote the incident and reflected velocity waves by $v^{+}(t)$ and $v^{-}(t)$ , respectively, and similarly denote the force-wave components by $f^{{+}}(t)$ and $f^{{-}}(t)$ . Finally, denote the velocity of the termination itself by $v_b(t)=v^{+}(t)+v^{-}(t)$ , and its force-wave reflectance by

$\displaystyle \hat{\rho}_f(s) \isdefs \frac{F^{-}(s)}{F^{+}(s)} \eqsp -\frac{V^{-}(s)}{V^{+}(s)} \eqsp \frac{R_b(s)-R_0}{R_b(s)+R_0},$

where

denotes the string wave impedance.

The bridge velocity is given by

$\displaystyle V_b(s) \eqsp V^{+}(s) + V^{-}(s),$

so that the bridge velocity transmittance is given by

$\displaystyle \hat{\tau}_v(s) \isdefs \frac{V_b(s)}{V^{+}(s)} \eqsp \frac{V^{+}(s)+V^{-}(s)}{V^{+}(s)} \eqsp 1+\hat{\rho}_v(s) \eqsp 1-\hat{\rho}_f(s),$

and the bridge force transmittance is given by

$\displaystyle \hat{\tau}_f(s) \isdefs \frac{F_b(s)}{F^{+}(s)} \eqsp \frac{R_0[V^{+}(s)-V^{-}(s)]}{R_0V^{+}(s)} \eqsp 1+\hat{\rho}_f(s)$

where the bridge force is defined as ``up'' so that it is given for small displacements by the string tension times minus the string slope at the bridge. (Recall from §C.7.2 that force waves on the string are defined by

where

and

denote the string tension and slope, respectively.

Power-Complementary Reflection and Transmission

We can show that the reflectance and transmittance of the yielding termination are power complementary. That is, the reflected and transmitted signal-power sum to yield the incident signal-power.

The average power incident at the bridge at frequency $\omega$ can be expressed in the frequency domain as $F^{+}(e^{j\omega T})\overline{V^{+}(e^{j\omega T})}$ . The reflected power is then $F^{-}\overline{V^{-}} = -\left\vert\hat{\rho}_f\right\vert^2F^{+}\overline{V^{+}}$ . Removing the minus sign, which can be associated with reversed direction of travel, we obtain that the power reflection frequency response is $\left\vert\hat{\rho}_f\right\vert^2$ , which generalizes by analytic continuation to $\hat{\rho}_f(s)\hat{\rho}_f(-s)$ . The power transmittance is given by

$\displaystyle F_b\overline{V_b} \eqsp (\hat{\tau}_fF^{+})\overline{(1-\hat{\rh... ...ine{V^{+}}) \eqsp (1-\left\vert\hat{\rho}_f\right\vert^2)F^{+}\overline{V^{+}}$

which generalizes to the

plane as

$\displaystyle F_b(s)V_b(-s) = \left[1-\hat{\rho}_f(s)\hat{\rho}_f(-s)\right]F^{+}(s)V^{+}(-s)$

Finally, we see that adding up the reflected and transmitted power yields the incident power:

$\displaystyle -F^{-}(s)V^{-}(-s) + F_b(s)V_b(-s) \eqsp F^{+}(s)V^{+}(-s)$

Positive Real Functions

Any passive driving-point impedance, such as the impedance of a violin bridge, is positive real. Positive real functions have been studied extensively in the continuous-time case in the context of network synthesis [68,524]. Very little, however, seems to be available in the discrete time case. This section (reprinted from [428]) summarizes the main properties of positive real function in the plane (i.e., the discrete-time case).

Definition. A complex valued function of a complex variable is said to be positive real (PR) if

$z\ real\,\implies f(z)\ real$
$\left\vert z\right\vert\geq 1\implies$ re $\left\{f(z)\right\}\geq 0$

We now specialize to the subset of functions representable as a ratio of finite-order polynomials in . This class of ``rational'' functions is the set of all transfer functions of finite-order time-invariant linear systems, and we write to denote a member of this class. We use the convention that stable, minimum phase systems are analytic and nonzero in the strict outer disk.^C.8 Condition (1) implies that for to be PR, the polynomial coefficients must be real, and therefore complex poles and zeros must exist in conjugate pairs. We assume from this point on that $H(z)\neq 0$ satisfies (1). From (2) we derive the facts below.

Property 1. A real rational function is PR iff $\left\vert z\right\vert\geq 1 \implies \left\vert\angle{H(z)}\right\vert\leq \frac{\pi}{2}$ .

Proof. Expressing in polar form gives

$\begin{eqnarray*} \mbox{re}\left\{H(z)\right\} &=& \mbox{re}\left\{\left\vert H... ...er} \quad \left\vert\angle{H(z)}\right\vert\leq \frac{\pi}{2}, \end{eqnarray*}$

since the zeros of are isolated. $\Box$

Property 2. is PR iff is PR.

Proof. Assuming is PR, we have by Property 1,

$\displaystyle \left\vert\angle{H^{-1}(z)}\right\vert = \left\vert-\angle{H(z)}\... ...angle{H(z)}\right\vert \leq \frac{\pi}{2},\quad \left\vert z\right\vert\geq 1.$

$\Box$

Property 3. A PR function is analytic and nonzero in the strict outer disk.

Proof. (By contradiction)

Without loss of generality, we treat only $n^{th}$ order polynomials

$\displaystyle \alpha_0z^n + \alpha_1z^{n-1} + \cdots + \alpha_{n-1}z + \alpha_n$

which are nondegenerate in the sense that $\alpha_0,\alpha_n\neq 0$ . Since facts about $\alpha_0 H(z)$ are readily deduced from facts about

, we set $\alpha_0 = 1$ at no great loss.

The general (normalized) causal, finite-order, linear, time-invariant transfer function may be written

$\displaystyle H(z)$	$\displaystyle =$	$\displaystyle z^{-\nu}\frac{b(z)}{a(z)}$
	$\displaystyle =$	$\displaystyle z^{-\nu}\frac{1 + b_1 z^{-1} + \cdots + b_M z^{-M} }{ 1 + a_1 z^{-1} + \cdots + a_N z^{-N} }$
	$\displaystyle =$	$\displaystyle z^{-\nu}\frac{ \prod_{i=1}^M (1 - q_i z^{-1}) }{ \prod_{i=1}^N (1 - p_i z^{-1}) }$
	$\displaystyle =$	$\displaystyle z^{-\nu}\sum_{i=1}^{N_d} \sum_{j=1}^{\mu_i} \frac{z\,K_{i,j} }{ (z-p_i)^j} , \qquad \nu\geq 0,$	(C.80)

where

is the number of distinct poles, each of multiplicity $\mu_i$ ,and

$\displaystyle \sum_{i=1}^{N_d}\mu_i = \max\{N,M\}.$

Suppose there is a pole of multiplicity outside the unit circle. Without loss of generality, we may set $\mu_1=m$ , and $p_1 = R\,e^{j\phi}$ with . Then for near , we have

$\begin{eqnarray*} z^\nu H(z) &=&\frac{ z\,K_{1,m} }{ (z-R\,e^{j\phi})^m} + \fr... ...+ \cdots\\ & \approx \frac{z\, K_{1,m} }{ (z-R\,e^{j\phi})^m}. \end{eqnarray*}$

Consider the circular neighborhood of radius $\rho$ described by $z=R\,e^{j\phi}+ \rho\,e^{j\psi},\;-\pi\leq \psi<\pi$ . Since we may choose $\rho < R-1$ so that all points in this neighborhood lie outside the unit circle. If we write the residue of the factor $(z-R\,e^{j\phi})^m$ in polar form as $K_{1,m}=C\,e^{j\xi}$ , then we have, for sufficiently small $\rho$ ,

$\displaystyle z^\nu H(z) \approx \frac{K_{1,m}R\,e^{j\phi}}{ (z-R\,e^{j\phi})^m... ...\phi}}{ \rho^m\,e^{jm\psi}} = \frac{C\,R}{ \rho^m}\, e^{j(\phi + \xi - m\psi)}.$

(C.81)

Therefore, approaching the pole $R\,e^{j\phi}$ at an angle $\psi$ gives

$\displaystyle \lim_{\rho\to 0}\left\vert\angle{H(R\,e^{j\phi}+\rho\,e^{j\psi})}... ... = \left\vert\phi(1-\nu) + \xi - m\psi \right\vert, \qquad -\pi \leq \psi <\pi$

which cannot be confined to satisfy Property 1 regardless of the value of the residue angle $\xi$ , or the pole angle $\phi$ (

cannot be zero by hypothesis). We thus conclude that a PR function

can have no poles in the outer disk. By Property 2, we conclude that positive real functions must be minimum phase. $\Box$

Corollary. In equation Eq. $\,$ (C.80), $\nu=0$ .

Proof. If $\nu>0$ , then there are $\nu$ poles at infinity. As $\left\vert z\right\vert\to \infty$ , $H(z)\to z^{-\nu} \implies \left\vert\angle{H(z)}\right\vert \to \left\vert\nu\angle{z}\right\vert$ , we must have $\nu=0$ . $\Box$

Corollary. The log-magnitude of a PR function has zero mean on the unit circle.

This is a general property of stable, minimum-phase transfer functions which follows immediately from the argument principle [297,326].

Corollary. A rational PR function has an equal number of poles and zeros all of which are in the unit disk.

This really a convention for numbering poles and zeros. In Eq. $\,$ (C.80), we have $\nu=0$ , and all poles and zeros inside the unit disk. Now, if then we have extra poles at induced by the numerator. If , then zeros at the origin appear from the denominator.

Corollary. Every pole on the unit circle of a positive real function must be simple with a real and positive residue.

Proof. We repeat the previous argument using a semicircular neighborhood of radius $\rho$ about the point $p_1 = e^{j\phi}$ to obtain

$\displaystyle \lim_{\rho\to 0}\left\vert\angle{H(e^{j\phi} + \rho e^{j\psi})}\r... ...m\psi \right\vert,\qquad \phi-\frac{\pi}{2} \leq \psi\leq \phi + \frac{\pi}{2}.$

(C.82)

In order to have $\left\vert\angle{H(z)}\right\vert\leq \pi/2$ near this pole, it is necessary that

and $\xi=0$ . $\Box$

Corollary. If is PR with a zero at $z=q_1=e^{j\phi}$ , then

$\displaystyle H^\prime(z) \isdef \frac{H(z)}{ (1-q_1 z^{-1})}$

must satisfy

$\begin{eqnarray*} H^\prime(q_1) &\neq& 0 \\ \angle{H^\prime(q_1)}&=& 0 . \end{eqnarray*}$

Proof. We may repeat the above for .

Property. Every PR function has a causal inverse z transform .

Proof. This follows immediately from analyticity in the outer disk [342, pp. 30-36] However, we may give a more concrete proof as follows. Suppose is non-causal. Then there exists such that $h(-k)\neq 0$ . We have,

$\begin{eqnarray*} H(z)\;&\isdef & \sum_{n=-\infty }^\infty h(n)\,z^{-n}\\ &=& h(-k)z^k + \sum_{n \neq -k} h(n)\,z^{-n}. \end{eqnarray*}$

Hence, has at least one pole at infinity and cannot be PR by Property 3. Note that this pole at infinity cannot be cancelled since otherwise

$\begin{eqnarray*} && h(-k)z^k = \sum_{l\neq -k}\alpha(l)\,z^{-l}\\ &&\implies... ... \sum_{m\neq -k}\alpha(m)\,\delta(m-n)\\ &&\implies h(-k) = 0 \end{eqnarray*}$

which contradicts the hypothesis that is non-causal. $\Box$

Property. is PR iff it is analytic for $\left\vert z\right\vert>1$ , poles on the unit circle are simple with real and positive residues, and re $\left\{\\ right\}{H(e^{j\theta})\}\geq 0$ for $0\leq \theta\leq \pi$ .

Proof. If is positive real, the conditions stated hold by virtue of Property 3 and the definition of positive real.

To prove the converse, we first show nonnegativity on the upper semicircle implies nonnegativity over the entire circle.

$\begin{eqnarray*} \mbox{re}\left\{H(e^{j\theta})\right\} &\geq& 0, \qquad 0 \leq... ...eft\{H(e^{j\theta})\right\}&\geq& 0, \qquad -\pi<\theta\leq \pi. \end{eqnarray*}$

Alternatively, we might simply state that real implies re $\left\{H(e^{j\theta})\right\}$ is even in $\theta$ .

Next, since the function is analytic everywhere except at $z=\infty$ , it follows that $f(z)=e^{-H(z)}$ is analytic wherever is finite. There are no poles of outside the unit circle due to the analyticity assumption, and poles on the unit circle have real and positive residues. Referring again to the limiting form Eq. $\,$ (C.81) of near a pole on the unit circle at $e^{j\phi}$ , we see that, as ${\rho \to 0}$ , we have

$\displaystyle H(e^{j\phi} + \rho\,e^{j\psi})$	$\displaystyle \to$	$\displaystyle \frac{C}{ \rho}\, e^{j(\phi - \psi)}, \qquad \phi- \frac{\pi}{2} \leq \psi\leq \phi + \frac{\pi}{2}$
	$\displaystyle \isdef$	$\displaystyle \;\; \frac{C}{ \rho}\,e^{j\theta}, \qquad \theta \isdefs \phi - \psi, \quad - \frac{\pi}{ 2} \leq \theta \leq \frac{\pi}{ 2}$
$\displaystyle \implies \quad f(z)$	$\displaystyle \to$	$\displaystyle e^{- \frac{C}{ \rho}\, e^{j\theta}}$
	$\displaystyle =$	$\displaystyle e^{- \frac{C}{ \rho}\, \cos\theta} e^{-j \frac{C}{ \rho}\,\sin\theta}$
	$\displaystyle \to$	$\displaystyle 0, \protect$	(C.83)

since the residue

is positive, and the net angle $\theta$ does not exceed $\pm\pi/2$ . From Eq. $\,$ (C.83) we can state that for points $z,z^\prime$ with modulus $\geq 1$ , we have For all $\epsilon>0$ , there exists $\delta>0$ such that $\left\vert z-z^\prime \right\vert<\delta \;\implies\; \left\vert f(z)-f(z^\prime )\right\vert<\epsilon$ . Thus

is analytic in the strict outer disk, and continuous up to the unit circle which forms its boundary. By the maximum modulus theorem [83],

$\displaystyle \sup_{\left\vert z\right\vert\geq 1} \left\vert f(z)\right\vert \... ...)\right\}} = \inf_{\left\vert z\right\vert\geq 1} \mbox{re}\left\{H(z)\right\}$

occurs on the unit circle. Consequently,

$\displaystyle \inf_{-\pi<\theta\leq \pi}$ re $\displaystyle \left\{H(e^{j\theta})\right\}\geq 0 \quad\implies\quad\inf_{\left\vert z\right\vert\geq 1}$ re $\displaystyle \left\{H(z)\right\}\geq 0 \quad\implies\quad H(z)\;PR\,.$

For example, if a transfer function is known to be asymptotically stable, then a frequency response with nonnegative real part implies that the transfer function is positive real.

Note that consideration of leads to analogous necessary and sufficient conditions for to be positive real in terms of its zeros instead of poles. $\Box$

Relation to Stochastic Processes

Property. If a stationary random process $\{x_n\}$ has a rational power spectral density $Re^{j\omega}$ corresponding to an autocorrelation function $r(k)={\cal E}\left\{x_nx_{n+k}\right\}$ , then

$\displaystyle R_+(z)\isdef \frac{r(0)}{ 2} + \sum_{n=1}^\infty r(n)z^{-n}$

is positive real.

Proof.

By the representation theorem [19, pp. 98-103] there exists an asymptotically stable filter which will produce a realization of $\{x_n\}$ when driven by white noise, and we have $Re^{j\omega} = H(e^{j\omega})H(e^{-j\omega})$ . We define the analytic continuation of $Re^{j\omega}$ by $R(z) = H(z)H(z^{-1})$ . Decomposing into a sum of causal and anti-causal components gives

$\begin{eqnarray*} R(z) = \frac{b(z)b(z^{-1})}{ a(z)a(z^{-1})} &=&R_+(z) + R_-(z) \\ &=& \frac{q(z)}{ a(z)}+\frac{q(z^{-1})}{ a(z^{-1})} \end{eqnarray*}$

where is found by equating coefficients of like powers of in

$\displaystyle b(z)b(z^{-1})=q(z)a(z^{-1}) + a(z)q(z^{-1}).$

Since the poles of and are the same, it only remains to be shown that re $\left\{R_+(e^{j\omega})\right\}\geq 0,\;0\leq \omega\leq \pi$ .

Since spectral power is nonnegative, $Re^{j\omega}\geq 0$ for all $\omega$ , and so

$\begin{eqnarray*} Re^{j\omega}&\isdef & \sum_{n=-\infty }^\infty r(n)\,e^{j\omeg... ...)\\ &=&2\mbox{re}\left\{R_+(e^{j\omega})\right\}\\ &\geq& 0. \end{eqnarray*}$

$\Box$

Relation to Schur Functions

Definition. A Schur function is defined as a complex function analytic and of modulus not exceeding unity in $\vert z\vert\leq 1$ .

Property. The function

$\displaystyle S(z) \isdef \frac{1-R(z)}{ 1+R(z)}$

(C.84)

is a Schur function if and only if

is positive real.

Proof.

Suppose is positive real. Then for $\vert z\vert\geq 1$ , re $\left\{R(z)\right\}\geq 0 \,\,\Rightarrow\,\,1+$ re $\left\{R(z)\right\}\geq 0 \,\,\Rightarrow\,\,1+R(z)$ is PR. Consequently, is minimum phase which implies all roots of lie in the unit circle. Thus is analytic in $\vert z\vert\leq 1$ . Also,

$\displaystyle \left\vert S(e^{j\omega})\right\vert = \frac{1-2\mbox{re}\left\{R... ...left\{R(e^{j\omega})\right\} + \left\vert R(e^{j\omega})\right\vert^2} \leq 1.$

By the maximum modulus theorem,

takes on its maximum value in $\vert z\vert\geq 1$ on the boundary. Thus

is Schur.

Conversely, suppose is Schur. Solving Eq. $\,$ (C.84) for and taking the real part on the unit circle yields

$\begin{eqnarray*} R(z) &=& \alpha\,\frac{1-S(z)}{ 1+S(z)} \\ \mbox{re}\left\{R(... ...ert^2 }{ \left\vert 1+S(e^{j\omega})\right\vert^2}\\ &\geq& 0. \end{eqnarray*}$

If $S(z)=\alpha$ is constant, then $R(z)=(1-\vert\alpha\vert^2)/\vert 1+\alpha\vert^2$ is PR. If is not constant, then by the maximum principle, for $\vert z\vert>1$ . By Rouche's theorem applied on a circle of radius $1+\epsilon$ , $\epsilon>0$ , on which $\vert S(z)\vert<1$ , the function has the same number of zeros as the function in $\vert z\vert\geq 1+\epsilon$ . Hence, is minimum phase which implies is analytic for $z\geq 1$ . Thus is PR. $\Box$

Relation to functions positive real in the right-half plane

Property. re $\left\{H(z)\right\}\geq 0$ for $\left\vert z\right\vert\geq 1$ whenever

re $\displaystyle \left\{H\left(\frac{\alpha+s}{ \alpha-s}\right)\right\}\geq 0$

for re $\left\{s\right\}\geq 0$ , where $\alpha$ is any positive real number.

Proof. We shall show that the change of variable $z\leftarrow (\alpha+s)/(\alpha-s),\; \alpha>0$ , provides a conformal map from the z-plane to the s-plane that takes the region $\left\vert z\right\vert\geq 1$ to the region re $\left\{s\right\}\geq 0$ . The general formula for a bilinear conformal mapping of functions of a complex variable is given by

$\displaystyle \frac{(z-z_1)(z_2-z_3)}{ (z-z_3)(z_2-z_1)} = \frac{(s-s_1)(s_2-s_3)}{ (s-s_3)(s_2-s_1)}.$

(C.85)

In general, a bilinear transformation maps circles and lines into circles and lines [83]. We see that the choice of three specific points and their images determines the mapping for all and . We must have that the imaginary axis in the s-plane maps to the unit circle in the z-plane. That is, we may determine the mapping by three points of the form $z_i=e^{j\theta_i}$ and $s_i=j\omega_i,\;i=1,2,3$ . If we predispose one such mapping by choosing the pairs $(s_1=\pm\infty ) \leftrightarrow (z_1= -1)$ and $(s_3=0)\leftrightarrow(z_3=1)$ , then we are left with transformations of the form

$\displaystyle s=\left(s_2 \frac{z_2+1}{ z_2-1}\right)\left( \frac{z-1}{ z+1}\right) =\alpha\left( \frac{z-1}{ z+1}\right)$

$\displaystyle z\leftarrow \frac{\alpha+s}{ \alpha-s},$

(C.86)

Letting

be some point $j\omega$ on the imaginary axis, and

be some point $e^{j\theta}$ on the unit circle, we find that

$\displaystyle \alpha = j\omega \frac{e^{j\theta} + 1}{ e^{j\theta} -1 } = \omega \frac{\sin\theta}{ 1-\cos\theta} = \omega\cot(\theta/2)$

which gives us that $\alpha$ is real. To avoid degeneracy, we require $s_2\neq 0,\infty ,\;z_2\neq \pm1$ , and this translates to $\alpha$ finite and nonzero. Finally, to make the unit disk map to the left-half s-plane, $\omega$ and $\theta$ must have the same sign in which case $\alpha>0$ . $\Box$

There is a bonus associated with the restriction that $\alpha$ be real which is that

$\displaystyle z= \frac{\alpha+s}{ \alpha-s}\in\Re \quad \leftrightarrow \quad s = \alpha \frac{z-1}{ z+1}\in\Re .$

(C.87)

We have therefore proven

Property.

$\displaystyle H(z)$ PR $\displaystyle \;\Leftrightarrow\;H\left(\frac{\alpha+s}{\alpha-s}\right)$ PR $\displaystyle , \protect$

where $\alpha$ is any positive real number.

The class of mappings of the form Eq. $\,$ (C.85) which take the exterior of the unit circle to the right-half plane is larger than the class Eq. $\,$ (C.86). For example, we may precede the transformation Eq. $\,$ (C.86) by any conformal map which takes the unit disk to the unit disk, and these mappings have the algebraic form of a first order complex allpass whose zero lies inside the unit circle.

$\displaystyle z\leftarrow e^{j\theta} \frac{w-w_0}{ \overline{w_0}\,w - 1},\qquad \left\vert w_0\right\vert<1$

(C.88)

where

is the zero of the allpass and the image (also pre-image) of the origin, and $\theta$ is an angle of pure rotation. Note that Eq. $\,$ (C.88) is equivalent to a pure rotation, followed by a real allpass substitution (

real), followed by a pure rotation. The general preservation of condition (2) in Def. 2 forces the real axis to map to the real axis. Thus rotations by other than $\pi$ are useless, except perhaps in some special cases. However, we may precede Eq. $\,$ (C.86) by the first order real allpass substitution

$\displaystyle z\leftarrow \frac{w-r}{ r\,w - 1},\qquad \left\vert r\right\vert<1,\; r$ real $\displaystyle , \protect$

which maps the real axis to the real axis. This leads only to the composite transformation,

$\displaystyle z \leftarrow \frac{s+ \left(\alpha \frac{1-r}{1+r}\right)}{ s- \left(\alpha \frac{1-r}{1+r}\right)}$

which is of the form Eq. $\,$ (C.86) up to a minus sign (rotation by $\pi$ ). By inspection of Eq. $\,$ (C.85), it is clear that sign negation corresponds to the swapping of points

and

, or

and

. Thus the only extension we have found by means of the general disk to disk pre-transform, is the ability to interchange two of the three points already tried. Consequently, we conclude that the largest class of bilinear transforms which convert functions positive real in the outer disk to functions positive real in the right-half plane is characterized by

$\displaystyle z \leftarrow \pm \frac{\alpha+s}{ \alpha-s}.$

(C.89)

Riemann's theorem may be used to show that Eq. $\,$ (C.89) is also the largest such class of conformal mappings. It is not essential, however, to restrict attention solely to conformal maps. The pre-transform $z\leftarrow \overline{z}$ , for example, is not conformal and yet PR is preserved.

The bilinear transform is one which is used to map analog filters into digital filters. Another such mapping is called the matched transform [362]. It also preserves the positive real property.

Property. is PR if $H(e^{sT})$ is positive real in the analog sense, where is interpreted as the sampling period.

Proof. The mapping $z \leftarrow e^{sT}$ takes the right-half -plane to the outer disk in the -plane. Also is real if is real. Hence $H(e^{sT})$ PR implies PR. (Note, however, that rational functions do not in general map to rational functions.) $\Box$

These transformations allow application of the large battery of tests which exist for functions positive real in the right-half plane [524].

Special cases and examples

The sum of positive real functions is positive real.
The difference of positive real functions is conditionally positive real.
The product or division of positive real functions is conditionally PR.
PR $\implies \alpha\,z^{\pm k}H(z)$ not PR for $\alpha>0,k\geq 2$ .

Minimum Phase (MP) polynomials in

All properties of MP polynomials apply without modification to marginally stable allpole transfer functions (cf. Property 2):

Every first-order MP polynomial is positive real.
Every first-order MP polynomial $b(z)=1+b_1\,z^{-1}$ is such that $\frac{1}{ b(z)} - \frac{1}{ 2}$ is positive real.
A PR second-order MP polynomial with complex-conjugate zeros,

$\begin{eqnarray*} H(z)&=& 1+b_1z^{-1}+b_2z^{-2}\\ &=& 1-(2R\cos\phi)z^{-1}+R^2z^{-2},\quad R\leq 1 \end{eqnarray*}$

satisfies

$\displaystyle R^2 + \frac{\cos^2\phi}{ 2} \leq 1.$
If $2R^2+\cos^2\phi=2$ , then re $\left\{H(e^{j\omega})\right\}$ has a double zero at

$\begin{eqnarray*} \omega &=& \cos^{-1}\left(\pm \sqrt{ \frac{1-R^2}{ 2R^2}}\righ... ...frac{1}{ \sqrt{2}} \frac{\cos\phi}{ \sqrt{1+\sin^2\phi}}\right). \end{eqnarray*}$
All polynomials of the form

$\displaystyle H(z)=1+R^nz^{-n},\quad R\leq 1$
are positive real. (These have zeros uniformly distributed on a circle of radius .)

Miscellaneous Properties

If all poles and zeros of a PR function are on the unit circle, then they alternate along the circle. Since this property is preserved by the bilinear transform, it is true in both the and planes. It can be viewed as a consequence of the $\pm 90\degrees$ phase bounds for positive-real functions.
If is PR, then so is $B^\prime(s)/A^\prime(s)$ , where the prime denotes differentiation in .

Loaded Waveguide Junctions

In this section, scattering relations will be derived for the general case of N waveguides meeting at a load. When a load is present, the scattering is no longer lossless, unless the load itself is lossless. (i.e., its impedance has a zero real part). For , $v^{+}_i$ will denote a velocity wave traveling into the junction, and will be called an ``incoming'' velocity wave as opposed to ``right-going.''^C.9

**Figure C.29:** Four ideal strings intersecting at a point to which a lumped impedance is attached. This is a series junction for transverse waves.
$\includegraphics[width=\twidth]{eps/fNstrings}$

Consider first the series junction of waveguides containing transverse force and velocity waves. At a series junction, there is a common velocity while the forces sum. For definiteness, we may think of ideal strings intersecting at a single point, and the intersection point can be attached to a lumped load impedance , as depicted in Fig.C.29 for . The presence of the lumped load means we need to look at the wave variables in the frequency domain, i.e., $V(s) = {\cal L}\{v\}$ for velocity waves and $F(s) = {\cal L}\{f\}$ for force waves, where ${\cal L}\{\cdot\}$ denotes the Laplace transform. In the discrete-time case, we use the transform instead, but otherwise the story is identical. The physical constraints at the junction are

$\displaystyle V_1(s) = V_2(s) = \cdots = V_N(s)$	$\displaystyle \isdef$	$\displaystyle V_J(s)$	(C.90)
$\displaystyle F_1(s) + F_2(s) + \cdots + F_N(s)$	$\displaystyle =$	$\displaystyle V_J(s) R_J(s) \isdef F_J(s)$	(C.91)

where the reference direction for the load force

is taken to be opposite that for the

. (It can be considered the ``equal and opposite reaction'' force at the junction.) For a wave traveling into the junction, force is positive pulling up, acting toward the junction. When the load impedance

is zero, giving a free intersection point, the junction reduces to the unloaded case, and signal scattering will be energy preserving. In general, the loaded junction is lossless if and only if re $\left\{R_J(j\omega)\right\}\equiv0$ , and it is memoryless if and only if im $\left\{R_J(j\omega)\right\}\equiv0$ .

The parallel junction is characterized by

$\displaystyle F_1(s) = F_2(s) = \cdots = F_N(s)$	$\displaystyle \isdef$	$\displaystyle F_J(s)$	(C.92)
$\displaystyle V_1(s) + V_2(s) + \cdots + V_N(s)$	$\displaystyle =$	$\displaystyle F_J(s)/R_J(s) \isdef V_J(s)$	(C.93)

For example,

could be pressure in an acoustic tube and

the corresponding volume velocity. In the parallel case, the junction reduces to the unloaded case when the load impedance

goes to infinity.

The scattering relations for the series junction are derived as follows, dropping the common argument `' for simplicity:

$\displaystyle R_J V_J$	$\displaystyle =$	$\displaystyle \sum_{i=1}^N F_i = \sum_{i=1}^N (F^+_i + F^-_i)$	(C.94)
	$\displaystyle =$	$\displaystyle \sum_{i=1}^N (R_i V^+_i - R_i \underbrace{V^-_i}_{V_J-V^+_i})$	(C.95)
	$\displaystyle =$	$\displaystyle \sum_{i=1}^N (2 R_i V^+_i - R_i V_J)$	(C.96)

where

is the wave impedance in the

th waveguide, a real, positive constant. Bringing all terms containing

to the left-hand side, and solving for the junction velocity gives

$\displaystyle V_J$	$\displaystyle =$	$\displaystyle 2\left(R_J + \sum_{i=1}^N R_i\right)^{-1} \sum_{i=1}^N R_i V^+_i$	(C.97)
	$\displaystyle \isdef$	$\displaystyle \sum_{i=1}^N{\cal A}_i(s) V^+_i(s)$	(C.98)

(written to be valid also in the multivariable case involving square impedance matrices

[433]), where

$\displaystyle {\cal A}_i(s) \isdef \frac{2R_i}{R_J(s) + R_1 + \cdots + R_N}$

(C.99)

Finally, from the basic relation

, the outgoing velocity waves can be computed from the junction velocity and incoming velocity waves as

$\displaystyle V^-_i(s) = V_J(s) - V^+_i(s)$

(C.100)

Similarly, the scattering relations for the loaded parallel junction are given by

$\displaystyle F_J(s)$	$\displaystyle =$	$\displaystyle \sum_{i=1}^N{\cal A}_i(s) F^+_i(s), \quad {\cal A}_i(s) \isdef \frac{2\Gamma _i}{\Gamma _J(s) + \Gamma _1 + \cdots + \Gamma _N}$
$\displaystyle F^-_i(s)$	$\displaystyle =$	$\displaystyle F_J(s) - F^+_i(s)$	(C.101)

where

is the Laplace transform of the force across all elements at the junction, $\Gamma _J(s)$ is the load admittance, and $\Gamma _i=1/R_i$ are the branch admittances.

It is interesting to note that the junction load is equivalent to an st waveguide having a (generalized) wave impedance given by the load impedance. This makes sense when one recalls that a transmission line can be ``perfectly terminated'' (i.e., suppressing all reflections from the termination) using a lumped resistor equal in value to the wave impedance of the transmission line. Thus, as far as a traveling wave is concerned, there is no difference between a wave impedance and a lumped impedance of the same value.

In the unloaded case, , and we can return to the time domain and define (for the series junction)

$\displaystyle \alpha_i = \frac{2R_i}{R_1 + \cdots + R_N}$

(C.102)

These we call the alpha parameters, and they are analogous to those used to characterize ``adaptors'' in wave digital filters (§F.2.2). For unloaded junctions, the alpha parameters obey

$\displaystyle 0\leq\alpha_i \leq 2$

(C.103)

and

$\displaystyle \sum_{i=1}^N\alpha_i = 2$

(C.104)

In the unloaded case, the series junction scattering relations are given (in the time domain) by

$\displaystyle v_J(t)$	$\displaystyle =$	$\displaystyle \sum_{i=1}^N \alpha_i v^+_i(t) \protect$	(C.105)
$\displaystyle v^-_i(t)$	$\displaystyle =$	$\displaystyle v_J(t) - v^+_i(t)$	(C.106)

The alpha parameters provide an interesting and useful parametrization of waveguide junctions. They are explicitly the coefficients of the incoming traveling waves needed to compute junction velocity for a series junction (or junction force or pressure at a parallel junction), and losslessness is assured provided only that the alpha parameters be nonnegative and sum to

. Having them sum to something less than

simulates a ``resistive load'' at the junction.

Note that in the lossless, equal-impedance case, in which all waveguide impedances have the same value , (C.102) reduces to

$\displaystyle \alpha_i = \frac{2}{N}$

(C.107)

When, furthermore,

is a power of two, we have that there are no multiplies in the scattering relations (C.105). This fact has been used to build multiply-free reverberators and other structures using digital waveguide meshes [430,518,396,520].

An elaborated discussion of strings intersection at a load is given in in §9.2.1. Further discussion of the digital waveguide mesh appears in §C.14.

Two Coupled Strings

**Figure C.30:** Two strings terminated at a common bridge impedance.
$\includegraphics{eps/f_two_strings}$

Two Ideal Strings Coupled at an Impedance

A diagram for the two-string case is shown in Fig. C.30. This situation is a special case of the loaded waveguide junction, Eq. $\,$ (C.98), with the number of waveguides being , and the junction load being the transverse driving-point impedance where the string drives the bridge. If the bridge is passive, then its impedance is a positive real function (see §C.11.2). For a direct derivation, we need only observe that (1) the string velocities of each string endpoint must each be equal to the velocity of the bridge, or , and (2) the sum of forces of both strings equals the force applied to the bridge: . The bridge impedance relates the force and velocity of the bridge via . Expanding into traveling wave components in the Laplace domain, we have

$\displaystyle R_b(s) V_b(s)$	$\displaystyle =$	$\displaystyle F_b(s) = F_1(s) + F_2(s)$	(C.108)
	$\displaystyle =$	$\displaystyle [F^+_1(s) + F^-_1(s)] + [F^+_2(s) + F^-_2(s)]$	(C.109)
	$\displaystyle =$	$\displaystyle R_1 \{V^+_1(s) - [V_b(s) - V^+_1(s)] \}$	(C.110)
	$\displaystyle \,+\,$	$\displaystyle R_2 \{V^+_2(s) - [V_b(s) - V^+_2(s)]\}$	(C.111)

$\displaystyle V_b(s) = H_b(s) [ R_1 V^+_1(s) + R_2 V^+_2(s) ]$

where

is the wave impedance of string

, and

$\displaystyle H_b(s)\isdef \frac{2}{R_b(s) + R_1 + R_2}$

Thus, in the time domain, the incoming velocity waves are scaled by their respective wave impedances, summed together, and filtered according to the transfer function

to obtain the velocity of the bridge

Given the filter output , the outgoing traveling velocity waves are given by

$\displaystyle v^-_1(t)$	$\displaystyle =$	$\displaystyle v_b(t) - v^+_1(t)$	(C.112)
$\displaystyle v^-_2(t)$	$\displaystyle =$	$\displaystyle v_b(t) - v^+_2(t) \$	(C.113)

Thus, the incoming waves are subtracted from the bridge velocity to get the outgoing waves.

Since when , and vice versa exchanging strings and , may be interpreted as the transmission admittance filter associated with the bridge coupling. It can also be interpreted as the bridge admittance transfer function from every string, since its output is the bridge velocity resulting from the sum of incident traveling force waves.

A general coupling matrix contains a filter transfer function in each entry of the matrix. For strings, each conveying a single type of wave (e.g., horizontally polarized), the general linear coupling matrix would have transfer-function entries. In the present formulation, only one transmission filter is needed, and it is shared by all the strings meeting at the bridge. It is easy to show that the shared transmission filter for two coupled strings generalizes to strings coupled at a common bridge impedance: From (C.98), we have

$\displaystyle V_b(s) = H_b(s) \sum_{i=1}^N R_i V^{+}_i(s)$

where

$\displaystyle H_b(s) = \frac{2}{R_b(s) + \sum_{i=1}^N R_i}$

Thus,

is the shared portion of the bridge filtering, leaving only a scaling according to relative impedance to be done in each branch.

The above sequence of operations is formally similar to the one multiply scattering junction frequently used in digital lattice filters [297]. In this context, it would be better termed the ``one-filter scattering termination.''

When the two strings are identical (as would be appropriate in a model for coupled piano strings), the computation of bridge velocity simplifies to

$\displaystyle V_b(s) = H_b(s) [V^+_1(s) + V^+_2(s)] \protect$

(C.114)

where $H_b(s)\isdef 2/[2 + R_b(s)/R]$ is the velocity transmission filter. In this case, the incoming velocities are simply summed and fed to the transmission filter which produces the bridge velocity at its output. A commuted simulation diagram appears in Fig. C.31.

**Figure C.31:** General linear coupling of two equal-impedance strings using a common bridge filter.
$\includegraphics{eps/fcouplednompy}$

Note that a yielding bridge introduces losses into all attached strings. Therefore, in a maximally simplified implementation, all string loop filters (labeled LPF and LPF in Fig.C.31) may be eliminated, resulting in only one filter--the transmission filter--serving to provide all losses in a coupled-string simulation. If that transmission filter has no multiplies, then neither does the entire multi-string simulation.

Coupled Strings Eigenanalysis

In §6.12.2, general coupling of horizontal and vertical planes of vibration in an ideal string was considered. This eigenanalysis will now be applied here to obtain formulas for the damping and mode tuning caused by the coupling of two identical strings at a bridge. This is the case that arises in pianos [543].

The general formula for linear, time-invariant coupling of two strings can be written, in the frequency domain, as

$\displaystyle \left[\begin{array}{c} V_1^-(s) \\ [2pt] V_2^-(s) \end{array}\rig... ... \left[\begin{array}{c} V_1^+(s) \\ [2pt] V_2^+(s) \end{array}\right]. \protect$

(C.115)

Filling in the elements of this coupling matrix $\mathbf{H}_c$ from the results of §C.13.1, we obtain

$\displaystyle \mathbf{H}_c(s) = \left[\begin{array}{cc} 1-H_b(s) & -H_b(s) \\ [2pt] -H_b(s) & 1-H_b(s) \end{array}\right]$

where

$\displaystyle H_b(s) = \frac{2}{2+\tilde{R}_b}.$

Here $\tilde{R}_b\isdef R_b/R$ is the bridge impedance divided by the string impedance. Treating $\mathbf{H}_c(s)$ as a constant complex matrix for each fixed

, the eigenvectors are found^C.10to be

$\displaystyle \underline{e}_1 = \left[\begin{array}{c} 1 \\ [2pt] 1 \end{array}... ...uad \underline{e}_2 = \left[\begin{array}{c} 1 \\ [2pt] -1 \end{array}\right],$

respectively, and the eigenvalues are

$\displaystyle \lambda_1(s) = 1 - 2H_b(s), \qquad \lambda_2 = 1.$

Note that only one eigenvalue depends on $s=j\omega$ , and neither eigenvector is a function of

We conclude that ``in-phase vibrations'' see a longer effective string length, lengthened by the phase delay of

$\displaystyle 1-2H_b = \frac{\tilde{R}_b(s)-2}{\tilde{R}_b(s)+2} = \frac{R_b(s)-2R}{R_b(s)+2R}$

which is the reflectance seen from two in-phase strings each having impedance

. This makes physical sense because the in-phase vibrations will move the bridge in the vertical direction, causing more rapid decay of the in-phase mode.

We similarly conclude that the ``anti-phase vibrations'' see no length correction at all, because the bridge point does not move at all in this case. In other words, any bridge termination at a point is rigid with respect to anti-phase vibration of the two strings connected to that point.

The above analysis predicts that, in ``stiffness controlled'' frequency intervals (in which the bridge ``looks like a damped spring''), the ``initial fast decay'' of a piano note should be a measurably flatter than the ``aftersound'' which should be exactly in tune as if the termination were rigid.

Digital Waveguide Mesh

In §C.12, the theory of multiport scattering was derived, i.e., the reflections and transmissions that occur when digital waveguides having wave impedances are connected together. It was noted that when is a power of two, there are no multiplies in the scattering relations Eq. $\,$ (C.105), and that this fact has been used to build multiply-free reverberators and other structures using digital waveguide meshes [430,518,146,396,520,521,398,399,401,55,202,321,320,322,422,33].

The Rectilinear 2D Mesh

**Figure C.32:** The 2D rectilinear digital waveguide mesh.
$\includegraphics[width=4in]{eps/SchematicWaveguideMesh}$

Figure C.32 shows the basic layout of the rectilinear 2D waveguide mesh. It can be thought of as simulating a plane using 1D digital waveguides in the same way that a tennis racket acts as a membrane composed of 1D strings.

At each node (string intersection), we have the following simple formula for the node velocity in terms of the four incoming traveling-wave components:

$\displaystyle v = \frac{\mbox{in}_{1} + \mbox{in}_{2} + \mbox{in}_{3} + \mbox{in}_{4}}{2}$

By continuity of velocity in a series connection (all string endpoints must move together at the node), the outgoing velocity-wave components must be given by

$\displaystyle \hbox{out}_k = v -$ in $\displaystyle _{k}, \qquad k=1,2,3,4.$

Dispersion

Since the digital waveguide mesh is lossless by construction (when modeling lossless membranes and volumes), and since it is also linear and time-invariant by construction, being made of ordinary digital filtering computations, there is only one type of error exhibited by the mesh: dispersion. Dispersion can be quantified as an error in propagation speed as a function of frequency and direction along the mesh. The mesh geometry (rectilinear, triangular, hexagonal, tetrahedral, etc.) strongly influences the dispersion properties. Many cases are analyzed in [55] using von Neumann analysis (see also Appendix D).

The triangular waveguide mesh [146] turns out to be the simplest mesh geometry in 2D having the least dispersion variation as a function of direction of propagation on the mesh. In other terms, the triangular mesh is closer to isotropic than all other known elementary geometries. The interpolated waveguide mesh [398] can also be configured to optimize isotropy, but at a somewhat higher compuational cost.

Recent Developments

An interesting approach to dispersion compensation is based on frequency-warping the signals going into the mesh [399]. Frequency warping can be used to compensate frequency-dependent dispersion, but it does not address angle-dependent dispersion. Therefore, frequency-warping is used in conjunction with an isotropic mesh.

The 3D waveguide mesh [518,521,399] is seeing more use for efficient simulation of acoustic spaces [396,182]. It has also been applied to statistical modeling of violin body resonators in [203,202,422,428], in which the digital waveguide mesh was used to efficiently model only the ``reverberant'' aspects of a violin body's impulse response in statistically matched fashion (but close to perceptually equivalent). The ``instantaneous'' filtering by the violin body is therefore modeled using a separate equalizer capturing the important low-frequency body and air modes explicitly. A unified view of the digital waveguide mesh and wave digital filters (§F.1) as particular classes of energy invariant finite difference schemes (Appendix D) appears in [54]. The problem of modeling diffusion at a mesh boundary was addressed in [268], and maximally diffusing boundaries, using quadratic residue sequences, was investigated in [279]; an introduction to this topic is given in §C.14.6 below.

2D Mesh and the Wave Equation

**Figure C.33:** Region of nodes in a rectilinear waveguide mesh.
$\includegraphics{eps/mesh}$

**Figure C.34:** Zoom-in about node at time in a rectilinear waveguide mesh, showing traveling-wave components entering and leaving the node. (All variables are at time ,)
$\includegraphics{eps/meshzoom}$

Consider the 2D rectilinear mesh, with nodes at positions and , where and are integers, and and denote the spatial sampling intervals along and , respectively (see Fig.C.33). Then from Eq. $\,$ (C.105) the junction velocity $v_{lm}$ at time is given by

$\displaystyle v_{lm}(n) = \frac{1}{2}\left[ v_{lm}^{+\textsc{n}}(n) + v_{lm}^{+\textsc{e}}(n) + v_{lm}^{+\textsc{s}}(n) + v_{lm}^{+\textsc{w}}(n)\right]$

where $v_{lm}^{+\textsc{n}}(n)$ is the ``incoming wave from the north'' to node

, and similarly for the waves coming from east, west, and south (see Fig.C.34).

These incoming traveling-wave components arrive from the four neighboring nodes after a one-sample propagation delay. For example, $v_{lm}^{+\textsc{n}}(n)$ , arriving from the north, departed from node at time , as $v_{l,m+1}^{-\textsc{s}}(n-1)$ . Furthermore, the outgoing components at time will arrive at the neighboring nodes one sample in the future at time . For example, $v_{lm}^{-\textsc{n}}(n)$ will become $v_{l,m+1}^{+\textsc{s}}(n+1)$ . Using these relations, we can write $v_{lm}(n+1)$ in terms of the four outgoing waves from its neighbors at time :

$\displaystyle v_{lm}(n+1) = \frac{1}{2}\left[ v_{l,m+1}^{-\textsc{s}}(n) + v_{l... ...}}(n) + v_{l,m-1}^{-\textsc{n}}(n) + v_{l-1,m}^{-\textsc{e}}(n)\right] \protect$

(C.116)

where, for instance, $v_{lm}^{-\textsc{n}}(n)$ is the ``outgoing wave to the north'' from node

. Similarly, the outgoing waves leaving $v_{lm}(n-1)$ become the incoming traveling-wave components of its neighbors at time

$\displaystyle v_{lm}(n-1) = \frac{1}{2}\left[ v_{l,m+1}^{+\textsc{s}}(n) + v_{l... ...}}(n) + v_{l,m-1}^{+\textsc{n}}(n) + v_{l-1,m}^{+\textsc{e}}(n)\right] \protect$

(C.117)

This may be shown in detail by writing

$\begin{eqnarray*} v_{lm}(n-1) &=& \frac{1}{2}[v_{lm}^{+\textsc{n}}(n-1) + \cdot... ...}^{-\textsc{n}}(n-1) + \cdots + v_{lm}^{-\textsc{w}}(n-1)\right] \end{eqnarray*}$

so that

$\begin{eqnarray*} v_{lm}(n-1) &=& \frac{1}{2}[v_{lm}^{-\textsc{n}}(n-1) + \cdot... ... v_{l,m-1}^{+\textsc{n}}(n) + v_{l-1,m}^{+\textsc{e}}(n)\right]. \end{eqnarray*}$

Adding Equations (C.116-C.116), replacing terms such as $v_{l,m+1}^{+\textsc{s}}(n) + v_{l,m+1}^{-\textsc{s}}(n)$ with $v_{l,m+1}(n)$ , yields a computation in terms of physical node velocities:

$\begin{eqnarray*} \lefteqn{v_{lm}(n+1) + v_{lm}(n-1) = } \\ & & \frac{1}{2}\left[ v_{l,m+1}(n) + v_{l+1,m}(n) + v_{l,m-1}(n) + v_{l-1,m}(n)\right] \end{eqnarray*}$

Thus, the rectangular waveguide mesh satisfies this equation giving a formula for the velocity at node , in terms of the velocity at its neighboring nodes one sample earlier, and itself two samples earlier. Subtracting $2v_{lm}(n)$ from both sides yields

$\begin{eqnarray*} \lefteqn{v_{lm}(n+1) - 2 v_{lm}(n) + v_{lm}(n-1)} \\ &=& \fra... .... \left[v_{l+1,m}(n) - 2 v_{lm}(n) + v_{l-1,m}(n)\right]\right\} \end{eqnarray*}$

Dividing by the respective sampling intervals, and assuming (square mesh-holes), we obtain

$\begin{eqnarray*} \lefteqn{\frac{v_{lm}(n+1) - 2 v_{lm}(n) + v_{lm}(n-1)}{T^2}} ... ...ft.\frac{v_{l+1,m}(n) - 2 v_{lm}(n) + v_{l-1,m}(n)}{X^2}\right]. \end{eqnarray*}$

In the limit, as the sampling intervals approach zero such that remains constant, we recognize these expressions as the definitions of the partial derivatives with respect to , , and , respectively, yielding

$\displaystyle \frac{\partial^2 v(t,x,y)}{\partial t^2} = \frac{X^2}{2T^2} \left... ...^2 v(t,x,y)}{\partial x^2} + \frac{\partial^2 v(t,x,y)}{\partial y^2} \right].$

This final result is the ideal 2D wave equation $\ddot v = c^2 \nabla^2 v$ , i.e.,

$\displaystyle \frac{\partial^2 v}{\partial t^2} = c^2 \left[ \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right]$

with

$\displaystyle c = \frac{1}{\sqrt{2}}\frac{X}{T} = \frac{\sqrt{2}}{2}\frac{X}{T}. \protect$

(C.118)

Discussion regarding solving the 2D wave equation subject to boundary conditions appears in §B.8.3. Interpreting this value for the wave propagation speed

, we see that every two time steps of

seconds corresponds to a spatial step of $\sqrt{2}X$ meters. This is the distance from one diagonal to the next in the square-hole mesh. We will show later that diagonal directions on the mesh support exact propagation (of plane waves traveling at 45-degree angles with respect to the

axes). In the

and

directions, propagation is highly dispersive, meaning that different frequencies travel at different speeds. The exactness of 45-degree angles can be appreciated by considering Huygens' principle on the mesh.

The Lossy 2D Mesh

Because the finite-difference form of the digital waveguide mesh is the more efficient computationally than explicitly computing scattering wave variables (too see this, count the multiplies required per node), it is of interest to consider the finite-difference form also in the case of frequency-dependent losses. The method of §C.5.5 extends also to the waveguide mesh, which can be shown by generalizing the results of §C.14.4 above using the technique of §C.5.5.

The basic idea is once again that wave propagation during one sampling interval (in time) is associated with linear filtering by . That is, is regarded as the per-sample wave propagation filter.

Diffuse Reflections in the Waveguide Mesh

In [416], Manfred Schroeder proposed the design of a diffuse reflector based on a quadratic residue sequence. A quadratic residue sequence corresponding to a prime number is the sequence mod , for all integers . The sequence is periodic with period , so it is determined by for $n=0,1,2,\ldots,p-1$ (i.e., one period of the infinite sequence).

For example, when , the first period of the quadratic residue sequence is given by

$\begin{eqnarray*} a_7 &=& [0^2,1^2,2^2,3^2,4^2,5^2,6^2] \quad (\mbox{mod }7)\\ &=& [0, 1, 4, 2, 2, 4, 1] \end{eqnarray*}$

An amazing property of these sequences is that their Fourier transforms have precisely constant magnitudes. That is, the sequence

$\displaystyle c_p(n) \isdef e^{j\frac{2\pi}{p} a_p(n)}$

has a DFT with exactly constant magnitude:^C.11

$\displaystyle \vert C_p(\omega_k)\vert \isdef \vert\dft _k(c_p)\vert \isdef \l... ...^{p-1} c_p(n) e^{-j2\pi nk/p}\right\vert = \sqrt{p}, \quad \forall k\in[0,p-1]$

This property can be used to give highly diffuse reflections for incident plane waves.

Figure C.35 presents a simple matlab script which demonstrates the constant-magnitude Fourier property for all odd integers from 1 to 99.

Figure C.35: Matlab script for demonstrating the Fourier property of an odd-length quadratic residue sequence.

function [c] = qrsfp(Ns)
%QRSFP Quadratic Residue Sequence Fourier Property demo
  if (nargin<1)
     Ns = 1:2:99; % Test all odd integers from 1 to 99
  end
  for N=Ns
    a = mod([0:N-1].^2,N);
    c = zeros(N-1,N);
    CM = zeros(N-1,N);
    c = exp(j*2*pi*a/N);
    CM = abs(fft(c))*sqrt(1/N);
    if (abs(max(CM)-1)>1E-10) || (abs(min(CM)-1)>1E-10)
       warn(sprintf("Failure for N=%d",N));
    end
  end
  r = exp(2i*pi*[0:100]/100); % a circle
  plot(real(r), imag(r),"k"); hold on;
  plot(c,"-*k"); % plot sequence in complex plane
end

Quadratic residue diffusers have been applied as boundaries of a 2D digital waveguide mesh in [279]. An article reviewing the history of room acoustic diffusers may be found in [94].

FDNs as Digital Waveguide Networks

This section supplements §2.7 on Feedback Delay Networks in the context of digital waveguide theory. Specifically, we review the interpretation of an FDN as a special case of a digital waveguide network, summarizing [463,464,385].

Figure C.36 illustrates an -branch DWN. It consists of a single scattering junction, indicated by a white circle, to which branches are connected. The far end of each branch is terminated by an ideal non-inverting reflection (black circle). The waves traveling into the junction are associated with the FDN delay line outputs , and the length of each waveguide is half the length of the corresponding FDN delay line (since a traveling wave must traverse the branch twice to complete a round trip from the junction to the termination and back). When is odd, we may replace the reflecting termination by a unit-sample delay.

**Figure C.36:** Waveguide network consisting of a single scattering junction, indicated by an open circle, to which branches are connected. The far end of each branch is terminated by an ideal, non-inverting reflection.
$\includegraphics[scale=0.5]{eps/DWN}$

Lossless Scattering

The delay-line inputs (outgoing traveling waves) are computed by multiplying the delay-line outputs (incoming traveling waves) by the $N\times N$ feedback matrix (scattering matrix) $\mathbf{A}= [a_{i,j}]$ . By defining , , we obtain the more usual DWN notation

$\displaystyle \mathbf{p}^- = \mathbf{A}\mathbf{p}^+$

(C.119)

where $\mathbf{p}^+$ is the vector of incoming traveling-wave samples arriving at the junction at time

, $\mathbf{p}^-$ is the vector of outgoing traveling-wave samples leaving the junction at time

, and $\mathbf{A}$ is the scattering matrix associated with the waveguide junction.

The junction of physical waveguides determines the structure of the matrix $\mathbf{A}$ according to the basic principles of physics.

Considering the parallel junction of lossless acoustic tubes, each having characteristic admittance $\Gamma_j=1/R_j$ , the continuity of pressure and conservation of volume velocity at the junction give us the following scattering matrix for the pressure waves [433]:

$\displaystyle {\bf A} = \left[ \begin{array}{rrrr} \frac{2 \Gamma_{1}}{\Gamma_J... ...{2}}{\Gamma_J} & \dots & \frac{2 \Gamma_{N}}{\Gamma_J} -1\\ \end{array} \right]$

(C.120)

where

$\displaystyle \Gamma_J = \sum_{i=1}^N\Gamma_{i}.$

(C.121)

Equation (C.121) can be derived by first writing the volume velocity at the

-th tube in terms of pressure waves as $v_j = (p_j^+ - p_j^-)\Gamma_j$ . Applying the conservation of velocity we can find the expression

$\displaystyle p = 2 \sum_{i=1}^{N}\Gamma_{i} p_i^+ / \Gamma_J$

for the junction pressure. Finally, if we express the junction pressure as the sum of incoming and outgoing pressure waves at any branch, we derive (C.121). See §C.12 for further details.

Normalized Scattering

For ideal numerical scaling in the sense, we may choose to propagate normalized waves which lead to normalized scattering junctions analogous to those encountered in normalized ladder filters [297]. Normalized waves may be either normalized pressure $\tilde{p}_j^+ = p_j^+\sqrt{\Gamma_i}$ or normalized velocity $\tilde{v}_j^+ = v_j^+/\sqrt{\Gamma_i}$ . Since the signal power associated with a traveling wave is simply ${\cal P_j^+} = (\tilde{p}_j^+)^2 = (\tilde{v}_j^+)^2$ , they may also be called root-power waves [432]. Appendix C develops this topic in more detail.

The scattering matrix for normalized pressure waves is given by

$\displaystyle \tilde{\mathbf{A}}= \left[ \begin{array}{llll} \frac{2 \Gamma_{1}... ..._{2}}}{\Gamma_J} & \dots & \frac{2 \Gamma_{n}}{\Gamma_J} -1 \end{array} \right]$

(C.122)

The normalized scattering matrix can be expressed as a negative Householder reflection

$\displaystyle \tilde{\mathbf{A}}= \frac{2}{ \vert\vert\,\tilde{{\bm \Gamma}}\,\vert\vert ^2}\tilde{{\bm \Gamma}}\tilde{{\bm \Gamma}}^T-\mathbf{I}$

(C.123)

where $\tilde{{\bm \Gamma}}^T= [\sqrt{\Gamma_1},\ldots,\sqrt{\Gamma_N}]$ , and $\Gamma_i$ is the wave admittance in the

th waveguide branch. To eliminate the sign inversion, the reflections at the far end of each waveguide can be chosen as -1 instead of 1. The geometric interpretation of (C.124) is that the incoming pressure waves are reflected about the vector $\tilde{{\bm \Gamma}}$ . Unnormalized scattering junctions can be expressed in the form of an ``oblique'' Householder reflection $\mathbf{A}= 2\mathbf{1}{\bm \Gamma}^T/\left<\mathbf{1},{{\bm \Gamma}}\right>-\mathbf{I}$ , where $\mathbf{1}^T=[1,\ldots,1]$ and ${\bm \Gamma}^T= [\Gamma_1,\ldots,\Gamma_N]$ .

General Conditions for Losslessness

The scattering matrices for lossless physical waveguide junctions give an apparently unexplored class of lossless FDN prototypes. However, this is just a subset of all possible lossless feedback matrices. We are therefore interested in the most general conditions for losslessness of an FDN feedback matrix. The results below are adapted from [463,385].

Consider the general case in which $\mathbf{A}$ is allowed to be any scattering matrix, i.e., it is associated with a not-necessarily-physical junction of physical waveguides. Following the definition of losslessness in classical network theory, we may say that a waveguide scattering matrix $\mathbf{A}$ is said to be lossless if the total complex power [35] at the junction is scattering invariant, i.e.,

$\displaystyle {\mathbf{p}^+}^\ast {\bm \Gamma}\mathbf{p}^+$	$\displaystyle =$	$\displaystyle {\mathbf{p}^-}^\ast {\bm \Gamma}\mathbf{p}^-$
$\displaystyle \implies \quad \mathbf{A}^\ast {\bm \Gamma}\mathbf{A}$	$\displaystyle =$	$\displaystyle {\bm \Gamma} \protect$	(C.124)

where ${\bm \Gamma}$ is any Hermitian, positive-definite matrix (which has an interpretation as a generalized junction admittance). The form $x^\ast {\bm \Gamma} x$ is by definition the square of the elliptic norm of

induced by ${\bm \Gamma}$ , or $\vert\vert\,x\,\vert\vert _{\bm \Gamma}^2 = x^\ast {\bm \Gamma}x$ . Setting ${\bm \Gamma}=\mathbf{I}$ , we obtain that $\mathbf{A}$ must be unitary. This is the case commonly used in current FDN practice.

The following theorem gives a general characterization of lossless scattering:

Theorem: A scattering matrix (FDN feedback matrix) $\mathbf{A}$ is lossless if and only if its eigenvalues lie on the unit circle and its eigenvectors are linearly independent.

Proof: Since ${\bm \Gamma}$ is positive definite, it can be factored (by the Cholesky factorization) into the form ${\bm \Gamma}= \mathbf{U}^\ast \mathbf{U}$ , where $\mathbf{U}$ is an upper triangular matrix, and $\mathbf{U}^\ast$ denotes the Hermitian transpose of $\mathbf{U}$ , i.e., $\mathbf{U}^\ast \isdef \overline{\mathbf{U}}^T$ . Since ${\bm \Gamma}$ is positive definite, $\mathbf{U}$ is nonsingular and can be used as a similarity transformation matrix. Applying the Cholesky decomposition ${\bm \Gamma}= \mathbf{U}^\ast \mathbf{U}$ in Eq. $\,$ (C.125) yields

$\begin{eqnarray*} & & \mathbf{A}^\ast {\bm \Gamma}\mathbf{A}= {\bm \Gamma}\\ &\... ...\implies& \tilde{\mathbf{A}}^\ast \tilde{\mathbf{A}}= \mathbf{I} \end{eqnarray*}$

where $\mathbf{U}^{-\ast}\isdef (\mathbf{U}^{-1})^\ast$ , and

$\displaystyle \tilde{\mathbf{A}}\isdef \mathbf{U}\mathbf{A}\mathbf{U}^{-1}$

is similar to $\mathbf{A}$ using $\mathbf{U}^{-1}$ as the similarity transform matrix. Since $\tilde{\mathbf{A}}$ is unitary, its eigenvalues have modulus 1. Hence, the eigenvalues of every lossless scattering matrix lie on the unit circle in the

plane. It readily follows from similarity to $\tilde{\mathbf{A}}$ that $\mathbf{A}$ admits

linearly independent eigenvectors. In fact, $\tilde{\mathbf{A}}$ is a normal matrix ( $\mathbf{A}\tilde{\mathbf{A}}= \tilde{\mathbf{A}}\mathbf{A}$ ), since every unitary matrix is normal, and normal matrices admit a basis of linearly independent eigenvectors [346].

Conversely, assume $\vert\lambda\vert = 1$ for each eigenvalue of $\mathbf{A}$ , and that there exists a matrix $\mathbf{E}$ of linearly independent eigenvectors of $\mathbf{A}$ . The matrix $\mathbf{E}$ diagonalizes $\mathbf{A}$ to give $\mathbf{E}^{-1}\mathbf{A}\mathbf{E}= \mathbf{D}$ , where $\mathbf{D}=$ diag $(\lambda_1,\dots,\lambda_N)$ . Taking the Hermitian transform of this equation gives $\mathbf{E}^\ast \mathbf{A}^\ast \mathbf{E}^{-\ast}= \mathbf{D}^\ast$ . Multiplying, we obtain $\mathbf{E}^\ast \mathbf{A}^\ast \mathbf{E}^{-\ast}\mathbf{E}^{-1}\mathbf{A}\ma... ...t \mathbf{E}^{-\ast}\mathbf{E}^{-1}\mathbf{A}=\mathbf{E}^{-\ast}\mathbf{E}^{-1}$ . Thus, (C.125) is satisfied for ${\bm \Gamma}=\mathbf{E}^{-\ast}\mathbf{E}^{-1}$ which is Hermitian and positive definite. $\Box$

Thus, lossless scattering matrices may be fully parametrized as $\mathbf{A} = \mathbf{E}^{-1}\mathbf{D}\mathbf{E}$ , where $\mathbf{D}$ is any unit-modulus diagonal matrix, and $\mathbf{E}$ is any invertible matrix. In the real case, we have $\mathbf{D}=$ diag $(\pm 1)$ and $\mathbf{E}\in\Re^{N\times N}$ .

Note that not all lossless scattering matrices have a simple physical interpretation as a scattering matrix for an intersection of lossless reflectively terminated waveguides. In addition to these cases (generated by all non-negative branch impedances), there are additional cases corresponding to sign flips and branch permutations at the junction. In terms of classical network theory [35], such additional cases can be seen as arising from the use of ``gyrators'' and/or ``circulators'' at the scattering junction [433]).

Waveguide Transformers and Gyrators

The ideal transformer, depicted in Fig. C.37 a, is a lossless two-port electric circuit element which scales up voltage by a constant [110,35]. In other words, the voltage at port 2 is always times the voltage at port 1. Since power is voltage times current, the current at port 2 must be times the current at port 1 in order for the transformer to be lossless. The scaling constant is called the turns ratio because transformers are built by coiling wire around two sides of a magnetically permeable torus, and the number of winds around the port 2 side divided by the winding count on the port 1 side gives the voltage stepping constant .

**Figure C.37:** a) Two-port description of the ideal transformer with ``turns ratio'' . b) Corresponding wave digital transformer.
$\includegraphics[width=\twidth]{eps/lTransformer}$

In the case of mechanical circuits, the two-port transformer relations appear as

$\begin{eqnarray*} F_2(s) &=& g F_1(s) \\ [5pt] V_2(s) &=& \frac{1}{g} V_1(s) \end{eqnarray*}$

where and denote force and velocity, respectively. We now convert these transformer describing equations to the wave variable formulation. Let and denote the wave impedances on the port 1 and port 2 sides, respectively, and define velocity as positive into the transformer. Then

$\begin{eqnarray*} f^{{+}}_1(t) &=& \frac{f_1(t) + R_1 v_1(t)}{2} \\ f^{{-}}_1(t... ...2(t)}{2} \\ &=& \frac{1}{g} \frac{f_2(t) + R_1 g^2 v_2(t)}{2}. \end{eqnarray*}$

Similarly,

$\begin{eqnarray*} f^{{+}}_2(t) &=& \frac{f_2(t) + R_2 v_2(t)}{2} \\ f^{{-}}_2(t... ...1(t)}{2} \\ &=& g \frac{f_1(t) + R_2 \frac{1}{g^2} v_1(t)}{2}. \end{eqnarray*}$

We see that choosing

$\displaystyle g^2 = \frac{R_2}{R_1}$

eliminates the scattering terms and gives the simple relations

$\begin{eqnarray*} f^{{-}}_2(t) &=& g f^{{+}}_1(t)\\ [5pt] f^{{-}}_1(t) &=& \frac{1}{g}f^{{+}}_2(t). \end{eqnarray*}$

The corresponding wave flow diagram is shown in Fig. C.37 b.

Thus, a transformer with a voltage gain corresponds to simply changing the wave impedance from to , where $g=\sqrt{R_2/R_1}$ . Note that the transformer implements a change in wave impedance without scattering as occurs in physical impedance steps (§C.8).

Gyrators

Another way to define the ideal waveguide transformer is to ask for a two-port element that joins two waveguide sections of differing wave impedance in such a way that signal power is preserved and no scattering occurs. From Ohm's Law for traveling waves (Eq. $\,$ (6.6)), and from the definition of power waves (§C.7.5), we see that to bridge an impedance discontinuity between $R_{i-1}$ and with no power change and no scattering requires the relations

$\displaystyle \frac{[f^{{+}}_i]^2}{R_i} = \frac{[f^{{+}}_{i-1}]^2}{R_{i-1}}, \qquad\qquad \frac{[f^{{-}}_i]^2}{R_i} = \frac{[f^{{-}}_{i-1}]^2}{R_{i-1}}.$

Therefore, the junction equations for a transformer can be chosen as

$\displaystyle f^{{+}}_i= g_i\, f^{{+}}_{i-1}\qquad\qquad f^{{-}}_{i-1}= g_i^{-1}\, f^{{-}}_i \protect$

(C.125)

where

$\displaystyle g_i \isdef \pm\sqrt{\frac{R_i}{R_{i-1}}} \protect$

(C.126)

Choosing the negative square root for $g_i^{-1}$ gives a gyrator [35]. Gyrators are often used in electronic circuits to replace inductors with capacitors. The gyrator can be interpreted as a transformer in cascade with a dualizer [433]. A dualizer converts one from wave variable type (such as force) to the other (such as velocity) in the waveguide.

The dualizer is readily derived from Ohm's Law for traveling waves:

$\begin{eqnarray*} f^{{+}}\eqsp Rv^{+}, \qquad f^{{-}}\eqsp -Rv^{-}\\ [5pt] \Lon... ...i\eqsp Rv^{+}_{i-1}, \qquad v^{-}_{i-1} \eqsp -R^{-1} f^{{-}}_i \end{eqnarray*}$

In this case, velocity waves in section are converted to force waves in section , and vice versa (all at wave impedance ). The wave impedance can be changed as well by cascading a transformer with the dualizer, which changes to $R\sqrt{R_i/R}=\sqrt{RR_i}$ (where we assume $R=R_{i-1}$ ). Finally, the velocity waves in section can be scaled to equal their corresponding force waves by introducing a transformer $g=\sqrt{1/R}$ on the left, which then coincides Eq. $\,$ (C.126) (but with a minus sign in the second equation).

The Digital Waveguide Oscillator

In this section, adapted from [460], a digital sinusoidal oscillator derived from digital waveguide theory is described which has good properties for VLSI implementation. Its main features are no wavetable and a computational complexity of only one multiply per sample when amplitude and frequency are constant. Three additions are required per sample. A piecewise exponential amplitude envelope is available for the cost of a second multiplication per sample, which need not be as expensive as the tuning multiply. In the presence of frequency modulation (FM), the amplitude coefficient can be varied to exactly cancel amplitude modulation (AM) caused by changing the frequency of oscillation.

Additive Synthesis

One of the very first computer music techniques introduced was additive synthesis [379]. It is based on Fourier's theorem which states that any sound can be constructed from elementary sinusoids, such as are approximately produced by carefully struck tuning forks. Additive synthesis attempts to apply this theorem to the synthesis of sound by employing large banks of sinusoidal oscillators, each having independent amplitude and frequency controls. Many analysis methods, e.g., the phase vocoder, have been developed to support additive synthesis. A summary is given in [424].

While additive synthesis is very powerful and general, it has been held back from widespread usage due to its computational expense. For example, on a single DSP56001 digital signal-processing chip, clocked at 33 MHz, only about sinusoidal partials can be synthesized in real time using non-interpolated, table-lookup oscillators. Interpolated table-lookup oscillators are much more expensive, and when all the bells and whistles are added, and system overhead is accounted for, only around fully general, high-quality partials are sustainable at KHz on a DSP56001 (based on analysis of implementations provided by the NeXT Music Kit).

At CD-quality sampling rates, the note A1 on the piano requires $22050/55\approx 400$ sinusoidal partials, and at least the low-frequency partials should use interpolated lookups. Assuming a worst-case average of partials per voice, providing 32-voice polyphony requires partials, or around DSP chips, assuming we can pack an average of partials into each DSP. A more reasonable complement of DSP chips would provide only -voice polyphony which is simply not enough for a piano synthesis. However, since DSP chips are getting faster and cheaper, DSP-based additive synthesis looks viable in the future.

The cost of additive synthesis can be greatly reduced by making special purpose VLSI optimized for sinusoidal synthesis. In a VLSI environment, major bottlenecks are wavetables and multiplications. Even if a single sinusoidal wavetable is shared, it must be accessed sequentially, inhibiting parallelism. The wavetable can be eliminated entirely if recursive algorithms are used to synthesize sinusoids directly.

Digital Sinusoid Generators

In [168], three techniques were examined for generating sinusoids digitally by means of recursive algorithms.^C.12 The recursions can be interpreted as implementations of second-order digital resonators in which the damping is set to zero. The three methods considered were (1) the 2D rotation (2DR), or complex multiply (also called the ``coupled form''), (2) the modified coupled form (MCF), or ``magic circle'' algorithm,^C.13which is similar to (1) but with better numerical behavior, and (3) the direct-form, second-order, digital resonator (DFR) with its poles set to the unit circle.

These three recursions may be defined as follows:

$\begin{displaymath} \begin{array}{crcll} (1) & x_1(n) &=& c_nx_1(n-1) + s_nx_2(n... ..._2(n-1) & \\ & x_2(n) &=& x_1(n-1) & \mbox{(DFR)} \end{array}\end{displaymath}$

where $c_n\isdef \cos(2\pi f_n T)$ , $s_n\isdef \sin(2\pi f_n T)$ ,

is the instantaneous frequency of oscillation (Hz) at time sample

, and

is the sampling period in seconds. The magic circle parameter is $\epsilon=2\sin(\pi f_n T)$ .

The digital waveguide oscillator appears to have the best overall properties yet seen for VLSI implementation. This structure, introduced in [460], may be derived from the theory of digital waveguides (see Appendix C, particularly §C.9, and [433,464]). Any second-order digital filter structure can be used as a starting point for developing a corresponding sinusoidal signal generator, so in this case we begin with the second-order waveguide filter.

The Second-Order Waveguide Filter

The first step is to make a second-order digital filter with zero damping by abutting two unit-sample sections of waveguide medium, and terminating on the left and right with perfect reflections, as shown in Fig.C.38. The wave impedance in section is given by $R_i=\rho c/A_i$ , where $\rho$ is air density, is the cross-sectional area of tube section , and is sound speed. The reflection coefficient is determined by the impedance discontinuity via . It turns out that to obtain sinusoidal oscillation, one of the terminations must provide an inverting reflection while the other is non-inverting.

**Figure C.38:** The second-order, lossless, digital waveguide oscillator, built using two *acoustic tube* sections.
$\includegraphics[width=\twidth]{eps/wgf}$

At the junction between sections and , the signal is partially transmitted and partially reflected such that energy is conserved, i.e., we have lossless scattering. The formula for the reflection coefficient can be derived from the physical constraints that (1) pressure is continuous across the junction, and (2) there is no net flow into or out of the junction. For traveling pressure waves $p^\pm(t)$ and volume-velocity waves $u^\pm(t)$ , we have and . The physical pressure and volume velocity are obtained by summing the traveling-wave components.

The discrete-time simulation for the physical system of Fig.C.38 is shown in Fig.C.39. The propagation time from the junction to a reflecting termination and back is one sample period. The half sample delay from the junction to the reflecting termination has been commuted with the termination and combined with the half sample delay to the termination. This is a special case of a ``half-rate'' waveguide filter [433].

**Figure C.39:** The second-order, lossless, waveguide filter.
$\includegraphics[width=\twidth]{eps/dwgf}$

Since only two samples of delay are present, the digital system is at most second order, and since the coefficients are real, at most one frequency of oscillation is possible in $(0,\pi)$ .

The scattering junction shown in the figure is called the Kelly-Lochbaum junction in the literature on lattice and ladder digital filters [173]. While it is the most natural from a physical point of view, it requires four multiplies and two additions for its implementation.

It is well known that lossless scattering junctions can be implemented in a variety of equivalent forms, such as the two-multiply and even one-multiply junctions. However, most have the disadvantage of not being normalized in the sense that changing the reflection coefficient changes the amplitude of oscillation. This can be understood physically by noting that a change in implies a change in . Since the signal power contained in a waveguide variable, say , is $\left[p_1^+(n)\right]^2/R_1$ , we find that modulating the reflection coefficient corresponds to modulating the signal energy represented by the signal sample in at least one of the two delay elements. Since energy is proportional to amplitude squared, energy modulation implies amplitude modulation.

The well-known normalization procedure is to replace the traveling pressure waves $p^\pm$ by ``root-power'' pressure waves $\tilde{p}^\pm = p^\pm/\sqrt{R}$ so that signal power is just the square of a signal sample $(\tilde{p}^\pm)^2$ . When this is done, the scattering junction transforms from the Kelly-Lochbaum or one-multiply form into the normalized ladder junction in which the reflection coefficients are again $\pm k$ , but the forward and reverse transmission coefficients become $\sqrt{1-k^2}$ . Defining $k=\sin(\theta)$ , the transmission coefficients can be seen as $\cos(\theta)$ , and we arrive essentially at the coupled form, or two-dimensional vector rotation considered in [168].

An alternative normalization technique is based on the digital waveguide transformer (§C.16). The purpose of a ``transformer'' is to ``step'' the force variable (pressure in our example) by some factor without scattering and without affecting signal energy. Since traveling signal power is proportional to pressure times velocity , it follows that velocity must be stepped by the inverse factor to keep power constant. This is the familiar behavior of transformers for analog electrical circuits: voltage is stepped up by the ``turns ratio'' and current is stepped down by the reciprocal factor. Now, since , traveling signal power is equal to . Therefore, stepping up pressure through a transformer by the factor corresponds to stepping up the wave impedance by the factor . In other words, the transformer raises pressure and decreases volume velocity by raising the wave impedance (narrowing the acoustic tube) like a converging cone.

If a transformer is inserted in a waveguide immediately to the left, say, of a scattering junction, it can be used to modulate the wave impedance ``seen'' to the left by the junction without having to use root-power waves in the simulation. As a result, the one-multiply junction can be used for the scattering junction, since the junction itself is not normalized. Since the transformer requires two multiplies, a total of three multiplies can effectively implement a normalized junction, where four were needed before. Finally, in just this special case, one of the transformer coefficients can be commuted with the delay element on the left and combined with the other transformer coefficient. For convenience, the coefficient on the left is commuted into the junction so it merely toggles the signs of inputs to existing summers. These transformations lead to the final form shown in Fig.C.40.

**Figure C.40:** The transformer-normalized, digital waveguide oscillator.
$\includegraphics[scale=0.9]{eps/tnwgoAI}$

The ``tuning coefficient'' is given by $C(n)=\cos(2\pi f_nT)$ , where is the desired oscillation frequency in Hz at sample (in the undamped case), and is the sampling period in seconds. The ``amplitude coefficient'' is $G(n) = r_n g_n/g_{n-1}$ , where growth or decay factor per sample ( $r_n\equiv 1$ for constant amplitude),^C.14 and is the normalizing transformer ``turns ratio'' given by

$\displaystyle g_n=\sqrt{\frac{1-C(n)}{1+C(n)}}.$

When both amplitude and frequency are constant, we have $G(n)\equiv 1$ , and only the tuning multiply is operational. When frequency changes, the amplitude coefficient deviates from unity for only one time sample to normalize the oscillation amplitude.

When amplitude and frequency are constant, there is no gradual exponential growth or decay due to round-off error. This happens because the only rounding is at the output of the tuning multiply, and all other computations are exact. Therefore, quantization in the tuning coefficient can only cause quantization in the frequency of oscillation. Note that any one-multiply digital oscillator should have this property. In contrast, the only other known normalized oscillator, the coupled form, does exhibit exponential amplitude drift because it has two coefficients $c=\cos(\theta)$ and $s=\sin(\theta)$ which, after quantization, no longer obey for most tunings.

Application to FM Synthesis

The properties of the new oscillator appear well suited for FM applications in VLSI because of the minimized computational expense. However, in this application there are issues to be resolved regarding conversion from modulator output to carrier coefficients. Preliminary experiments indicate that FM indices less than are well behaved when the output of a modulating oscillator simply adds to the coefficient of the carrier oscillator (bypassing the exact FM formulas). Approximate amplitude normalizing coefficients have also been derived which provide a first-order approximation to the exact AM compensation at low cost. For music synthesis applications, we believe a distortion in the details of the FM instantaneous frequency trajectory and a moderate amount of incidental AM can be tolerated since they produce only second-order timbral effects in many situations.

Digital Waveguide Resonator

Converting a second-order oscillator into a second-order filter requires merely introducing damping and defining the input and output signals. In Fig.C.40, damping is provided by the coefficient , which we will take to be a constant

$\displaystyle G(n)\equiv\message{CHANGE eqv TO equiv IN SOURCE}g.$

When

, the oscillator decays exponentially to zero from any initial conditions. The two delay elements constitute the state of the resonator. Let us denote by

the output of the delay element on the left in Fig.C.40 and let

be the delay-element output on the right. In general, an output signal

may be formed as any linear combination of the state variables:

$\displaystyle y(n) = c_1 x_1(n) + c_2 x_2(n)$

Similarly, input signals

may be summed into the state variables

scaled by arbitrary gain factors

The foregoing modifications to the digital waveguide oscillator result in the so-called digital waveguide resonator (DWR) [304]:

$\displaystyle \tilde{x}_1(n)$	$\displaystyle =$	$\displaystyle g x_1(n)\protect$	(C.127)
$\displaystyle v(n)$	$\displaystyle =$	$\displaystyle c[\tilde{x}_1(n) + x_2(n)]$	(C.128)
$\displaystyle x_1(n+1)$	$\displaystyle =$	$\displaystyle v(n) - x_2(n) + \tilde{b}_1 u(n)$	(C.129)
$\displaystyle x_2(n+1)$	$\displaystyle =$	$\displaystyle \tilde{x}_1(n) + v(n) + b_2 u(n)$	(C.130)
$\displaystyle y(n)$	$\displaystyle =$	$\displaystyle c_1 x_1(n) + c_2 x_2(n) \protect$	(C.131)

where, as derived in the next section, the coefficients are given by

$\displaystyle g$	$\displaystyle =$	$\displaystyle r^2\protect$	(C.132)
$\displaystyle \tilde{b}_1$	$\displaystyle =$	$\displaystyle b_1 \sqrt{\frac{1-c}{1+c}}$	(C.133)
$\displaystyle c$	$\displaystyle =$	$\displaystyle \sqrt{\frac{1}{1 + \frac{\tan^2(\theta)(1+g)^2+(1-g)^2}{4g}}}$	(C.134)
	$\displaystyle \approx$	$\displaystyle 1 - \frac{\tan^2(\theta)(1+g)^2 + (1-g)^2}{8g}. \protect$	(C.135)

where $\lambda=r e^{j\theta}$ denotes one desired pole (the other being at $\overline{\lambda}$ ). Note that $c=\cos(\theta)$ when

(undamped case). The DWR requires only two multiplies per sample. As seen earlier, when the decay time is set to $\infty$ (

), one of the multiplies disappears, leaving only one multiply per sample for sinusoidal oscillation.

Figure C.41 shows an overlay of initial impulse responses for the three resonators discussed above. The decay factor was set to , and the output of each multiplication was quantized to 16 bits, as were all coefficients. The three waveforms sound and look identical. (There are small differences, however, which can be seen by plotting the differences of pairs of waveforms.)

**Figure:** Overlay of three resonator impulse responses, with $\mathbf {B}^T=[1,0]$ and $\mathbf {C}=[0,1]$ , for the (1) complex-multiply resonator (labeled ``2DR'' for ``2D rotation''), (2) modified coupled form (MCF), and (3) second-order digital waveguide resonator (DWR).
$\includegraphics[width=\twidth]{eps/tosc16}$

Figure C.42 shows the same impulse-response overlay but with and only 4 significant bits in the coefficients and signals. The complex multiply oscillator can be seen to decay toward zero due to coefficient quantization ( ). The MCF and DWR remain steady at their initial amplitude. All three suffer some amount of tuning perturbation.

**Figure:** Overlay of three resonator impulse responses, as in Fig.C.41, but with and quantization of coefficients and signals to 4 significant bits.
$\includegraphics[width=\twidth]{eps/tosc4}$

State-Space Analysis

We will now use state-space analysis^C.15[449] to determine Equations (C.133-C.136).

From Equations (C.128-C.132),

$\displaystyle x_1(n+1) = c[g x_1(n) + x_2(n)] - x_2(n) = c\,g x_1(n) + (c-1) x_2(n)$

and

$\displaystyle x_2(n+1) = g x_1(n) + c[g x_1(n) + x_2(n)] = (1+c) g x_1(n) + c\,x_2(n)$

In matrix form, the state time-update can be written

$\displaystyle \left[\begin{array}{c} x_1(n+1) \\ [2pt] x_2(n+1) \end{array}\rig... ...bf{A} \left[\begin{array}{c} x_1(n) \\ [2pt] x_2(n) \end{array}\right] \protect$

(C.136)

or, in vector notation,

$\displaystyle \underline{x}(n+1)$	$\displaystyle =$	$\displaystyle \mathbf{A}\underline{x}(n) + \mathbf{B}u(n)$	(C.137)
$\displaystyle y(n)$	$\displaystyle =$	$\displaystyle \mathbf{C}\underline{x}(n)$	(C.138)

where we have introduced an input signal

, which sums into the state vector via the $2\times 1$ (or $2\times N_u$ ) vector $\mathbf{B}$ . The output signal is defined as the $1\times 2$ vector $\mathbf{C}$ times the state vector $\underline{x}(n)$ . Multiple outputs may be defined by choosing $\mathbf{C}$ to be an $N_y \times 2$ matrix.

A basic fact from linear algebra is that the determinant of a matrix is equal to the product of its eigenvalues. As a quick check, we find that the determinant of is

$\displaystyle \det{\mathbf{A}} = c^2g - g(c+1)(c-1) = c^2g - g(c^2-1) = c^2g - gc^2+g = g. \protect$

(C.139)

When the eigenvalues ${\lambda_i}$ of $\mathbf{A}$ (system poles) are complex, then they must form a complex conjugate pair (since $\mathbf{A}$ is real), and we have $\det{\mathbf{A}} = \vert{\lambda_i}\vert^2 = g$ . Therefore, the system is stable if and only if $\vert g\vert<1$ . When making a digital sinusoidal oscillator from the system impulse response, we have $\vert g\vert=1$ , and the system can be said to be ``marginally stable''. Since an undriven sinusoidal oscillator must not lose energy, and since every lossless state-space system has unit-modulus eigenvalues (consider the modal representation), we expect $\left\vert\det{A}\right\vert=1$ , which occurs for

Note that $\underline{x}(n) = \mathbf{A}^n\underline{x}(0)$ . If we diagonalize this system to obtain $\tilde{\mathbf{A}}= \mathbf{E}^{-1}\mathbf{A}\mathbf{E}$ , where $\tilde{\mathbf{A}}=$ diag $[\lambda_1,\lambda_2]$ , and $\mathbf{E}$ is the matrix of eigenvectors of $\mathbf{A}$ , then we have

$\displaystyle \tilde{\underline{x}}(n) = \tilde{A}^n\,\tilde{\underline{x}}(0) ... ...eft[\begin{array}{c} \tilde{x}_1(0) \\ [2pt] \tilde{x}_2(0) \end{array}\right]$

where $\tilde{\underline{x}}(n) \isdef \mathbf{E}^{-1}\underline{x}(n)$ denotes the state vector in these new ``modal coordinates''. Since $\tilde{A}$ is diagonal, the modes are decoupled, and we can write

$\begin{eqnarray*} \tilde{x}_1(n) &=& \lambda_1^n\,\tilde{x}_1(0)\\ \tilde{x}_2(n) &=& \lambda_2^n\,\tilde{x}_2(0) \end{eqnarray*}$

If this system is to generate a real sampled sinusoid at radian frequency $\omega$ , the eigenvalues $\lambda_1$ and $\lambda_2$ must be of the form

$\begin{eqnarray*} \lambda_1 &=& e^{j\omega T}\\ \lambda_2 &=& e^{-j\omega T}, \end{eqnarray*}$

(in either order) where $\omega$ is real, and denotes the sampling interval in seconds.

Thus, we can determine the frequency of oscillation $\omega$ (and verify that the system actually oscillates) by determining the eigenvalues $\lambda_i$ of . Note that, as a prerequisite, it will also be necessary to find two linearly independent eigenvectors of (columns of $\mathbf{E}$ ).

Eigenstructure

Starting with the defining equation for an eigenvector $\underline{e}$ and its corresponding eigenvalue $\lambda$ ,

$\displaystyle \mathbf{A}\underline{e}_i = {\lambda_i}\underline{e}_i,\quad i=1,2$

we get, using Eq. $\,$ (C.137),

$\displaystyle \left[\begin{array}{cc} gc & c-1 \\ [2pt] gc+g & c \end{array}\ri... ...n{array}{c} {\lambda_i} \\ [2pt] {\lambda_i}\eta_i \end{array}\right]. \protect$

(C.140)

We normalized the first element of $\underline{e}_i$ to 1 since $g\underline{e}_i$ is an eigenvector whenever $\underline{e}_i$ is. (If there is a missing solution because its first element happens to be zero, we can repeat the analysis normalizing the second element to 1 instead.)

Equation (C.141) gives us two equations in two unknowns:

$\displaystyle gc+\eta_i(c-1)$	$\displaystyle =$	$\displaystyle {\lambda_i} \protect$	(C.141)
$\displaystyle g(1+c) +c\eta_i$	$\displaystyle =$	$\displaystyle {\lambda_i}\eta_i$	(C.142)

Substituting the first into the second to eliminate ${\lambda_i}$ , we get

$\begin{eqnarray*} g+gc+c\eta_i &=& [gc+\eta_i(c-1)]\eta_i = gc\eta_i + \eta_i^2 ... ...{g\left(\frac{1+c}{1-c}\right) - \frac{c^2(1-g)^2}{4(1-c)^2}}. \end{eqnarray*}$

As approaches (no damping), we obtain

$\displaystyle \eta_i = \pm j\sqrt{\frac{1+c}{1-c}} \qquad \hbox{(when $g=1$)}.$

Thus, we have found both eigenvectors:

$\begin{eqnarray*} \underline{e}_1&=&\left[\begin{array}{c} 1 \\ [2pt] \eta \end{... ...t{g\left(\frac{1+c}{1-c}\right) - \frac{c^2(1-g)^2}{4(1-c)^2}} \end{eqnarray*}$

They are linearly independent provided $\eta\neq0$ . In the undamped case (), this holds whenever $c\neq -1$ . The eigenvectors are finite when $c\neq 1$ . Thus, the nominal range for is the interval $c\in(-1,1)$ .

We can now use Eq. $\,$ (C.142) to find the eigenvalues:

$\begin{eqnarray*} {\lambda_i}&=& gc+ \eta_i(c-1)\\ &=& gc+ \frac{(1-g)c}{2}\pm ... ... \pm j\sqrt{g(1-c^2) - \left[\frac{c(1-g)}{2}\right]^2} \protect \end{eqnarray*}$

Damping and Tuning Parameters

The tuning and damping of the resonator impulse response are governed by the relation

$\displaystyle {\lambda_i}= e^{\frac{T}{\tau}} e^{\pm j\omega T}$

where

denotes the sampling interval, $\tau$ is the time constant of decay, and $\omega$ is the frequency of oscillation in radians per second. The eigenvalues are presumed to be complex, which requires, from Eq. $\,$ (C.144),

$\displaystyle g(1-c^2) \geq\frac{c^2(1-g)^2}{4} \,\,\Rightarrow\,\,c^2 \leq \frac{4g}{(1+g)^2}$

To obtain a specific decay time-constant $\tau$ , we must have

$\begin{eqnarray*} e^{-2T/\tau} &=& \left\vert{\lambda_i}\right\vert^2 = c^2\left... ...left[g(1-c^2) - c^2\left(\frac{1-g}{2}\right)^2\right]\\ &=& g \end{eqnarray*}$

Therefore, given a desired decay time-constant $\tau$ (and the sampling interval ), we may compute the damping parameter for the digital waveguide resonator as

$\displaystyle \zbox {g = e^{-2T/\tau}.}$

Note that this conclusion follows directly from the determinant analysis of Eq. $\,$ (C.140), provided it is known that the poles form a complex-conjugate pair.

To obtain a desired frequency of oscillation, we must solve

$\begin{eqnarray*} \theta = \omega T &=& \tan^{-1}\left[\frac{\sqrt{g(1-c^2) - [... ...,\tan^2{\theta} &=& \frac{g(1-c^2) - [c(1-g)/2]^2}{[c(1+g)/2]^2} \end{eqnarray*}$

for , which yields

$\displaystyle \zbox { c= \sqrt{\frac{1}{1 + \frac{\tan^2(\theta)(1+g)^2+(1-g)^2}{4g}}} \approx 1 - \frac{\tan^2(\theta)(1+g)^2 + (1-g)^2}{8g}. }$

Note that this reduces to $c=\cos(\theta)$ when

(undamped case).

Eigenvalues in the Undamped Case

When , the eigenvalues reduce to

$\displaystyle {\lambda_i}= c\pm j\sqrt{1-c^2}$

Assuming $\left\vert c\right\vert<1$ , the eigenvalues can be expressed as

$\displaystyle {\lambda_i}= c\pm j\sqrt{1-c^2} = \cos(\theta) \pm j\sin(\theta) = e^{\pm j\theta} \protect$

(C.143)

where $\theta=\omega T$ denotes the angular advance per sample of the oscillator. Since $c\in(-1,1)$ corresponds to the range $\theta\in(-\pi,\pi)$ , we see that

in this range can produce oscillation at any digital frequency.

For $\left\vert c\right\vert>1$ , the eigenvalues are real, corresponding to exponential growth and/or decay. (The values $c=\pm 1$ were excluded above in deriving Eq. $\,$ (C.144).)

In summary, the coefficient in the digital waveguide oscillator () and the frequency of sinusoidal oscillation $\omega$ is simply

$\displaystyle \fbox{$\displaystyle c= \cos(\omega T)$}.$

Summary

This section introduced and analyzed the digital waveguide oscillator and resonator, as well as some related algorithms. As a recursive algorithm for digital sinusoid generation, it has excellent properties for VLSI implementation. It is like the 2D rotation (complex multiply) in that it offers instantaneous amplitude from its state and constant amplitude in the presence of frequency modulation. However, its implementation requires only two multiplies per sample instead of four. When used as a constant-frequency oscillator, it requires only one multiply per sample.

Matlab Sinusoidal Oscillator Implementations

This section provides Matlab/Octave program listings for the sinusoidal resonator/oscillator algorithms discussed above:

Planar 2D Rotation (2DR) (``complex multiply'')
Modified Coupled Form (MCF) (``Magic Circle'')
Digital Waveguide Resonator (DWR)

The test program computes an impulse response of each resonator, and plots them overlaid for comparison. Next, it optionally plays each one as a sound and plots them individually.

% Filter test program in matlab

%N = 300;   % Number of samples to generate
N = 3000;   % Number of samples to generate
f = 100;   % Desired oscillation frequency (Hz)
fs = 8192; % Audio sampling rate (Hz)
%R = .99;   % Decay factor (1=never, 0=instant)
R = 1;   % Decay factor (1=never, 0=instant)
b1 = 1;    % Input gain to state variable x(n)
b2 = 0;    % Input gain to state variable y(n)

%nd = 16;   % Number of significant digits to use
nd = 4;   % Number of significant digits to use
base = 2;  % Mantissa base (two significant figures)
           % (see 'help chop' in Matlab)

u = [1,zeros(1,N-1)]; % driving input test signal
theta = 2*pi*f/fs; % resonance frequency, rad/sample

% ================================================
% 2D PLANAR ROTATION (COMPLEX MULTIPLY)

x1 = zeros(1,N); % Predeclare saved-state arrays
y1 = zeros(1,N);

x1(1) = 0;      % Initial condition
y1(1) = 0;      % Initial condition

c = chop(R*cos(theta),nd,base); % coefficient 1
s = chop(R*sin(theta),nd,base); % coefficient 2

for n=1:N-1,
  x1(n+1) = chop(   c*x1(n) - s*y1(n) + b1*u(n), nd,base);
  y1(n+1) = chop(   s*x1(n) + c*y1(n) + b2*u(n), nd,base);
end

% ================================================
% MODIFIED COUPLED FORM ("MAGIC CIRCLE")
% (ref: http://ccrma.stanford.edu/~jos/wgo/ )

x2 = zeros(1,N); % Predeclare saved-state arrays
y2 = zeros(1,N);

x2(1) = 0.0;     % Initial condition
y2(1) = 0.0;     % Initial condition

e = chop(2*sin(theta/2),nd,base); % tuning coefficient

for n=1:N-1,
  x2(n+1) = chop(R*(x2(n)-e*y2(n))+b1*u(n),nd,base);
  y2(n+1) = chop(R*(e*x2(n+1)+y2(n))+b2*u(n),nd,base);
end

% ================================================
% DIGITAL WAVEGUIDE RESONATOR (DWR)

x3 = zeros(1,N); % Predeclare saved-state arrays
y3 = zeros(1,N);

x3(1) = 0;      % Initial condition
y3(1) = 0;      % Initial condition

g = R*R;        % decay coefficient
t = tan(theta); % toward tuning coefficient
cp = sqrt(g/(g + t^2*(1+g)^2/4 + (1-g)^2/4)); % exact
%cp = 1 - (t^2*(1+g)^2 + (1-g)^2)/(8*g); % good approx
b1n = b1*sqrt((1-cp)/(1+cp)); % input scaling

% Quantize coefficients:
cp = chop(cp,nd,base);
g = chop(g,nd,base);
b1n = chop(b1n,nd,base);

for n=1:N-1,
  gx3 = chop(g*x3(n), nd,base); % mpy 1
  y3n = y3(n);
  temp = chop(cp*(gx3 + y3n), nd,base); % mpy 2
  x3(n+1) = temp - y3n + b1n*u(n);
  y3(n+1) = gx3 + temp + b2*u(n);
end

% ================================================
% playandplot.m

figure(4);
title('Impulse Response Overlay');
ylabel('Amplitude');
xlabel('Time (samples)');
alt=1;%alt = (-1).^(0:N-1); % for visibility
plot([y1',y2',(alt.*y3)']); grid('on');
legend('2DR','MCF','WGR');
title('Impulse Response Overlay');
ylabel('Amplitude');
xlabel('Time (samples)');

saveplot(sprintf('eps/tosc%d.ps',nd));

if 1
  playandplot(y1,u,fs,'2D rotation',1);
  playandplot(y2,u,fs,'Magic circle',2);
  playandplot(y3,u,fs,'WGR',3);
end

function playandplot(y,u,fs,ttl,fnum);
% sound(y,fs,16);
  figure(fnum);

  unwind_protect # protect graph state
  subplot(211);
  axis("labely");
  axis("autoy");
  xlabel("");
  title(ttl);
  ylabel('Amplitude');
  xlabel('Time (ms)');
  t = 1000*(0:length(u)-1)/fs;
  timeplot(t,u,'-');
  legend('input');
  subplot(212);
  ylabel('Amplitude');
  xlabel('Time (ms)');
  timeplot(t,y,'-');
  legend('output');
  unwind_protect_cleanup # restore graph state
    grid("off");
    axis("auto","label");
    oneplot();
  end_unwind_protect
endfunction

Note that the chop utility only exists in Matlab. In Octave, the following ``compatibility stub'' can be used:

function [y] = chop(x,nd,base)
y = x;

Faust Implementation

The function oscw (file osc.lib) implements a digital waveguide sinusoidal oscillator in the Faust language [453,154,170]. There is also oscr implementing the 2D rotation case, and oscs implementing the modified coupled form (magic circle).

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)

where

$\begin{displaymath} \begin{array}{rclrcl} c & \isdef & \mbox{sound speed} & \qqu... ...'_x & \isdef & \dfrac{\partial}{\partial x}p_x(t,x) \end{array}\end{displaymath}$

and

is the pressure at time

and radial position

along the cone axis (or wall). In terms of

(rather than

), 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

(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 is replaced by . Thus, the solution is a superposition of left- and right-going traveling wave components, scaled by :

$\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

may be interpreted as a fixed waveshape traveling to the right, (i.e., in the positive

direction), with speed

. Similarly, a function of

may be seen as a wave traveling to the left (negative

direction) at

meters per second. The point

corresponds to the tip of the cone (center of the sphere), and

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 seconds, which implies a spatial sampling interval of $X\isdeftext cT$ meters. Sampling is carried out mathematically by the change of variables

$\begin{displaymath} \begin{array}{rclrl} x& \to & x_m&=& x_0+ mX\\ t & \to & t_n&=& nT \end{array}\end{displaymath}$

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

$\begin{eqnarray*} 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]. \end{eqnarray*}$

Define

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

where

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

and

). 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 and . The symbol `` $z^{-1}$ '' denotes a one-sample delay.
$\includegraphics[width=\twidth]{eps/fcideal}$

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.
$\includegraphics[width=\twidth]{eps/fcone}$

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 denote the area of the surface of constant phase at radial coordinate in the tube. Then the total force acting on the surface due to pressure is , as shown in Fig.C.45.

**Figure C.45:** Differential volume element for the conical acoustic tube.
$\includegraphics[width=3in]{eps/fconic}$

The net force to the right across the volume element between and 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

and

, respectively. To apply Newton's second law equating net force to mass times acceleration, we need the mass of the volume element

$\begin{eqnarray*} 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, \end{eqnarray*}$

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 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

, this can be written

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

(C.148)

Note that

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, 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 and , 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.
$\includegraphics[width=\twidth]{eps/fconeparams}$

Referring to Fig.C.46, the area function can be written for any cone in