The General Linear Time-Invariant Case

Let $y(t,x)$ denote the transverse displacement, in one plane, of a vibrating string [60]. The complete linear time-invariant generalization of the lossy, stiff string is described by the differential equation

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

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

$$\sum_{k=0}^\infty \alpha_k s^k = \sum_{l=0}^\infty \beta_l v^l.$$ (21)

Solving for $v$ in terms of $s$ 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 $s=0$, we normally also have $v=0$ (corresponding to the absence of static deformation in the medium). Stepping $s$ forward by a small differential $j{\scriptstyle\Delta}\omega$, the left-hand side can be approximated by $\alpha_0+\alpha_1j{\scriptstyle\Delta}\omega$. Requiring the generalized wave velocity $s/v(s)$ 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 $s$ steps forward, higher order terms become important one by one on both sides of the equation. Each new term in $v$ spawns a new solution for $v$ in terms of $s$, since the order of the polynomial in $v$ is incremented. For each solution $v(s)$, 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){\tiny\stackrel{\triangle}{=}}\omega/{v_i}(\omega)$, and sampling according to $x\to x_m{\tiny\stackrel{\triangle}{=}}mX$ and $t\to t_n{\tiny\stackrel{\triangle}{=}}nT(\omega)$, with $X{\tiny\stackrel{\triangle}{=}}c(\omega)T(\omega)$, (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

$$e^{j\omega t_n\pm v(j\omega)x_m} = e^{\pm{v_r}(\omega)x_m}\cdot e^{ j\omega\left[t_n\pm x_m/c(\omega)\right]}$$

$$= G^m(\omega)\cdot e^{ j\omega\left(n \pm m\right)T(\omega)}$$ (22)

where $G(\omega){\tiny\stackrel{\triangle}{=}}e^{\pm{v_r}(\omega)X}$. Thus, a general map of $v$ versus $s$, corresponding to a partial differential equation of any order in the form (20), 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 is typical in, say, the stiff string [19].

In summary, 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. Higher order spatial derivatives can be approximated by higher order time derivatives and treated similarly [10]. The corresponding digital simulation of an arbitrarily long (undriven and unobserved) section of 1D medium (such as a string or acoustic tube) can be simplified via commutativity to at most two pure delays and at most two linear, time-invariant filters. In higher dimensions, such as for the 2D mesh, the per-sample filtering cannot in general be exactly commuted to the boundaries of the mesh. However, an approximation problem can be solved which matches the observed modal frequencies and decay rates using sparsely distributed low-order filters in an otherwise lossless mesh (e.g., around the rim).

In practical physical simulation scenarios, such as for real-world strings or acoustic tubes, it is generally most effective to identify experimentally the attenuation and dispersion associated with wave propagation at each frequency over the band of interest [108,49]. These data can be used to design a digital filter which gives an optimal approximation over the propagation distance desired. If needed, the per-sample filter identified in this way can be translated into higher-order terms in the wave equation for the medium. Thus, the wave equation itself can be ``identified'' from measured input-output behavior of the medium (assuming it is linear and uniform) rather than being derived from physical principles and physical constants of the medium as is classically done [61,11].