DSPRelated.com
Free Books
   Free Books    Physical Audio Signal Processing  

Wave Digital Filters

A Wave Digital Filter (WDF) [136] is a particular kind of digital filter based on physical modeling principles. Unlike most digital filter types, every delay element in a WDF can be interpreted physically as holding the current state of a mass or spring (or capacitor or inductor). WDFs can also be viewed as a particular kind of finite difference scheme having unusually good numerical properties [55]. (See Appendix D for an introduction to finite difference schemes.) WDFs have been applied often in music signal processing [394,338,556,361,348,554,555,56,527,523,484].

Wave digital filters were developed initially by Alfred Fettweis [135] in the late 1960s for digitizing lumped electrical circuits composed of inductors, capacitors, resistors, transformers, gyrators, circulators, and other elements of classical network theory [136]. The WDF approach is based on the traveling-wave formulation of lumped electrical elements introduced by Belevitch [34].

A WDF is constructed by interconnecting simple discrete-time models of individual masses, springs, and dashpots (or inductors, capacitors, and resistors). The rules for interconnecting the elementary models are based on scattering theory (discussed in §C.8). This is a direct result of the fact that all signals explicitly computed may be physically interpreted as traveling wave components of physical variables.

Wave Digital Elements

When modeling mechanical systems composed of masses, springs, and dashpots, it is best to begin with an electrical equivalent circuit. Equivalent circuits make clear the network-theoretic structure of the system, clearly indicating, for example, whether interacting elements should be connected in series or parallel. Each element of the equivalent circuit can then be replaced by a first-order wave digital element, and the elements are finally parallel or series connected by means of scattering-junction interfaces known as adaptors.

Wave digital elements may be derived from their describing differential equations (in continuous time) as follows:

  • First express all physical quantities (such as force and velocity) in terms of traveling-wave components. The traveling wave components are called wave variables. For example, the force $ f(n)$ on a mass is decomposed as $ f(n) = f^{{+}}(n)+f^{{-}}(n)$, where $ f^{{+}}(n)$ is regarded as a traveling wave propagating toward the mass, while $ f^{{-}}(n)$ is seen as the traveling component propagating away from the mass. A ``traveling wave'' view of force mediation is actually much closer to physical reality than any instantaneous model.

  • Second, digitize the resulting traveling-wave system using the bilinear transform. The bilinear transform is equivalent in the time domain to the trapezoidal rule for numerical integration (see §7.3.2).

  • Connect $ N$ elementary units together by means of $ N$-port scattering junctions. There are two basic types of scattering junction, one for parallel, and one for series connection. (See §C.8 for the theory of scattering junctions.)

The next section will examine the above steps in greater detail.

An important benefit of introducing wave variables prior to bilinear transformation is the elimination of delay-free loops when connecting elementary building blocks. In other words, any number of elementary models can be interconnected, in series or in parallel, and the resulting finite-difference scheme remains explicit (free of delay-free loops).

A Physical Derivation of Wave Digital Elements

This section provides a ``physical'' derivation of Wave Digital Filters (WDF), which contrasts somewhat with the more formal derivation common in the literature. The derivation is presented as a numbered series of steps (some with rather long discussions):

  1. To each element, such as a capacitor or inductor, attach a length of waveguide (electrical transmission line) having wave impedance $ R_0$, and make it infinitesimally long. (Take the limit as its length goes to zero.) A schematic depiction of this is shown in Fig.F.1a. For consistency, all signals are Laplace transforms of their respective time-domain signals. The length must approach zero in order not to introduce propagation delays into the signal path.

    Figure F.1: a) Physical schematic for the derivation of a wave digital model of driving-point impedance $ R(s)$. The inserted waveguide impedance $ R_0$ is real and positive, but otherwise arbitrary. b) Expanded view of the interior of the infinitesimal waveguide section, also representing the termination impedance $ R(s)$ as an impedance-step within the waveguide.
    \includegraphics{eps/wdelt}

    Points to note:

    • The infinitesimal waveguide is terminated by the element. The element reflects waves as if it were a new waveguide section at impedance $ R(s)$, as depicted in Fig.F.1b.

    • The interface to the element is recast as traveling-wave components $ F^{+}(s)$ and $ F^{-}(s)$ at impedance $ R_0$. In terms of these components, the physical force on the element is obtained by adding them together: $ F(s) = F^{+}(s)+F^{-}(s)$.

    • The waveguide impedance $ R_0$ is arbitrary because it has been physically introduced. We will need to know it when we connect this element to other elements. The element's interface to other elements is now a waveguide (transmission line) at real impedance $ R_0$.

    • The junction is ``parallel'' (cf. §7.2):

      • Force (voltage) must be continuous across the junction, since otherwise there would be a finite force across a zero mass, producing infinite acceleration.

      • The sum of velocities (currents) into the junction must be zero by conservation of mass (charge).


Reflectance of a General Lumped Waveguide Termination

Calculate the reflectance of the terminated waveguide. That is, find the Laplace transform of the return wave divided by the Laplace transform of the input wave going into the waveguide. In general, the reflectance of an impedance step for force waves (voltage waves in the electrical case) is

$\displaystyle \fbox{$\displaystyle \hat{\rho}_R(s) \isdef \frac{F^{-}(s)}{F^{+}(s)} = \frac{R(s)-R_0}{R(s)+R_0}$} \protect$ (F.1)

This is easily derived from continuity constraints across the junction. Specifically, referring to Fig.F.1b, let $ F_R(s) = F^{+}_R(s) + F^{-}_R(s)$ denote the physical force and its traveling-wave components within the ``pseudo-infinitesimal-generalized-waveguide'' defined by the element impedance $ R(s)$, with the `$ +$' superscript denoting a right-going wave.F.1 Similarly, let $ V(s) = V^{+}(s)+V^{-}(s)$ denote the velocity and its component wave variables on the side of the junction at impedance $ R_0$, and let $ V_R(s) = V^{+}_R(s)+V^{-}_R(s)$ denote the corresponding quantities on the element-side of the junction at impedance $ R(s)$. Again, the `$ +$' superscript denotes travel to the right. Then the physical continuity constraints imply

\begin{eqnarray*} F(s) &=& F_R(s)\\ 0 &=& V(s) + V_R(s) \end{eqnarray*}

By the definition of wave impedance in a waveguide, we have

\begin{eqnarray*} F^{+}(s) &=& \quad\! R_0 V^{+}(s)\\ F^{-}(s) &=& - R_0 V^{-}(... ... &=& \quad\! R(s) V^{+}_R(s)\\ F^{-}_R(s) &=& - R(s) V^{-}_R(s) \end{eqnarray*}

Thus,

\begin{eqnarray*} 0 &=& V(s) + V_R(s)\\ &=& \left[V^{+}(s)+V^{-}(s)\right] + ... ...s)}\right] &=& \frac{2}{R_0}F^{+}(s) + \frac{2}{R(s)}F^{+}_R(s) \end{eqnarray*}

Defining $ \Gamma_0 \isdef 1/R_0$ and $ \Gamma(s) \isdef 1/R(s)$, we have

$\displaystyle F(s) = \frac{2\Gamma_0}{\Gamma_0+\Gamma(s)} F^{+}(s) + \frac{2\Ga... ...a(s)} F^{+}_R(s) \isdef {\cal A}(s) F^{+}(s) + {\cal A}_R(s)F^{+}_R(s) \protect$ (F.2)

Now that we've solved for the junction force $ F(s)$, the outgoing waves are simply obtained from the force continuity constraint, $ F(s) = F^{+}(s)+F^{-}(s) = F^{+}_R(s)+F^{-}_R(s)$:
$\displaystyle F^{-}(s)$ $\displaystyle =$ $\displaystyle F(s) - F^{+}(s) \protect$ (F.3)
$\displaystyle F^{-}_R(s)$ $\displaystyle =$ $\displaystyle F(s) - F^{+}_R(s) \protect$ (F.4)

Finally, the force-wave reflectance of an impedance step from $ R_0$ to $ R(s)$ can be found by solving Eq.$ \,$(F.3) and (F.2) for $ F^{-}(s)/F^{+}(s)$ with $ F^{+}_R(s)$ set to zero:

\begin{eqnarray*} \hat{\rho}(s) &\isdef & \frac{F^{-}(s)}{F^{+}(s)} = \frac{F(s)... ...mma_0-\Gamma(s)}{\Gamma_0+\Gamma(s)} = \frac{R(s)-R_0}{R(s)+R_0} \end{eqnarray*}

as claimed.

Reflectances of Elementary Impedances

We now derive the reflectances of the elements used in LTI analog electric circuits, viz., the capacitor, inductor, and resistor.

Capacitor Reflectance

For a capacitor of $ C$ Farads, the driving-point impedance is (see §7.1.3)

$\displaystyle R_C(s)=\frac{1}{Cs} $

(or $ k/s$ for a spring with constant $ k$). Substituting into Eq.$ \,$(F.1) gives the reflectance

$\displaystyle \hat{\rho}_C(s) = \frac{R_C(s)-R_0}{R_C(s)+R_0} = \frac{1 - R_0 C s}{1 + R_0 C s} \protect$ (F.5)

Inductor Reflectance

For an inductor of $ L$ Henrys, we have

$\displaystyle R_L(s)$ $\displaystyle =$ $\displaystyle Ls$  
$\displaystyle \,\,\Rightarrow\,\,\hat{\rho}_L(s)$ $\displaystyle =$ $\displaystyle \frac{Ls-R_0}{Ls+R_0} = \frac{ s - R_0/L }{ s + R_0/L} \protect$ (F.6)

Resistor Reflectance

Finally, for a resistor of $ R$ Ohms, we get

$\displaystyle \hat{\rho}_R(s) = \frac{R-R_0}{R+R_0} = \frac{1 - R_0/R }{ 1 + R_0/R } \protect$ (F.7)

Note that both the capacitor and inductor reflectances are stable allpass filters, as they must be. Also, the resistor reflectance is always less than 1, no matter what waveguide impedance $ R_0>0$ we choose.

Choosing Impedance to Simplify Element Reflectance

Observe that there is a natural choice for each waveguide impedance which will give us a normalized, ``universal reflectance'' for each element:

  • For the capacitor, setting $ R_0 = 1/C$ gives

    $\displaystyle \fbox{$\displaystyle \hat{\rho}_C(s) = \frac{1 - s}{1 + s}$} \protect$ (F.8)

  • For the inductor, setting $ R_0=L$ gives

    $\displaystyle \fbox{$\displaystyle \hat{\rho}_L(s) = - \frac{1 - s}{1 + s}$} \protect$ (F.9)

  • And for the resistor, we set $ R_0 = R$ to obtain

    $\displaystyle \fbox{$\displaystyle \hat{\rho}_R(s) = 0$} \protect$ (F.10)

Digitizing Elementary Reflectances by Bilinear Transform

Going to discrete time via the bilinear transform means making the substitution

$\displaystyle s = c \frac{1-z^{-1}}{1+z^{-1}}$ (F.11)

where $ c>0$ is an arbitrary real constant, usually taken to be $ c=2/T$.

Solving for $ z^{-1}$ gives us the inverse bilinear transform:

$\displaystyle z^{-1}= \frac{1-s/c}{1+s/c} \protect$ (F.12)

In this case, we see that setting $ c=1$ further simplifies our universal reflectances in the digital domain:

Note that this choice of $ c$ is also the only one that eliminates delay-free paths in the fundamental elements. This allows them to be used as building blocks for explicit finite difference schemes.

We may still obtain the above results using the more typical value $ c=2/T$ (instead of $ c=1$) in the bilinear transform. From Eq.$ \,$(F.12), it is clear that changing $ c$ amounts to a linear frequency scaling of $ s=j\omega$. Such a scaling may be compensated by choosing the waveguide (port) impedances to be $ R_L = Lc = 2L/T$ (instead of $ R_L=L$) for the inductor, and $ R_C = T/(2C)$ (instead of $ 1/C$) for the capacitor.

Summary of Wave Digital Elements

From Eq.$ \,$(F.1), we have that the general reflectance of impedance $ R(s)$ with respect to the reference impedance $ R_0$ in the wave variable formulation is given by

$\displaystyle \fbox{$\displaystyle \hat{\rho}(s) \isdef \frac{R(s)-R_0}{R(s)+R_0}$} \protect$ (F.13)

In WDF construction, the free constant in the bilinear transform is taken to be $ c=1$. Thus we obtain $ \hat{\rho}_d(z) = \hat{\rho}[(1-z^{-1})/(1+z^{-1})]$. When $ R(s)$ is first order, it is possible to choose the reference impedance $ R_0$ so as to eliminate the delay-free path in the digital reflectance $ \hat{\rho}_d(z)$, and so its value depends on the actual physical element being digitized.

Wave Digital Mass

In the case of a mass $ m$, we have

$\displaystyle R(s) = ms $

which implies its reflectance is, from Eq.$ \,$(F.13),

$\displaystyle \hat{\rho}_m(s) = \frac{ms - R_0}{ms + R_0} $

Setting $ R_0= m$ gives

$\displaystyle \hat{\rho}_m(s) = \frac{s - 1 }{s + 1} $

and this choice also turns out to eliminate the delay-free path in the digital version. In view of the expression for the inverse bilinear transform in Eq.$ \,$(F.12), i.e., $ z=(1+s)/(1-s)$, the bilinear transform of $ \hat{\rho}(s)$ is immediately seen to be

$\displaystyle \fbox{$\displaystyle \hat{\tilde{\rho}}_m(z) = -z^{-1}$} $

where we defined $ \hat{\tilde{\rho}}_m(z) \isdef \hat{\rho}_m\left(\frac{z-1}{z+1}\right)$. The corresponding difference equation for the wave digital mass is

$\displaystyle f^{{-}}(n) = - f^{{+}}(n-1) $

and its wave flow diagram is drawn in Fig.F.2.

Figure F.2: Wave flow diagram for the wave digital mass. Note that the wave variables are written in the time domain as is customary in digital filter diagrams, while it would be more consistent (with the $ z^{-1}$ block) to keep them in the frequency domain as $ F^{+}(z)$ and $ F^{-}(z)$.
\includegraphics{eps/lWaveDigitalMass}

Thus, the wave digital mass is simply a unit-sample delay and a negation. The fact that the value of the mass has been canceled out will be addressed below in the subsection on ``adaptors,'' i.e., it only affects interconnection with other elements. For now, just remember that the reference impedance was chosen to be equal to the mass in order to get this simple wave flow diagram. Also note that the WDF mass simulator has no delay-free path from input to output.

Wave Digital Spring

In the case of a spring with stiffness $ k$, we have the impedance

$\displaystyle R(s) = k/s $

which gives the reflectance

$\displaystyle \hat{\rho}_k(s) = \frac{k/s - R_0}{k/s + R_0} $

As before, we may eliminate $ k$ by choosing $ R_0=k$ to get

$\displaystyle \hat{\rho}_k(s) = \frac{1 - s }{1 + s} = z^{-1} $

under the bilinear transform. So we have the digital reflectance

$\displaystyle \fbox{$\displaystyle \hat{\tilde{\rho}}_k(z) = z^{-1}$} \qquad\makebox[0pt][l]{(Wave Digital Spring)} $

and corresponding difference equation

$\displaystyle f^{{-}}(n) = f^{{+}}(n-1). $

Again the delay-free path has been eliminated. The wave flow diagram is shown in Fig.F.3.

Figure F.3: Wave flow diagram for the Wave Digital Spring.
\includegraphics{eps/lWaveDigitalSpring}

Thus, the WDF of a spring is simply a unit-sample delay, which is just the negative of the WDF mass. If we were to switch to velocity waves instead of force waves, both masses and springs would again correspond to unit-sample delays, but the spring would become inverting and the mass non-inverting.

Wave Digital Dashpot

Starting with a dashpot with coefficient $ \mu $, we have

$\displaystyle R(s) = \mu $

and reflectance

$\displaystyle \hat{\rho}_\mu(s) = \frac{\mu - R_0}{\mu + R_0} $

This time, choosing $ R_0$ equal to the element value gives

$\displaystyle \hat{\rho}_\mu(s) = 0 $

Conformally mapping the zero function yields the zero function so that

$\displaystyle \fbox{$\displaystyle \hat{\tilde{\rho}}_\mu(z) = 0$} $

as well. Thus, the WDF of a dashpot is a ``wave sink,'' as diagrammed in Fig.F.4.

Figure F.4: Wave flow diagram for the Wave Digital Dashpot.
\includegraphics{eps/lWaveDigitalDashpot}

In the context of waveguide theory, a zero reflectance corresponds to a matched impedance, i.e., the terminating transmission-line impedance equals the characteristic impedance of the line.

The difference equation for the wave digital dashpot is simply $ f^{{-}}(n)=0$. While this may appear overly degenerate at first, remember that the interface to the element is a port at impedance $ R_0=\mu$. Thus, in this particular case only, the infinitesimal waveguide interface is the element itself.

Limiting Cases

The force-wave reflectance of an infinite impedance (rigid wall or ``open circuit'') is

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

Similarly, the force-wave reflectance of a zero impedance (free termination, frictionless surface, or ``short circuit'') is

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

For velocity waves, we obtain the opposite results: rigid terminations are inverting, and free terminations are non-inverting.

Unit Elements

The unit element two-port is simply a bidirectional delay line with half a sample delay in each direction. As a result, it really belongs under the topic of distributed modeling. To avoid delay-free loops, Fettweis noted [135] that every pair of adaptors must be separated by at least one unit element. More recently, this objective is accomplished instead using ``reflection-free ports'' [136] (see also §F.2.2).

Adaptors for Wave Digital Elements

An adaptor is an $ N$-port memoryless interface which interconnects wave digital elements. Since each element's ``port'' is a connection to an infinitesimal waveguide section at some real wave impedance $ R_i$, and since the input/output signals are wave variables (traveling-waves within the waveguide), the adaptor must implement signal scattering appropriate for the connection of such waveguides. In other words, an $ N$-port adaptor in a wave digital filter performs exactly the same computation as an $ N$-port scattering junction in a digital waveguide network.F.2

This section first addresses the simpler two-port case, followed by a derivation of the general $ N$-port adaptor, for both parallel and series connections of wave digital elements.

As discussed in §7.2, a physical connection of two or more ports can either be in parallel (forces are equal and the velocities sum to zero) or in series (velocities equal and forces sum to zero). Combinations of parallel and series connections are also of course possible.

Two-Port Parallel Adaptor for Force Waves

Figure F.5a illustrates a generic parallel two-port connection in terms of forces and velocities.

Figure: a) Two-port description of the adaptor implementing a parallel connection between reference impedances $ R_1$ and $ R_2$. b) Corresponding parallel force scattering junction (adaptor wave flow diagram) in Kelly-Lochbaum form. Compare with Fig.F.7.
\includegraphics[width=\twidth]{eps/lAdaptorParallel}

As discussed in §7.2, a parallel connection is characterized by a common force and velocities which sum to zero:

\begin{eqnarray*} && f_1(n) = f_2(n) \isdef f_J(n)\\ && v_1(n) + v_2(n) = 0 \end{eqnarray*}

Following the same derivation leading to Eq.$ \,$(F.2), and defining $ \Gamma _i=1/R_i$ for notational convenience, we obtain

\begin{eqnarray*} 0 &=& v_1+v_2 \\ &=& \frac{f^{{+}}_1-f^{{-}}_1}{R_1} + \frac... ...amma _1 f^{{+}}_1 + \Gamma _2 f^{{+}}_2 }{\Gamma _1+\Gamma _2} . \end{eqnarray*}

The outgoing wave variables are given by

\begin{eqnarray*} f^{{-}}_1(n) &=& f_J(n) - f^{{+}}_1(n) \\ f^{{-}}_2(n) &=& f_J(n) - f^{{+}}_2(n) \end{eqnarray*}

Defining the reflection coefficient as

$\displaystyle \rho \isdef \frac{R_2-R_1}{R_2+R_1} $

we have that the scattering relations for the two-port parallel adaptor are

\begin{eqnarray*} f^{{-}}_1 &=& \rho f^{{+}}_1 + (1-\rho) f^{{+}}_2 \protect \\ f^{{-}}_2 &=& (1+\rho)f^{{+}}_1 - \rho f^{{+}}_2 \protect \end{eqnarray*}

as diagrammed in Fig.F.5b. This can be called the Kelly-Lochbaum implementation of the two-port force-wave adaptor.

Now that we have a proper scattering interface between two reference impedances, we may connect two wave digital elements together, setting $ R_1$ to the port impedance of element 1, and $ R_2$ to the port impedance of element 2. An example is shown in Fig.F.35.

The Kelly-Lochbaum adaptor in Fig.F.5b evidently requires four multiplies and two additions. Note that we can factor out the reflection coefficient in each equation to obtain

\begin{eqnarray*} f^{{-}}_1 &=& f^{{+}}_2 + \rho(f^{{+}}_1 - f^{{+}}_2)\\ f^{{-}}_2 &=& f^{{+}}_1 + \rho(f^{{+}}_1 - f^{{+}}_2) \end{eqnarray*}

which requires only one multiplication and three additions. This can be called the one-multiply form. The one-multiply form is most efficient in custom VLSI. The Kelly-Lochbaum form, on the other hand, may be more efficient in software, and slightly faster (by one addition) in parallel hardware.

Compatible Port Connections

Note carefully that to connect a wave digital element to port $ i$ of the adaptor, we route the signal $ f^{{-}}(n)$ coming out of the element to become $ f^{{+}}_i(n)$ on the adaptor port, and the signal $ f^{{-}}_i(n)$ coming out of port $ i$ of the adaptor goes into the element as $ f^{{+}}(n)$. Such a connection is said to be a compatible port connection. In other words, the connections must be made such that the arrows go in the same direction in the wave flow diagram.

General Parallel Adaptor for Force Waves

In the more general case of $ N$ wave digital element ports being connected in parallel, we have the physical constraints

    $\displaystyle f_1(n) = f_2(n) = \cdots = f_N(n) \isdef f_J(n)$ (F.14)
    $\displaystyle v_1(n) + v_2(n) + \cdots + v_N(n) = 0$ (F.15)

The derivation for the two-port case extends to the $ N$-port case without modification:
0 $\displaystyle =$ $\displaystyle \sum_{i=1}^N v_i$  
  $\displaystyle =$ $\displaystyle \sum_{i=1}^N\frac{f^{{+}}_i-f^{{-}}_i}{R_i}$  
  $\displaystyle =$ $\displaystyle \sum_{i=1}^N\frac{2f^{{+}}_i-f_J}{R_i}$  
  $\displaystyle \isdef$ $\displaystyle \sum_{i=1}^N \left(2\Gamma _if^{{+}}_i-\Gamma _i f_J \right)$  
$\displaystyle \,\,\Rightarrow\,\, \sum_{j=1}^N \Gamma _j f_J$ $\displaystyle =$ $\displaystyle \sum_{i=1}^N 2\Gamma _i f^{{+}}_i$  
$\displaystyle \,\,\Rightarrow\,\, f_J$ $\displaystyle =$ $\displaystyle \frac{\sum_{i=1}^N 2\Gamma _i f^{{+}}_i}{\sum_{j=1}^N \Gamma _j} . \protect$ (F.16)

The outgoing wave variables are given by

$\displaystyle f^{{-}}_i(n) = f_J(n) - f^{{+}}_i(n) $

Alpha Parameters

It is customary in the wave digital filter literature to define the alpha parameters as

$\displaystyle \fbox{$\displaystyle \alpha_i \isdef \frac{2\Gamma _i}{\sum_{j=1}^N \Gamma _j}$} \protect$ (F.17)

where $ \Gamma _i \isdef 1/R_i$ are the admittances of the wave digital element interfaces (or ``reference admittances,'' in WDF terminology). In terms of the alpha parameters, the force-wave parallel adaptor performs the following computations:
$\displaystyle f_J(n)$ $\displaystyle =$ $\displaystyle \sum_{i=1}^N \alpha_i f^{{+}}_i(n)\protect$ (F.18)
$\displaystyle f^{{-}}_i(n)$ $\displaystyle =$ $\displaystyle f_J(n) - f^{{+}}_i(n)\protect$ (F.19)

We see that $ N$ multiplies and $ 2N-1$ additions are required. However, by observing from Eq.$ \,$(F.17) that

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

we may implement $ \alpha_1$, say, as $ 2-\sum_{i=2}^N\alpha_i$ in order to eliminate one multiply. In WDF terminology, port 1 is then a dependent port.

Reflection Coefficient, Parallel Case

The reflection coefficient seen at port $ i$ is defined as

$\displaystyle \rho_i \isdef \left. \frac{f^{{-}}_i(n)}{f^{{+}}_i(n)} \right\vert _{f^{{+}}_j(n)=0, \forall j\neq i} \protect$ (F.20)

In other words, the reflection coefficient specifies what portion of the incoming wave $ f^{{+}}_i(n)$ is reflected back to port $ i$ as part of the outgoing wave $ f^{{-}}_i(n)$. The total outgoing wave on port $ i$ is the superposition of the reflected wave and the $ N-1$ transmitted waves from the other ports:

$\displaystyle f^{{-}}_i(n) = \rho_i f^{{+}}_i + \sum_{j\neq i} \tau_{ji} f^{{+}}_j \protect$ (F.21)

where $ \tau_{ji}$ denotes the transmission coefficient from port $ j$ to port $ i$. Starting with Eq.$ \,$(F.19) and substituting Eq.$ \,$(F.18) gives

\begin{eqnarray*} f^{{-}}_i(n) &=& f_J(n) - f^{{+}}_i(n)\\ &=& \left(\sum_{j=1... ...\alpha_i - 1)f^{{+}}_i(n) + \sum_{j\neq i} \alpha_j f^{{+}}_j(n) \end{eqnarray*}

Equating like terms with Eq.$ \,$(F.21), we obtain

$\displaystyle \rho_i$ $\displaystyle =$ $\displaystyle \alpha_i - 1 \protect$ (F.22)
$\displaystyle \tau_{ji}$ $\displaystyle =$ $\displaystyle \alpha_j, \quad (i\neq j) \protect$ (F.23)

Thus, the $ j$th alpha parameter is the force transmission coefficient from $ j$th port to any other port (besides the $ i$th). To convert the transmission coefficient from the $ i$th port to the reflection coefficient for that port, we simply subtract 1. This general relationship is specific to force waves at a parallel junction, as we will soon see.

Physical Derivation of Reflection Coefficient

Physically, the reflection coefficient seen at port $ i$ is due to an impedance step from $ R_i$, that of the port interface, to a new impedance consisting of the parallel combination of all other port impedances meeting at the junction. Let

$\displaystyle \Gamma _J(i) \isdef \sum_{i\neq j} \Gamma _i \protect$ (F.24)

denote this parallel combination, in admittance form. Then we must have

$\displaystyle \rho_i = \frac{R_J(i)-R_i}{R_J(i)+R_i} = \frac{\Gamma _i-\Gamma _J(i)}{\Gamma _i+\Gamma _J(i)} \protect$ (F.25)

Let's check this ``physical'' derivation against the formal definition Eq.$ \,$(F.20) leading to $ \rho_i = \alpha_i - 1$ in Eq.$ \,$(F.22). Toward this goal, let

$\displaystyle \Gamma _J \isdef \sum_{j=1}^N \Gamma _j $

denote the parallel combination of all admittances connected to the junction. Then by Eq.$ \,$(F.24), we have $ \Gamma _J = \Gamma _i + \Gamma _J(i)$ for all $ i$. Now, from Eq.$ \,$(F.17),

\begin{eqnarray*} \rho_i &\isdef & \alpha_i - 1 \;\isdef \; \frac{2\Gamma _i}{\... ..._i + \Gamma _J(i)} \;=\; \frac{R_J(i) - R_i}{\Gamma _J(i)-R_i} \end{eqnarray*}

and the result is verified.

Reflection Free Port

It is useful in practice, such as when connecting two adaptors together, to make one port reflection free. A reflection-free port is defined to have a zero reflection coefficient. For port $ i$ of a parallel adaptor to be reflection free, we must have, from Eq.$ \,$(F.25),

$\displaystyle R_i = R_J(i) \isdef \frac{1}{\sum_{i\neq j} \Gamma _i} $

Thus, the port's impedance must equal the parallel combination of the other port impedances at the junction. In this case, the junction as a whole ``perfectly terminates'' the reflection free port, so no reflections come back from it.

Connecting two adaptors at a reflection-free port prevents the formation of a delay-free loop which would otherwise occur [136]. As a result, multi-port junctions can be joined without having to insert unit elements (see §F.1.7) to avoid creating delay-free loops. Only one of the two ports participating in the connection needs to be reflection free.

We can always make a reflection-free port at the connection of two adaptors because the ports used for this connection (one on each adaptor) were created only for purposes of this connection. They can be set to any impedance, and only one of them needs to be reflection free.

To interconnect three adaptors, labeled $ A$, $ B$, and $ C$, we may proceed as follows: Let $ A$ be augmented with two unconstrained ports, having impedances $ R_1$ and $ R_2$. Add a reflection-free port to $ B$, and suppose its impedance has to be $ R_B$. Add a reflection-free port to $ C$, and suppose its impedance has to be $ R_C$. Now set $ R_1=R_B$ and connect $ B$ to $ A$ via the corresponding ports. Similarly, set $ R_2=R_C$ and connect $ C$ to $ A$ accordingly. This adaptor-connection protocol clearly extends to any number of adaptors.

Two-Port Series Adaptor for Force Waves

Figure F.6a illustrates a generic two-port description of the series adaptor.

Figure F.6: a) Two-port description of the adaptor implementing a series connection between reference impedances $ R_1$ and $ R_2$. b) Corresponding series force scattering junction (adaptor wave flow diagram) in Kelly-Lochbaum form.
\includegraphics[width=\twidth]{eps/lAdaptorSeries}

As discussed in §7.2, a series connection is characterized by a common velocity and forces which sum to zero at the junction:

\begin{eqnarray*} && f_1(n) + f_2(n) = 0\\ && v_1(n) = v_2(n) \isdef v_J(n) \end{eqnarray*}

The derivation can proceed exactly as for the parallel junction in §F.2.1, but with force and velocity interchanged, i.e., $ f\leftrightarrow v$, and with impedance and admittance interchanged, i.e., $ R\leftrightarrow \Gamma $. In this way, we may take the dual of Eq.$ \,$(F.14) to get

\begin{eqnarray*} v^{-}_1 &=& -\rho v^{+}_1 + (1+\rho) v^{+}_2\\ v^{-}_2 &=& (1-\rho)v^{+}_1 + \rho v^{+}_2 \end{eqnarray*}

diagrammed in Fig.F.7. Converting back to force wave variables via $ f^{{+}}_i=R_iv^{+}_i$ and $ f^{{-}}_i=-R_iv^{-}_i$, and noting that $ (1+\rho)R_1/R_2 = 1-\rho$, we obtain, finally,

\begin{eqnarray*} f^{{-}}_1 &=& \rho f^{{+}}_1 - (1-\rho) f^{{+}}_2\\ f^{{-}}_2 &=& -(1+\rho)f^{{+}}_1 - \rho f^{{+}}_2 \end{eqnarray*}

as diagrammed in Fig.F.6b. The one-multiply form is now

\begin{eqnarray*} f^{{-}}_1 &=& -f^{{+}}_2 + \rho(f^{{+}}_1 + f^{{+}}_2)\\ f^{{-}}_2 &=& -f^{{+}}_1 - \rho(f^{{+}}_1 + f^{{+}}_2). \end{eqnarray*}

Figure F.7: Series velocity scattering junction in Kelly-Lochbaum form.
\includegraphics[scale=0.9]{eps/lscat_vel_series_renum}

General Series Adaptor for Force Waves

In the more general case of $ N$ ports being connected in series, we have the physical constraints

\begin{eqnarray*} && v_1(n) = v_2(n) = \cdots = v_N(n) \isdef v_J(n)\\ && f_1(n) + f_2(n) + \cdots + f_N(n) = 0 \end{eqnarray*}

The derivation is the dual of that in the parallel case (cf. Eq.$ \,$(F.16)), i.e., force and velocity are interchanged, and impedance and admittance are interchanged:

\begin{eqnarray*} 0 &=& \sum_{i=1}^N f_i \\ &=& \sum_{i=1}^NR_i\left(v^{+}_i-v... ...uad v_J &=& \frac{\sum_{i=1}^N 2R_i v^{+}_i}{\sum_{j=1}^N R_j} . \end{eqnarray*}

The outgoing wave variables are given by

$\displaystyle v^{-}_i(n) = v_J(n) - v^{+}_i(n) $

Beta Parameters

It is customary in the wave digital filter literature to define the beta parameters as

$\displaystyle \fbox{$\displaystyle \beta_i \isdef \frac{2R_i}{\sum_{j=1}^N R_j}$} \protect$ (F.26)

where $ R_i$ are the port impedances (attached element reference impedances). In terms of the beta parameters, the force-wave series adaptor performs the following computations:
$\displaystyle v_J(n)$ $\displaystyle =$ $\displaystyle \sum_{i=1}^N \beta_i v^{+}_i(n)\protect$ (F.27)
$\displaystyle v^{-}_i(n)$ $\displaystyle =$ $\displaystyle v_J(n) - v^{+}_i(n)\protect$ (F.28)

However, we normally employ a mixture of parallel and series adaptors, while keeping a force-wave simulation. Since $ f^{{+}}_i(n) = R_i v^{+}_i(n)$, we obtain, after a small amount of algebra, the following recipe for the series force-wave adaptor:

$\displaystyle f^{{+}}_J(n)$ $\displaystyle =$ $\displaystyle \sum_{i=1}^N f^{{+}}_i(n)\protect$ (F.29)
$\displaystyle f^{{-}}_i(n)$ $\displaystyle =$ $\displaystyle f^{{+}}_i(n) - \beta_if^{{+}}_J(n)\protect$ (F.30)

We see that we have $ N$ multiplies and $ 2N-1$ additions as in the parallel-adaptor case. However, we again have from Eq.$ \,$(F.26) that

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

so that we may implement one beta parameter as 2 minus the sum of the rest, thus eliminating a multiplication by creating a dependent port.

Reflection Coefficient, Series Case

The velocity reflection coefficient seen at port $ i$ is defined as

$\displaystyle \rho^v_i \isdef \left. \frac{v^{-}_i(n)}{v^{+}_i(n)} \right\vert _{v^{+}_j(n)=0, \forall j\neq i} \protect$ (F.31)

Representing the outgoing velocity wave $ v^{-}_i(n)$ as the superposition of the reflected wave $ \rho^v_iv^{+}_i(n)$ plus the $ N-1$ transmitted waves from the other ports, we have

$\displaystyle v^{-}_i(n) = \rho^v_i v^{+}_i + \sum_{j\neq i} \tau^v_{ji} v^{+}_j \protect$ (F.32)

where $ \tau^v_{ji}$ denotes the velocity transmission coefficientvelocity!transmission coefficient from port $ j$ to port $ i$. Substituting Eq.$ \,$(F.29) into Eq.$ \,$(F.30) yields

\begin{eqnarray*} v^{-}_i(n) &=& v_J(n) - v^{+}_i(n)\\ &=& \left(\sum_{j=1}^N ... ... &=& (\beta_i - 1)v^{+}_i(n) + \sum_{j\neq i} \beta_j v^{+}_j(n) \end{eqnarray*}

Equating like terms with Eq.$ \,$(F.32) gives

$\displaystyle \rho^v_i$ $\displaystyle =$ $\displaystyle \beta_i - 1 \protect$ (F.33)
$\displaystyle \tau^v_{ji}$ $\displaystyle =$ $\displaystyle \beta_j, \quad (i\neq j)$ (F.34)

Thus, the $ j$th beta parameter is the velocity transmission coefficient from $ j$th port to any other port (besides the $ i$th). To convert the transmission coefficient from the $ i$th port to the reflection coefficient for that port, we simply subtract 1. These relationships are specific to velocity waves at a series junction (cf. Eq.$ \,$(F.22)). They are exactly the dual of Equations (F.22-F.23) for force waves at a parallel junction.

Physical Derivation of Series Reflection Coefficient

Physically, the force-wave reflection coefficient seen at port $ i$ of a series adaptor is due to an impedance step from $ R_i$, that of the port interface, to a new impedance consisting of the series combination of all other port impedances meeting at the junction. Let

$\displaystyle R_J(i) \isdef \sum_{i\neq j} R_i \protect$ (F.35)

denote this series combination. Then we must have, as in Eq.$ \,$(F.25),

$\displaystyle \rho_i = \frac{R_J(i)-R_i}{R_J(i)+R_i}$ (F.36)

Let's check this ``physical'' derivation against the formal definition Eq.$ \,$(F.31) leading to $ \rho^v_i = \beta_i - 1$ in Eq.$ \,$(F.33). Define the total junction impedance as

$\displaystyle R_J \isdef \sum_{j=1}^N R_j $

This is the series combination of all impedances connected to the junction. Then by Eq.$ \,$(F.35), $ R_J = R_i + R_J(i)$ for all $ i$. From Eq.$ \,$(F.26), the velocity reflection coefficient is given by

\begin{eqnarray*} \rho^v_i &\isdef & \beta_i - 1 \;\isdef \; \frac{2R_i}{R_J} -... ..._J(i)}\\ &=& \frac{R_i - R_J(i)}{R_i + R_J(i)}\\ &=& -\rho_i \end{eqnarray*}

Since

$\displaystyle \rho^v_i\isdef \frac{v^{-}_i(n)}{v^{+}_i(n)} = \frac{-f^{{-}}_i(n)/R_i}{f^{{+}}_i(n)/R_i} = - \frac{f^{{-}}_i(n)}{f^{{+}}_i(n)} \isdef -\rho_i $

the result follows.

Series Reflection Free Port

For port $ i$ to be reflection free in a series adaptor, we require

$\displaystyle R_i = R_J(i) \isdef \sum_{i\neq j} R_i \protect$ (F.37)

That is, the port's impedance must equal the series combination of the other port impedances at the junction. This result can be compared with that for the parallel junction in §F.2.2.

The series adaptor has now been derived in a way which emphasizes its duality with respect to the parallel adaptor.

Wave Digital Modeling Examples

This section presents a series of examples, working up from very simple to more complicated situations.

``Piano hammer in flight''

Suppose we wish to model a situation in which a mass of size $ m$ kilograms is traveling with a constant velocity. This is an appropriate model for a piano hammer after its key has been pressed and before the hammer has reached the string.

Figure F.2 shows the ``wave digital mass'' derived previously. The derivation consisted of inserting an infinitesimal waveguideF.3 having (real) impedance $ m$, solving for the force-wave reflectance of the mass as seen from the waveguide, and then mapping it to the discrete time domain using the bilinear transform.

We now need to attach the other end of the transmission line to a ``force source'' which applies a force of zero newtons to the mass. In other words, we need to terminate the line in a way that corresponds to zero force.

Let the force-wave components entering and leaving the mass be denoted $ f^{{+}}$ and $ f^{{-}}$, respectively (i.e., we are dropping the subscript `d' in Fig.F.2). The physical force associated with the mass is

$\displaystyle f(n) = f^{{+}}(n) + f^{{-}}(n) = f^{{+}}(n) - f^{{+}}(n-1) $

The zero-force case is therefore obtained when $ f^{{+}}(n) = -f^{{-}}(n) = f^{{+}}(n-1)$. This is illustrated in Fig.F.8.

Figure F.8: Wave digital mass in flight at a constant velocity.
\includegraphics{eps/wdhammer}

Figure F.8a (left portion) illustrates what we derived by physical reasoning, and as such, it is most appropriate as a physical model of the constant-velocity mass. However, for actual implementation, Fig.F.8b would be more typical in practice. This is because we can always negate the state variable $ x(n)$ if needed to convert it from $ f^{{+}}(n-1)$ to $ f^{{-}}(n)$. It is very common to see final simplifications like this to maximize efficiency.

Note that Fig.F.8b can be interpreted physically as a wave digital spring displaced by a constant force $ f(n) = 2x(n)$.

Extracting Physical Quantities

Since we are using a force-wave simulation, the state variable $ x(n)$ (delay element output) is in units of physical force (newtons). Specifically, $ x(n) = f^{{+}}(n-1)$. (The physical force is, of course, 0, while its traveling-wave components are not 0 unless the mass is at rest.) Using the fundamental relations relating traveling force and velocity waves

\begin{eqnarray*} f^{{+}}(n) &\isdef & \quad\! R_0 v^{+}(n)\\ f^{{-}}(n) &\isdef & - R_0 v^{-}(n)\\ \end{eqnarray*}

where $ R_0= m$ here, it is easy to convert the state variable $ x(n)$ to other physical units, as we now demonstrate.

The velocity of the mass, for example, is given by

$\displaystyle v(n) = v^{+}(n) + v^{-}(n) = \frac{f^{{+}}(n)}{m} - \frac{f^{{-}}(n)}{m} = \frac{2f^{{+}}(n)}{m} = \frac{2}{m}x(n) $

Thus, the state variable $ x(n)$ can be scaled by $ 2/m$ to ``read out'' the mass velocity.

The kinetic energy of the mass is given by

$\displaystyle {\cal E}_m = \frac{1}{2}mv^2(n) = \frac{2}{m}x^2(n) $

I.e., the square of the state variable $ x(n)$ can be scaled by $ 2/m$ to produce the physical kinetic energy associated with the mass.

Force Driving a Mass

Suppose now that we wish to drive the mass along a frictionless surface using a variable force $ f(n)$. This is similar to the previous example, except that we now want the traveling-wave components of the force on the mass to sum to $ f(n)$ instead of 0:

$\displaystyle f(n) = f^{{+}}(n) + f^{{-}}(n) $

Since $ f(n)$ and $ f^{{-}}(n)$ are given, $ f^{{+}}(n)$ must be computed as $ f^{{+}}(n) = f(n) - f^{{-}}(n)$. This is shown in Fig.F.9.

Figure F.9: Wave digital mass driven by external force $ f(n)$.
\includegraphics{eps/wdhf}

The simplified form in Fig.F.9b can be interpreted as a wave digital spring with applied force $ f(n)$ delivered from an infinite source impedance. That is, when the applied force goes to zero, the termination remains rigid at the current displacement.

A More Formal Derivation of the Wave Digital Force-Driven Mass

Above we derived how to handle the external force by direct physical reasoning. In this section, we'll derive it using a more general step-by-step procedure which can be applied systematically to more complicated situations.

Figure F.10 gives the physical picture of a free mass driven by an external force in one dimension. Figure F.11 shows the electrical equivalent circuit for this scenario in which the external force is represented by a voltage source emitting $ f(t)$ volts, and the mass is modeled by an inductor having the value $ L=m$ Henrys.

Figure F.10: Physical diagram of an external force driving a mass sliding on a frictionless surface.
\includegraphics{eps/forcemass}

Figure: Electrical equivalent circuit of the force-driven mass in Fig.F.10.
\includegraphics{eps/forcemassec}

The next step is to convert the voltages and currents in the electrical equivalent circuit to wave variables. Figure F.12 gives an intermediate equivalent circuit in which an infinitesimal transmission line section with real impedance $ R_0= m$ has been inserted to facilitate the computation of the wave-variable reflectance, as we did in §F.1.1 to derive Eq.$ \,$(F.1).

Figure F.12: Intermediate equivalent circuit for the force-driven mass in which an infinitesimal transmission line section has been inserted to facilitate conversion of the mass impedance $ ms$ into a wave-variable reflectance.
\includegraphics{eps/forcemassscat}

Figure: Intermediate wave-variable model of the force-driven mass of Fig.F.11.
\includegraphics{eps/forcemassdt}

Figure F.13 depicts a next intermediate equivalent circuit in which the mass has been replaced by its reflectance (using ``$ z^{-1}$'' to denote the continuous-time reflectance $ (1-s)/(1+s)$, as derived in §F.1.1). The infinitesimal transmission-line section is now represented by a ``resistor'' since, when the voltage source is initially ``switched on'', it only ``sees'' a real resistance having the value $ m$ Ohms (the waveguide interface). After a short period of time determined by the reflectance of the mass,F.4 ``return waves'' from the mass result in an ultimately reactive impedance. This of course must be the case because the mass does not dissipate energy. Therefore, the ``resistor'' of $ m$ Ohms is not a resistor in the usual sense since it does not convert potential energy (the voltage drop across it) into heat. Instead, it converts potential energy into propagating waves with 100% efficiency. Since all of this wave energy is ultimately reflected by the terminating element (mass, spring, or any combination of masses and springs), the net effect is a purely reactive impedance, as we know it must be.

Figure F.14: Interconnection of the wave digital mass with an ideal force source by means of a two-port parallel adaptor. The symbol ``$ \vert\vert$'' is used in the WDF literature to signify a parallel adaptor.
\includegraphics{eps/forcemassjunc}

To complete the wave digital model, we need to connect our wave digital mass to an ideal force source which asserts the value $ f(n)$ each sample time. Since an ideal force source has a zero internal impedance, we desire a parallel two-port junction which connects the impedances $ R_1 = 0$ ( $ \Gamma _1=\infty$) and $ R_2=m$ ( $ \Gamma _2=1/m$), as shown in Fig.F.14. From Eq.$ \,$(F.18) we have that the common junction force is equal to

\begin{eqnarray*} f(n) &=& \alpha_1 f^{{+}}_1(n) + \alpha_2 f^{{+}}_2(n)\\ &=&... ...c{2}{m}}{\infty+\frac{1}{m}} f^{{+}}_2(n)\\ &=& 2 f^{{+}}_1(n) \end{eqnarray*}

from which we conclude that

$\displaystyle \fbox{$\displaystyle f^{{+}}_1(n) = \frac{f(n)}{2}$} $

The outgoing waves are, by Eq.$ \,$(F.19),

\begin{eqnarray*} f^{{-}}_1(n) &=& f(n) - f^{{+}}_1(n) = \frac{f(n)}{2}\\ f^{{-}}_2(n) &=& f(n) - f^{{+}}_2(n) = f(n) - f^{{-}}_m(n) \end{eqnarray*}

Since $ \alpha_1=1$ and $ \alpha_2=0$ for this model, the reflection coefficient seen on port 1 is $ \rho = \alpha_1 - 1 = 1$. The transmission coefficient from port 1 is $ 1+\rho=2$. In the opposite direction, the reflection coefficient on port 2 is $ -\rho= -1$, and the transmission coefficient from port 2 is $ 1-\rho = 0$. The final result, drawn in Kelly-Lochbaum form (see §F.2.1), is diagrammed in Fig.F.15, as well as the result of some elementary simplifications. The final model is the same as in Fig.F.9, as it should be.

Figure F.15: Wave digital mass driven by external force $ f(n)$.
\includegraphics{eps/forcemasswdf}

Force Driving a Spring against a Wall

For this example, we have an external force $ f(n)$ driving a spring $ k$ which is terminated on the other end at a rigid wall. Figure F.16 shows the physical diagram and the electrical equivalent circuit is given in Fig.F.17.

Figure F.16: External force driving a spring terminated by a rigid wall.
\includegraphics{eps/forcespring}

Figure: Electrical equivalent circuit of the compressed spring of Fig.F.16.
\includegraphics{eps/forcespringec}

Figure F.18 depicts the insertion of an infinitesimal transmission line, and Fig.F.19 shows the result of converting the spring impedance to wave variable form.

Figure F.18: Intermediate equivalent circuit for the force-driven spring in which an infinitesimal transmission line section has been inserted to facilitate conversion of the spring impedance $ k/s$ into a wave-variable reflectance.
\includegraphics{eps/forcespringscat}

Figure: Intermediate wave-variable model of Fig.F.17.
\includegraphics{eps/forcespringdt}

The two-port adaptor needed for this problem is the same as that for the force-driven mass, and the final result is shown in Fig.F.20.

Figure F.20: Wave digital spring driven by external force $ f(n)$.
\includegraphics{eps/forcespringwdf}

Note that the spring model is being driven by a force from a zero source impedance, in contrast with the infinite source impedance interpretation of Fig.F.8b as a compressed spring. In this case, if the driving force goes to zero, the spring force goes immediately to zero (``free termination'') rather than remaining fixed.

Spring and Free Mass

For this example, we have an external force $ f(n)$ driving a spring $ k$ which in turn drives a free mass $ m$. Since the force on the spring and the mass are always the same, they are formally ``parallel'' impedances.

This problem is easier than it may first appear since an ideal ``force source'' (i.e., one with a zero source impedance) driving impedances in parallel can be analyzed separately for each parallel branch. That is, the system is equivalent to one in which the mass and spring are not connected at all, and each has its own copy of the force source. With this insight in mind, one can immediately write down the final wave-digital model shown in Fig.F.25. However, we will go ahead and analyze this case more formally since it has some interesting features.

Figure F.21 shows the physical diagram of the spring-mass system driven by an external force. The electrical equivalent circuit appears in Fig.F.22, and the first stage of a wave-variable conversion is shown in Fig.F.23.

Figure F.21: External force driving a spring which in turn drives a free mass sliding on a frictionless surface.
\includegraphics{eps/springmass}

Figure: Electrical equivalent circuit of the spring/mass system of Fig.F.21.
\includegraphics{eps/springmassec}

Figure: Intermediate wave-variable model of the mass and dashpot of Fig.F.22.
\includegraphics{eps/springmassdt}

Figure F.24: Wave digital model for the parallel combination of a wave digital mass $ m$, wave digital spring $ k$, and driving force $ f(n)$. The corresponding port impedances are $ m$, $ k$, and 0, respectively.
\includegraphics{eps/springmassjunc}

For this example we need a three-port parallel adaptor, as shown in Fig.F.24 (along with its attached mass and spring). The port impedances are 0, $ k$, and $ m$, yielding alpha parameters $ \alpha_1 =2$ and $ \alpha_2=\alpha_3=0$. The final result, after the same sorts of elementary simplifications as in the previous example, is shown in Fig.F.25. As predicted, a force source driving elements in parallel is equivalent to a set of independent force-driven elements.

Figure F.25: Wave digital filter model of an external force driving a mass through a spring. The mass force-wave components are denoted $ f^\pm _m$, while those for the spring are denoted $ f^\pm _k$.
\includegraphics{eps/wdspringmass}

From this and the preceding example, we can see a general pattern: Whenever there is an ideal force source driving a parallel junction, then $ \Gamma _1=\infty$ and all other port admittances are finite. In this case, we always obtain $ \alpha_1 =2$ and $ \alpha_i=0$, $ i\neq 1$.

Mass and Dashpot in Series

This is our first example illustrating a series connection of wave digital elements. Figure F.26 gives the physical scenario of a simple mass-dashpot system, and Fig.F.27 shows the equivalent circuit. Replacing element voltages and currents in the equivalent circuit by wave variables in an infinitesimal waveguides produces Fig.F.28.

Figure F.26: External force driving a mass which in turn drives a dashpot terminated on the other end by a rigid wall.
\includegraphics{eps/massdash}

Figure: Electrical equivalent circuit of the mass and dashpot system of Fig.F.26.
\includegraphics{eps/massdashec}

Figure: Intermediate wave-variable model of the mass and dashpot of Fig.F.27.
\includegraphics{eps/massdashdt}

Figure F.29: Wave digital filter for an ideal force source in parallel with the series combination of a mass $ m$ and dashpot $ \mu $. The parallel and series adaptors are joined at an impedance $ R$ which is calculated to suppress reflection from port 1 of the series adaptor.
\includegraphics{eps/massdashjunc}

The system can be described as an ideal force source $ f(t)$ connected in parallel with the series connection of mass $ m$ and dashpot $ \mu $. Figure F.29 illustrates the resulting wave digital filter. Note that the ports are now numbered for reference. Two more symbols are introduced in this figure: (1) the horizontal line with a dot in the middle indicating a series adaptor, and (2) the indication of a reflection-free port on input 1 of the series adaptor (signal $ f^{{+}}_1(n)$). Recall that a reflection-free port is always necessary when connecting two adaptors together, to avoid creating a delay-free loop.

Let's first calculate the impedance $ R$ necessary to make input 1 of the series adaptor reflection free. From Eq.$ \,$(F.37), we require

$\displaystyle R = m + \mu $

That is, the impedance of the reflection-free port must equal the series combination of all other port impedances meeting at the junction.

The parallel adaptor, viewed alone, is equivalent to a force source driving impedance $ R=m+\mu$. It is therefore realizable as in Fig.F.20 with the wave digital spring replaced by the mass-dashpot assembly in Fig.F.29. However, we can also carry out a quick analysis to verify this: The alpha parameters are

\begin{eqnarray*} \alpha_1 &\isdef & \frac{2\Gamma _1}{\Gamma _1+\Gamma _2} = \... ...\mu}\right)}{\infty+\left(\frac{1}{m}+\frac{1}{\mu}\right)} = 0 \end{eqnarray*}

Therefore, the reflection coefficient seen at port 1 of the parallel adaptor is $ \rho = \alpha_1 - 1 = 1$, and the Kelly-Lochbaum scattering junction depicted in Fig.F.20 is verified.

Let's now calculate the internals of the series adaptor in Fig.F.29. From Eq.$ \,$(F.26), the beta parameters are

\begin{eqnarray*} \beta_1 &\isdef & \frac{2R_1}{R_1+R_2+R_3} = \frac{2(m+\mu)}{... ...R_3}{R_1+R_2+R_3} = \frac{2m}{(m+\mu)+m+\mu} = \frac{m}{m+\mu} \end{eqnarray*}

Following Eq.$ \,$(F.30), the series adaptor computes

\begin{eqnarray*} f^{{+}}_J(n) &=& f^{{+}}_1(n)+f^{{+}}_2(n)+f^{{+}}_3(n) = f(... ...(n)\\ &=& \frac{\mu}{m+\mu}f^{{+}}_3(n) - \frac{m}{m+\mu} f(n) \end{eqnarray*}

We do not need to explicitly compute $ f^{{-}}_2(n)$ because it goes into a purely resistive impedance $ \mu $ and produces no return wave. For the same reason, $ f^{{+}}_2(n)\equiv\message{CHANGE eqv TO equiv IN SOURCE}0$. Figure F.30 shows a wave flow diagram of the computations derived, together with the result of elementary simplifications.

Figure F.30: Wave flow diagram for the WDF modeling an ideal force source in parallel with the series combination of a mass $ m$ and dashpot $ \mu $.
\includegraphics{eps/massdashwdf}

Because the difference of the two coefficients in Fig.F.30 is 1, we can easily derive the one-multiply form in Fig.F.31.

Figure: One-multiply form of the WDF in Fig.F.30.
\includegraphics{eps/massdashwdfom}

Checking the WDF against the Analog Equivalent Circuit

Let's check our result by comparing the transfer function from the input force to the force on the mass in both the discrete- and continuous-time cases.

For the discrete-time case, we have

$\displaystyle H_m(z) \isdef \frac{F_3(z)}{F(z)} = \frac{F^{+}_3(z) + F^{-}_3(z)}{F(z)} = (1-z^{-1}) \frac{F^{-}_3(z)}{F(z)} $

where the last simplification comes from the mass reflectance relation $ F^{+}_3(z) = -z^{-1}F^{-}_3(z)$. (Note that we are using the standard traveling-wave notation for the adaptor, so that the $ \pm$ signs are swapped relative to element-centric notation.)

We now need $ F^{-}_3(z)/F(z)$. To simplify notation, define the two coefficients as

\begin{eqnarray*} a &=& \frac{m}{m+\mu}\\ b &=& \frac{\mu}{m+\mu} \end{eqnarray*}

From Figure F.30, we can write

\begin{eqnarray*} F^{-}_3(z) &=& -a\left[F(z)-z^{-1}F^{-}_3(z)\right] + b\left[-... ...\,\,\quad F^{-}_3(z) &=& -a\frac{F(z)}{1-(a-b)z^{-1}F^{-}_3(z)} \end{eqnarray*}

Thus, the desired transfer function is

$\displaystyle H_m(z) = -a \frac{1-z^{-1}}{1-(a-b)z^{-1}} = -\frac{m}{m+\mu} \frac{1-z^{-1}}{1-\left(\frac{m-\mu}{m+\mu}\right)z^{-1}} $

We now wish to compare this result to the bilinear transform of the corresponding analog transfer function. From Figure F.27, we can recognize the mass and dashpot as voltage divider:

$\displaystyle H^a_m(s) = \frac{ms}{ms+\mu} $

Applying the bilinear transform yields

\begin{eqnarray*} H^a_m\left(\frac{1-z^{-1}}{1+z^{-1}}\right) &=& \frac{m\left(\... ...{-1}}{1 - \left(\frac{m-\mu}{m+\mu}\right)z^{-1}}\\ &=& H_m(z) \end{eqnarray*}

Thus, we have verified that the force transfer-function from the driving force to the mass is identical in the discrete- and continuous-time models, except for the bilinear transform frequency warping in the discrete-time case.

Wave Digital Mass-Spring Oscillator

Let's look again at the mass-spring oscillator of §F.3.4, but this time without the driving force (which effectively decouples the mass and spring into separate first-order systems). The physical diagram and equivalent circuit are shown in Fig.F.32 and Fig.F.33, respectively.

Figure F.32: Elementary mass-spring oscillator.
\includegraphics{eps/tank}

Figure F.33: Equivalent circuit for the mass-spring oscillator.
\includegraphics{eps/tankec}

Note that the mass and spring can be regarded as being in parallel or in series. Under the parallel interpretation, we have the WDF shown in Fig.F.34 and Fig.F.35.F.5 The reflection coefficient $ \rho$ can be computed, as usual, from the first alpha parameter:

$\displaystyle \rho = \alpha_1 - 1 = \frac{2\Gamma _1}{\Gamma _1+\Gamma _2} - 1 ... ..._1-\Gamma _2}{\Gamma _1+\Gamma _2} = \frac{R_2-R_1}{R_2+R_1} = \frac{m-k}{m+k} $

This result, $ \rho=(m-k)/(m+k)$, is just the ``impedance step over impedance sum'', so no calculation was really necessary.

Figure F.34: Wave digital mass-spring oscillator.
\includegraphics{eps/tankjunc}

Figure F.35: Detailed wave-flow diagram for the wave digital mass-spring oscillator.
\includegraphics{eps/tankwdf}

Oscillation Frequency

From Fig.F.33, we can see that the impedance of the parallel combination of the mass and spring is given by

$\displaystyle R_{m\vert\vert k}(s) \isdef \left.\frac{k}{s} \right\Vert ms = \frac{\frac{k}{s}ms}{\frac{k}{s}+ms} = \frac{ks}{s^2+\frac{k}{m}} \protect$ (F.38)

(using the product-over-sum rule for combining impedances in parallel). The poles of this impedance are given by the roots of the denominator polynomial in $ s$:

$\displaystyle s = \pm j\sqrt{\frac{k}{m}} \protect$ (F.39)

The resonance frequency of the mass-spring oscillator is therefore

$\displaystyle \omega_0 = \sqrt{\frac{k}{m}} \protect$ (F.40)

Since the poles $ s=\pm j\omega_0$ are on the $ j\omega $ axis, there is no damping, as we expect.

We can now write reflection coefficient $ \rho$ (see Fig.F.35) as

$\displaystyle \rho = \frac{m-k}{m+k} = \frac{1-\frac{k}{m}}{1+\frac{k}{m}} = \frac{1-\omega_0^2}{1+\omega_0^2} $

We see that dc ( $ \omega_0=0$) corresponds to $ \rho=1$, and $ \omega_0=\infty$ corresponds to $ \rho=-1$.

DC Analysis of the WD Mass-Spring Oscillator

Considering the dc case first ($ \rho=1$), we see from Fig.F.35 that the state variable $ x_1(n)$ will circulate unchanged in the isolated loop on the left. Let's call this value $ x_1(n)\equiv\message{CHANGE eqv TO equiv IN SOURCE} x_0$. Then the physical force on the spring is always equal to

$\displaystyle f_k(n) = f^{{+}}_k(n) + f^{{-}}_k(n) = 2x_1(n) = 2 x_0. \qquad\hbox{(spring force, dc case)} \protect$ (F.41)

The loop on the right in Fig.F.35 receives $ 2 x_0$ and adds $ x_2(n)$ to that. Since $ x_2(n+1) = 2x_1(n)+x_2(n)$, we see it is linearly growing in amplitude. For example, if $ x_2(0)=0$ (with $ x_1(0)=x_0$), we obtain $ x_2=[0, 2x_0, 4x_0, 6x_0,\ldots]$, or

$\displaystyle x_2(n) = 2 n x_0, \quad n=0,1,2,3,\ldots\,. \protect$ (F.42)

At first, this result might appear to contradict conservation of energy, since the state amplitude seems to be growing without bound. However, the physical force is fortunately better behaved:

$\displaystyle f_m(n) = f^{{+}}_m(n) + f^{{-}}_m(n) = x_2(n+1) - x_2(n) = 2x_0. \protect$ (F.43)

Since the spring and mass are connected in parallel, it must be the true that they are subjected to the same physical force at all times. Comparing Equations (F.41-F.43) verifies this to be the case.

WD Mass-Spring Oscillator at Half the Sampling Rate

Under the bilinear transform, the $ s=\infty$ maps to $ z=-1$ (half the sampling rate). It is therefore no surprise that given $ \omega_0=\infty$ ($ \rho=-1$), inspection of Fig.F.35 reveals that any alternating sequence (sinusoid sampled at half the sampling rate) will circulate unchanged in the loop on the right, which is now isolated. Let $ x_2(n) = (-1)^n x_0$ denote this alternating sequence. The loop on the left receives $ - 2 x_2(n)$ and adds $ - x_1(n-1)$ to it, i.e., $ x_1(n+1)= - x_1(n) - 2 x_2(n) = -x_1(n) - 2x_2(0)(-1)^n$. If we start out with $ x_1(0)=0$ and $ x_2(0)=x_0$, we obtain $ x_1 = [0,-2x_0, 4x_0, -6x_0, \ldots]$, or

$\displaystyle x_1(n) = (-1)^n 2 n x_0, \quad n=0,1,2,3,\ldots\,. $

However, the physical spring force is well behaved, since

$\displaystyle f_k(n) = f^{{+}}_k(n) + f^{{-}}_k(n) = x_1(n+1) + x_1(n) = (-1)^{n+1}2 x_0 $

As a check, the mass force is found to be

\begin{eqnarray*} f_m(n) &=& x_2(n+1) - x_2(n)\\ &=& (-1)^{n+1}x_0 - (-1)^n x_0\\ &=& (-1)^{n+1}x_0 + (-1)^{n+1}x_0\\ &=& (-1)^{n+1}2 x_0, \end{eqnarray*}

which agrees with the spring, as it must.

Linearly Growing State Variables in WD Mass-Spring Oscillator

It may seem disturbing that such a simple, passive, physically rigorous simulation of a mass-spring oscillator should have to make use of state variables which grow without bound for the limiting cases of simple harmonic motion at frequencies zero and half the sampling rate. This is obviously a valid concern in practice as well. However, it is easy to show that this only happens at dc and $ f_s/2$, and that there is a true degeneracy at these frequencies, even in the physics. For all frequencies in the audio range (e.g., for typical sampling rates), such state variable growth cannot occur. Let's take closer look at this phenomenon, first from a signal processing point of view, and second from a physical point of view.

A Signal Processing Perspective on Repeated Mass-Spring Poles

Going back to the poles of the mass-spring system in Eq.$ \,$(F.39), we see that, as the imaginary part of the two poles, $ \pm\omega_0 = \pm\sqrt{k/m}$, approach zero, they come together at $ s=0$ to create a repeated pole. The same thing happens at $ \omega_0=\infty$ since both poles go to ``the point at infinity''.

It is a well known fact from linear systems theory that two poles at the same point $ s=s_0=\sigma_0$ in the $ s$ plane can correspond to an impulse-response component of the form $ te^{\sigma_0 t}$, in addition to the component $ e^{\sigma_0 t}$ produced by a single pole at $ s=\sigma_0$. In the discrete-time case, a double pole at $ z=r_0$ can give rise to an impulse-response component of the form $ n r_0^n$. This is the fundamental source of the linearly growing internal states of the wave digital sine oscillator at dc and $ f_s/2$. It is interesting to note, however, that such modes are always unobservable at any physical output such as the mass force or spring force that is not actually linearly growing.

Physical Perspective on Repeated Poles in the Mass-Spring System

In the physical system, dc and infinite frequency are in fact strange cases. In the case of dc, for example, a nonzero constant force implies that the mass $ m$ is under constant acceleration. It is therefore the case that its velocity is linearly growing. Our simulation predicts this, since, using Eq.$ \,$(F.43) and Eq.$ \,$(F.42),

\begin{eqnarray*} v_m(n) &=& \frac{f^{{+}}_m(n)}{m} - \frac{f^{{-}}_m(n)}{m} =... ...m} \left[2(n+1) + 2n\right]x_0 = \frac{1}{m} (4 n x_0 + 2 x_0). \end{eqnarray*}

The dc term $ 2x_0/m$ is therefore accompanied by a linearly growing term $ 2nx_0/m$ in the physical mass velocity. It is therefore unavoidable that we have some means of producing an unbounded, linearly growing output variable.

Mass-Spring Boundedness in Reality

To approach the limit of $ \omega_0 = \sqrt{k/m} = 0$, we must either take the spring constant $ k$ to zero, or the mass $ m$ to infinity, or both.

In the case of $ k\to0$, the constant force must approach zero, and we are left with at most a constant mass velocity in the limit (not a linearly growing one, since there can be no dc force at the limit). When the spring force reaches zero, $ x_1(n)=0$, so that only zeros will feed into the loop on the right in Fig.F.35, thus avoiding a linearly growing velocity, as demanded by the physics. (A constant velocity is free to circulate in the loop on the right, but the loop on the left must be zeroed out in the limit.)

In the case of $ m\to\infty$, the mass becomes unaffected by the spring force, so its final velocity must be zero. Otherwise, the attached spring would keep compressing or stretching forever, and this would take infinite energy. (Another way to arrive at this conclusion is to note that the final kinetic energy of the mass would be $ mv^2/2=\infty$.) Since the total energy in an undriven mass-spring oscillator is always constant, the infinite-mass limit must be accompanied by a zero-velocity limit.F.6 This means the mass's state variable $ x_2(n)$ in Fig.F.35 must be forced to zero in the limit so that there will be no linearly growing solution at dc.

In summary, when two or more system poles approach each other to form a repeated pole, care must be taken to ensure that the limit is approached in a physically meaningful way. In the case of the mass-spring oscillator, for example, any change in the spring constant $ k$ or mass $ m$ must be accompanied by the physically appropriate change in the state variables $ x_1(n)$ and/or $ x_2(n)$. It is obviously incorrect, for example, to suddenly set $ k=0$ in the simulation without simultaneously clearing the spring's state variable $ x_1(n)$, since the force across an infinitely compliant spring can only be zero.

Similar remarks apply to repeated poles corresponding to $ \omega_0=\infty$. In this case, the mass and spring basically change places.

Energy-Preserving Parameter Changes (Mass-Spring Oscillator)

If the change in $ k$ or $ m$ is deemed to be ``internal'', that is, involving no external interactions, the appropriate accompanying change in the internal state variables is that which conserves energy. For the mass and its velocity, for example, we must have

$\displaystyle \frac{1}{2} m_1 v_1^2 =\frac{1}{2} m_2 v_2^2 $

where $ m_1,m_2$ denote the mass values before and after the change, respectively, and $ v_1,v_2$ denote the corresponding velocities. The velocity must therefore be scaled according to

$\displaystyle v_2 = v_1\sqrt{\frac{m_1}{m_2}}, $

since this holds the kinetic energy of the mass constant. Note that the momentum of the mass is changed, however, since

$\displaystyle m_2v_2 = m_2 v_1\sqrt{\frac{m_1}{m_2}} =v_1\sqrt{m_1m_2} =m_1v_1\sqrt{\frac{m_2}{m_1}} $

If the spring constant $ k$ is to change from $ k_1$ to $ k_2$, the instantaneous spring displacement $ x$ must satisfy

$\displaystyle \frac{1}{2} k_1 x_1^2 =\frac{1}{2} k_2 x_2^2 $

In a velocity-wave simulation, displacement is the integral of velocity. Therefore, the energy-conserving velocity correction is impulsive in this case.

Exercises in Wave Digital Modeling

  1. Comparing digital and analog frequency formulas. This first exercise verifies that the elementary ``tank circuit'' always resonates at exactly the frequency it should, according to the bilinear transform frequency mapping $ \omega_a = \tan(\omega_d T / 2)$, where $ \omega_a$ denotes ``analog frequency'' and $ \omega_d$ denotes ``digital frequency''.
    1. Find the poles of Fig.F.35 in terms of $ \rho$.

    2. Show that the resonance frequency is given by

      $\displaystyle f_s\arccos\left(\rho\right) $

      where $ f_s$ denotes the sampling rate.

    3. Recall that the mass-spring oscillator resonates at $ \omega_0=\sqrt{k/m}$. Relate these two resonance frequency formulas via the analog-digital frequency map $ \omega_a = \tan(\omega_d T / 2)$.

    4. Show that the trig identity you discovered in this way is true. I.e., show that

      $\displaystyle f_s \arccos\left[\frac{k-m}{k+m}\right] = 2f_s \arctan\left[\sqrt{\frac{m}{k}}\right]. $

Next Section:
Resources on the Internet
Previous Section:
Equivalence of Digital Waveguide and Finite Difference Schemes