Free Books
   Free Books    Physical Audio Signal Processing  

Digitization of Lumped Models

Since lumped models are described by differential equations, they are digitized (brought into the digital-signal domain) by converting them to corresponding finite-difference equations (or simply ``difference equations''). General aspects of finite difference schemes are discussed in Appendix D. This chapter introduces a couple of elementary methods in common use:

  1. Finite difference approximation: $ sT\leftarrow 1-z^{-1}$
  2. Bilinear transform: $ sT\leftarrow 2\,(1-z^{-1})/(1+z^{-1})$
Note that digitization by the bilinear transform is closely related to the Wave Digital Filter (WDF) method introduced in Appendix F. Section 9.3.1 discusses a bilinearly transformed mass colliding with a digital waveguide string (an idealized struck-string example).

Finite Difference Approximation

A finite difference approximation (FDA) approximates derivatives with finite differences, i.e.,

$\displaystyle \frac{d}{dt} x(t) \isdefs \lim_{\delta\to 0} \frac{x(t) - x(t-\delta)}{\delta} \;\approx\; \frac{x(n T)-x[(n-1)T]}{T} \protect$ (8.2)

for sufficiently small $ T$.8.5

Equation (7.2) is also known as the backward difference approximation of differentiation.

See §C.2.1 for a discussion of using the FDA to model ideal vibrating strings.

FDA in the Frequency Domain

Viewing Eq.$ \,$(7.2) in the frequency domain, the ideal differentiator transfer-function is $ H(s)=s$, which can be viewed as the Laplace transform of the operator $ d/dt$ (left-hand side of Eq.$ \,$(7.2)). Moving to the right-hand side, the z transform of the first-order difference operator is $ (1-z^{-1})/T$. Thus, in the frequency domain, the finite-difference approximation may be performed by making the substitution

$\displaystyle s \;\leftarrow\; \frac{1-z^{-1}}{T} \protect$ (8.3)

in any continuous-time transfer function (Laplace transform of an integro-differential operator) to obtain a discrete-time transfer function (z transform of a finite-difference operator).

The inverse of substitution Eq.$ \,$(7.3) is

$\displaystyle z \eqsp \frac{1}{1 - sT} \eqsp 1 + sT+ (sT)^2 + \cdots \, . $

As discussed in §8.3.1, the FDA is a special case of the matched $ z$ transformation applied to the point $ s=0$.

Note that the FDA does not alias, since the conformal mapping $ s = {1 - z^{-1}}$ is one to one. However, it does warp the poles and zeros in a way which may not be desirable, as discussed further on p. [*] below.

Delay Operator Notation

It is convenient to think of the FDA in terms of time-domain difference operators using a delay operator notation. The delay operator $ d$ is defined by

$\displaystyle d^k x(n) \eqsp x(n-k). $

Thus, the first-order difference (derivative approximation) is represented in the time domain by $ (1-d)/T$. We can think of $ d$ as $ z^{-1}$ since, by the shift theorem for $ z$ transforms, $ z^{-k} X(z)$ is the $ z$ transform of $ x$ delayed (right shifted) by $ k$ samples.

The obvious definition for the second derivative is

$\displaystyle {\hat{\ddot x}}(n) \eqsp \frac{(1-d)^2}{T^2} x(n).$ (8.4)

However, a better definition is the centered finite difference

$\displaystyle {\hat{\ddot x}}(n) \isdefs \frac{(d^{-1}-1)(1-d)}{T^2} x(n) \eqsp \frac{d^{-1}-2+d}{T^2}x(n) \protect$ (8.5)

where $ d^{-1}$ denotes a unit-sample advance. This definition is preferable as long as one sample of look-ahead is available, since it avoids an operator delay of one sample. Equation (7.5) is a zero phase filter, meaning it has no delay at any frequency, while (7.4) is a linear phase filter having a delay of $ 1$ sample at all frequencies.

Bilinear Transformation

The bilinear transform is defined by the substitution

$\displaystyle s$ $\displaystyle =$ $\displaystyle c\frac{1-z^{-1}}{1+z^{-1}}, \quad c>0, \; c=2/T\;$   (typically) (8.6)
$\displaystyle \,\,\Rightarrow\,\, z$ $\displaystyle =$ $\displaystyle \frac{1+s/c}{1-s/c} \eqsp 1 + 2(s/c) + 2(s/c)^2 + 2(s/c)^3 + \cdots$ (8.7)

where $ c$ is some positive constant [83,326]. That is, given a continuous-time transfer function $ H_a(s)$, we apply the bilinear transform by defining

$\displaystyle H_d(z) = H_a\left(c\frac{1-z^{-1}}{1+z^{-1}}\right) $

where the ``$ d$'' subscript denotes ``digital,'' and ``$ a$'' denotes ``analog.''

It can be seen that analog dc ($ s=0$) maps to digital dc ($ z=1$) and the highest analog frequency ($ s=\infty$) maps to the highest digital frequency ($ z=-1$). It is easy to show that the entire $ j\omega $ axis in the $ s$ plane (where $ s\isdeftext \sigma+j\omega$) is mapped exactly once around the unit circle in the $ z$ plane (rather than summing around it infinitely many times, or ``aliasing'' as it does in ordinary sampling). With $ c$ real and positive, the left-half $ s$ plane maps to the interior of the unit circle, and the right-half $ s$ plane maps outside the unit circle. This means stability is preserved when mapping a continuous-time transfer function to discrete time.

Another valuable property of the bilinear transform is that order is preserved. That is, an $ N$th-order $ s$-plane transfer function carries over to an $ N$th-order $ z$-plane transfer function. (Order in both cases equals the maximum of the degrees of the numerator and denominator polynomials [449]).8.6

The constant $ c$ provides one remaining degree of freedom which can be used to map any particular finite frequency from the $ j\omega_a$ axis in the $ s$ plane to a particular desired location on the unit circle $ e^{j\omega_d}$ in the $ z$ plane. All other frequencies will be warped. In particular, approaching half the sampling rate, the frequency axis compresses more and more. Note that at most one resonant frequency can be preserved under the bilinear transformation of a mass-spring-dashpot system. On the other hand, filters having a single transition frequency, such as lowpass or highpass filters, map beautifully under the bilinear transform; one simply uses $ c$ to map the cut-off frequency where it belongs, and the response looks great. In particular, ``equal ripple'' is preserved for optimal filters of the elliptic and Chebyshev types because the values taken on by the frequency response are identical in both cases; only the frequency axis is warped.

The bilinear transform is often used to design digital filters from analog prototype filters [343]. An on-line introduction is given in [449].

Finite Difference Approximation vs. Bilinear Transform

Recall that the Finite Difference Approximation (FDA) defines the elementary differentiator by $ y(n) = x(n) - x(n-1)$ (ignoring the scale factor $ T$ for now), and this approximates the ideal transfer function $ H(s)=s$ by $ H_d(z)=1-z^{-1}$. The bilinear transform calls instead for the transfer function $ H'_d(z)=(1-z^{-1})/(1+z^{-1})$ (again dropping the scale factor) which introduces a pole at $ z=-1$ and gives us the recursion $ y(n) = x(n) - x(n-1) - y(n-1)$. Note that this new pole is right on the unit circle and is therefore undamped. Any signal energy at half the sampling rate will circulate forever in the recursion, and due to round-off error, it will tend to grow. This is therefore a potentially problematic revision of the differentiator. To get something more practical, we need to specify that the filter frequency response approximate $ H(j\omega)=j\omega$ over a finite range of frequencies $ [-\omega_c,\omega_c]$, where $ \omega_c\ll\pi f_s$, above which we allow the response to ``roll off'' to zero. This is how we pose the differentiator problem in terms of general purpose filter design (see §8.6) [362].

To understand the properties of the finite difference approximation in the frequency domain, we may look at the properties of its $ s$-plane to $ z$-plane mapping

$\displaystyle s = \frac{1 - z^{-1}}{T}. $

We see the FDA is actually a portion of the bilinear transform, since following the FDA mapping by the mapping $ s = (c/T)/(1+z^{-1})$ would convert it to the bilinear transform. Like the bilinear transform, the FDA does not alias, since the mapping $ s = 1 - z^{-1}$ is one-to-one.

Setting $ T$ to 1 for simplicity and solving the FDA mapping for z gives

$\displaystyle z = \frac{1 }{1-s}. $

We see that analog dc ($ s=0$) maps to digital dc ($ z=1$) as desired, but higher frequencies unfortunately map inside the unit circle rather than onto the unit circle in the $ z$ plane. Solving for the image in the z plane of the $ j\omega $ axis in the s plane gives

$\displaystyle z = \frac{1}{1-j \omega } = \frac{1 + j \omega }{1+\omega^2} $

From this it can be checked that the FDA maps the $ j\omega $ axis in the $ s$ plane to the circle of radius $ 1/2$ centered at the point $ z = 1/2$ in the $ z$ plane, as shown in Fig. 7.15

Figure 7.15: Image of the $ j\omega $ axis in the $ z$ plane: a circle of radius $ 1/2$ centered at the point $ z = 1/2$.
\includegraphics[width=3in]{eps/lfdacirc}

Under the FDA, analog and digital frequency axes coincide well enough at very low frequencies (high sampling rates), but at high frequencies relative to the sampling rate, artificial damping is introduced as the image of the $ j\omega $ axis diverges away from the unit circle.

While the bilinear transform ``warps'' the frequency axis, we can say the FDA ``doubly warps'' the frequency axis: It has a progressive, compressive warping in the direction of increasing frequency, like the bilinear transform, but unlike the bilinear transform, it also warps normal to the frequency axis.

Consider a point traversing the upper half of the unit circle in the z plane, starting at $ z=1$ and ending at $ z=-1$. At dc, the FDA is perfect, but as we proceed out along the unit circle, we diverge from the $ j\omega $ axis image and carve an arc somewhere out in the image of the right-half $ s$ plane. This has the effect of introducing an artificial damping.

Consider, for example, an undamped mass-spring system. There will be a complex conjugate pair of poles on the $ j\omega $ axis in the $ s$ plane. After the FDA, those poles will be inside the unit circle, and therefore damped in the digital counterpart. The higher the resonance frequency, the larger the damping. It is even possible for unstable $ s$-plane poles to be mapped to stable $ z$-plane poles.

In summary, both the bilinear transform and the FDA preserve order, stability, and positive realness. They are both free of aliasing, high frequencies are compressively warped, and both become ideal at dc, or as $ f_s$ approaches $ \infty$. However, at frequencies significantly above zero relative to the sampling rate, only the FDA introduces artificial damping. The bilinear transform maps the continuous-time frequency axis in the $ s$ (the $ j\omega $ axis) plane precisely to the discrete-time frequency axis in the $ z$ plane (the unit circle).

Application of the Bilinear Transform

The impedance of a mass in the frequency domain is

$\displaystyle R_M(s) = Ms. $

In the $ s$ plane, we have

$\displaystyle F_a(s) = (Ms) V_a(s) $

where the ``a'' subscript denotes ``analog''. For simplicity, let's choose the free constant $ c$ in the bilinear transform such that $ 1$ rad/sec maps to one fourth the sampling rate, i.e., $ s=j$ maps to $ z=j$ which implies $ c=1$. Then the impedance relation maps across as

$\displaystyle F_d(z) = \left(M\frac{1-z^{-1}}{1+z^{-1}}\right) V_d(z) $

where the ``d'' subscript denotes ``digital. Multiplying through by the denominator and applying the shift theorem for $ z$ transforms gives the corresponding difference equation

\begin{eqnarray*} (1+z^{-1})F_d(z) &=& M (1-z^{-1}) V_d(z) \\ \;\longleftrighta... ... \,\,\Rightarrow\,\,f_d(n) &=& M[v_d(n) - v_d(n-1)] - f_d(n-1). \end{eqnarray*}

This difference equation is diagrammed in Fig. 7.16. We recognize this recursive digital filter as the direct form I structure. The direct-form II structure is obtained by commuting the feedforward and feedback portions and noting that the two delay elements contain the same value and can therefore be shared [449]. The two other major filter-section forms are obtained by transposing the two direct forms by exchanging the input and output, and reversing all arrows. (This is a special case of Mason's Gain Formula which works for the single-input, single-output case.) When a filter structure is transposed, its summers become branching nodes and vice versa. Further discussion of the four basic filter section forms can be found in [333].

Figure 7.16: A direct-form-I digital filter simulating a mass $ M$ created using the bilinear transform $ s=(1-z^{-1})/(1+z^{-1})$.
\includegraphics[width=4in]{eps/lmassFilterDF1}

Practical Considerations

While the digital mass simulator has the desirable properties of the bilinear transform, it is also not perfect from a practical point of view: (1) There is a pole at half the sampling rate ($ z=-1$). (2) There is a delay-free path from input to output.

The first problem can easily be circumvented by introducing a loss factor $ g$, moving the pole from $ z=-1$ to $ z=-g$, where $ g\in[0,1)$ and $ g\approx1$. This is sometimes called the ``leaky integrator''.

The second problem arises when making loops of elements (e.g., a mass-spring chain which forms a loop). Since the individual elements have no delay from input to output, a loop of elements is not computable using standard signal processing methods. The solution proposed by Alfred Fettweis is known as ``wave digital filters,'' a topic taken up in §F.1.

Limitations of Lumped Element Digitization

Model discretization by the FDA (§7.3.1) and bilinear transform (§7.3.2) methods are order preserving. As a result, they suffer from significant approximation error, especially at high frequencies relative to half the sampling rate. By allowing a larger order in the digital model, we may obtain arbitrarily accurate transfer-function models of LTI subsystems, as discussed in Chapter 8. Of course, in higher-order approximations, the state variables of the simulation no longer have a direct physical intepretation, and this can have implications, particularly when trying to extend to the nonlinear case. The benefits of a physical interpretation should not be given up lightly. For example, one may consider oversampling in place of going to higher-order element approximations.

Next Section:
More General Finite-Difference Methods
Previous Section:
One-Port Network Theory