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 , we obtain the
modified wave equation
Thus, the wave equation has been extended by a ``first-order'' term, i.e., a term proportional to the first derivative of
![$ y$](http://www.dsprelated.com/josimages_new/pasp/img184.png)
![$ {\dddot y}$](http://www.dsprelated.com/josimages_new/pasp/img3323.png)
![$ {\partial^5 y/\partial t^5}$](http://www.dsprelated.com/josimages_new/pasp/img3324.png)
Setting
in the wave equation to find the relationship
between temporal and spatial frequencies in the eigensolution, the wave
equation becomes
![]() |
![]() |
![]() |
|
![]() |
![]() |
![]() |
|
![]() |
![]() |
![]() |
where
![$ c \isdef \sqrt{K/\epsilon }$](http://www.dsprelated.com/josimages_new/pasp/img3331.png)
![$ \vert s\vert$](http://www.dsprelated.com/josimages_new/pasp/img3332.png)
![$ \mu $](http://www.dsprelated.com/josimages_new/pasp/img60.png)
![$ \epsilon $](http://www.dsprelated.com/josimages_new/pasp/img186.png)
![]() |
(C.22) |
by the binomial theorem. For this small-loss approximation, we obtain the following relationship between temporal and spatial frequency:
![]() |
(C.23) |
The eigensolution is then
![]() |
(C.24) |
The right-going part of the eigensolution is
![]() |
(C.25) |
which decays exponentially in the direction of propagation, and the left-going solution is
![]() |
(C.26) |
which also decays exponentially in the direction of travel.
Setting and using superposition to build up arbitrary traveling
wave shapes, we obtain the general class of solutions
![]() |
(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*}](http://www.dsprelated.com/josimages_new/pasp/img3339.png)
The simulation diagram for the lossy digital waveguide is shown in Fig.C.5.
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
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
), 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.
![]() |
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]:
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
in which the third-order partial derivative with respect to time,
![$ {\dddot y}$](http://www.dsprelated.com/josimages_new/pasp/img3323.png)
![$ x$](http://www.dsprelated.com/josimages_new/pasp/img179.png)
![$ t$](http://www.dsprelated.com/josimages_new/pasp/img122.png)
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
. Because
propagation is passive, we will always have
.
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,
$](http://www.dsprelated.com/josimages_new/pasp/img3348.png)
![$ n>1$](http://www.dsprelated.com/josimages_new/pasp/img3349.png)
![$\displaystyle G(\omega) \propto \frac{1}{\omega^{n-1}},$](http://www.dsprelated.com/josimages_new/pasp/img3350.png)
![$\displaystyle \mbox{($n$\ odd)}$](http://www.dsprelated.com/josimages_new/pasp/img3351.png)
![$\displaystyle .
$](http://www.dsprelated.com/josimages_new/pasp/img162.png)
![$ \mu_3{\dddot y}$](http://www.dsprelated.com/josimages_new/pasp/img3352.png)
![$ \mu_5{\partial^5 y/\partial t^5}$](http://www.dsprelated.com/josimages_new/pasp/img3353.png)
![$ G(\omega)$](http://www.dsprelated.com/josimages_new/pasp/img3346.png)
![$\displaystyle G(\omega) = g_0 + g_2 \omega^2 + g_4 \omega^4
$](http://www.dsprelated.com/josimages_new/pasp/img3354.png)
![$ g_i$](http://www.dsprelated.com/josimages_new/pasp/img717.png)
![$ \mu_3$](http://www.dsprelated.com/josimages_new/pasp/img3355.png)
![$ \mu_5$](http://www.dsprelated.com/josimages_new/pasp/img3356.png)
![$ g$](http://www.dsprelated.com/josimages_new/pasp/img37.png)
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]
Thus, to the ideal string wave equation Eq.
![$ \,$](http://www.dsprelated.com/josimages_new/pasp/img196.png)
![$ x$](http://www.dsprelated.com/josimages_new/pasp/img179.png)
![$ x$](http://www.dsprelated.com/josimages_new/pasp/img179.png)
![$ t$](http://www.dsprelated.com/josimages_new/pasp/img122.png)
![$ t$](http://www.dsprelated.com/josimages_new/pasp/img122.png)
Digital Filter Models of Damped Strings
In an efficient digital simulation, lumped loss factors of the form
are approximated by a rational frequency response
. In general, the coefficients of the optimal rational
loss filter are obtained by minimizing
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
to
be finite in length (i.e., an FIR filter) and it must be symmetric
about time zero, i.e.,
. 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].
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*}](http://www.dsprelated.com/josimages_new/pasp/img3364.png)
Adding these equations gives
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},
$](http://www.dsprelated.com/josimages_new/pasp/img3368.png)
![$ g^{+}_m= g^{-}_m= g$](http://www.dsprelated.com/josimages_new/pasp/img3369.png)
![$ a_m= -g^2$](http://www.dsprelated.com/josimages_new/pasp/img3370.png)
![$ \vert g\vert\leq 1$](http://www.dsprelated.com/josimages_new/pasp/img3371.png)
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.
$](http://www.dsprelated.com/josimages_new/pasp/img3372.png)
![$ g(n)$](http://www.dsprelated.com/josimages_new/pasp/img451.png)
![$ G(z)$](http://www.dsprelated.com/josimages_new/pasp/img13.png)
![$ \,$](http://www.dsprelated.com/josimages_new/pasp/img196.png)
![\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*}](http://www.dsprelated.com/josimages_new/pasp/img3373.png)
where 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*}](http://www.dsprelated.com/josimages_new/pasp/img3374.png)
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}.
$](http://www.dsprelated.com/josimages_new/pasp/img3375.png)
![$ y_{n,m}$](http://www.dsprelated.com/josimages_new/pasp/img3376.png)
![$ G(z)$](http://www.dsprelated.com/josimages_new/pasp/img13.png)
![$ n$](http://www.dsprelated.com/josimages_new/pasp/img146.png)
![$ y^f_{n,m}$](http://www.dsprelated.com/josimages_new/pasp/img3377.png)
![$ y^{ff}_{n-1,m}$](http://www.dsprelated.com/josimages_new/pasp/img3378.png)
![$ m$](http://www.dsprelated.com/josimages_new/pasp/img6.png)
![$ y_{n+1,m}$](http://www.dsprelated.com/josimages_new/pasp/img3379.png)
![$ y^f_{n,m}$](http://www.dsprelated.com/josimages_new/pasp/img3377.png)
![$ \forall m$](http://www.dsprelated.com/josimages_new/pasp/img3380.png)
![$ G(z)$](http://www.dsprelated.com/josimages_new/pasp/img13.png)
![$ y^{ff}_{n-1,m}$](http://www.dsprelated.com/josimages_new/pasp/img3378.png)
![$ y_{n+1,m}$](http://www.dsprelated.com/josimages_new/pasp/img3379.png)
![$ G(z)$](http://www.dsprelated.com/josimages_new/pasp/img13.png)
![$ y^f_{n,m}$](http://www.dsprelated.com/josimages_new/pasp/img3377.png)
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.
Next Section:
The Dispersive 1D Wave Equation
Previous Section:
Sampled Traveling Waves