# The simulation of piano string vibration: From physical models to finite difference schemes and digital waveguides

Julien Bensa, Stefan Bilbao, Richard Kronland-Martinet, Julius O. Smith Iii

<figure><img src="image_1.png" /></figure>

## To cite this version:

Julien Bensa, Stefan Bilbao, Richard Kronland-Martinet, Julius O. Smith Iii. The simulation of piano string vibration: From physical models to finite difference schemes and digital waveguides. Journal of the Acoustical Society of America, 2003, 114 (2), pp.1095-1107. 10.1121/1.1587146. hal-00088329

HAL Id: hal-00088329

https://hal.science/hal-00088329v1

Submitted on 1 Aug 2006

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# The simulation of piano string vibration: From physical models to finite difference schemes and digital waveguides

Julien Bensa, Stefan Bilbao, Richard Kronland-Martinet, and Julius O. Smith III

Center for Computer Research in Music and Acoustics (CCRMA), Department of Music, Stanford University, Stanford, California 94305-8180

A model of transverse piano string vibration, second order in time, which models frequency-dependent loss and dispersion effects is presented here. This model has many desirable properties, in particular that it can be written as a well-posed initial-boundary value problem (permitting stable finite difference schemes) and that it may be directly related to a digital waveguide model, a digital filter-based algorithm which can be used for musical sound synthesis. Techniques for the extraction of model parameters from experimental data over the full range of the grand piano are discussed, as is the link between the model parameters and the filter responses in a digital waveguide. Simulations are performed. Finally, the waveguide model is extended to the case of several coupled strings.

## I. INTRODUCTION

Several models of transverse wave propagation on a stiff string, of varying degrees of complexity, have appeared in the literature.¹⁻⁵ These models, intended for the synthesis of musical tones, are always framed in terms of a partial differential equation (PDE), or system of PDEs; usually, the simplified starting point for such a model is the *one-dimensional wave equation*,⁶ and the more realistic features, such as dispersion, frequency-dependent loss and nonlinear hammer excitation (in the case of the piano string), are incorporated through several extra terms. Chaigne and Askenfelt³ have proposed such a model (see the Appendix for a concise description of this model), and used it as the basis for a synthesis technique, through the use of finite differences—the time waveform on a struck piano string is simulated in this way to a remarkable degree of fidelity.⁴ Frequency-dependent loss, the feature of primary interest in this paper, is modeled through the use of a third-order time derivative perturbation to the dispersive wave equation; while a physical justification for the use of such a term is tenuous, it does give rise to perceptually important variations in damping rates.

In Sec. II, we introduce a model of string vibration, which is of second order in time differentiation; frequency-dependent loss is introduced via mixed time-space derivative terms. As it turns out, the model discussed here is a substantial improvement in several different ways. First, the frequency domain analysis of a second-order system is quite straightforward, and it is quite easy to obtain explicit formulas for dispersion and loss curves; this is considerably more complicated for systems which are higher order in time, essentially requiring the factorization of a higher order polynomial dispersion relation. Second, it is easy to prove that our model, when complemented by initial and boundary conditions, is *well posed*.⁷,⁸ Though we do not give a complete description of this condition here, to say that such an initial-boundary value problem is well posed is to say, generally speaking, that solutions may not grow faster than exponen-

tially; for linear and shift-invariant systems such as simple stiff string models, the condition can be conveniently expressed in the frequency domain. We show in the Appendix that the PDE first proposed by Ruiz,⁹ and later popularized by Chaigne and Askenfelt³,⁴ is in fact not well-posed, and possesses a spurious unstable solution. Third, it becomes easy to develop *finite difference schemes*, for which precise numerical stability conditions may be easily derived. Finite difference schemes are discussed briefly in Sec. III. Fourth, it is possible to extend the model described here to a more realistic representation of dispersion and loss as a function of frequency through additional terms in the PDE, without compromising well-posedness, or the stability of resulting numerical schemes.

Finally, it is possible to identify the model PDE with a *digital waveguide*—this filterlike structure models one-dimensional wave propagation as purely lossless throughout the length of the string, with loss and dispersion lumped in terminating filters. It thus performs a simulation of a *modified* physical system.

We show, in Sec. IV A, how one can relate the PDE model presented in Sec. II to a digital waveguide structure, paying particular attention to the relationship between the lumped filters used to model loss and dispersion and the model parameters which define our PDE. In Sec. V, we perform several numerical simulations in order to compare the finite difference and waveguide approaches for this particular problem. In particular, in Sec. VI, we examine in detail a procedure allowing the resynthesis of natural string vibration. Using experimental data obtained from a grand piano, both the terminating filters of a digital waveguide and the parameters of the physical model are estimated over most of the keyboard range. This leads to a simple description of the variation of some of these parameters (namely loss parameters and string stiffness) over the piano's range which can be used for the convenient synthesis string vibrations at a given excitation point. Finally, in Sec. VII, interstring coupling is discussed and modeled using coupled digital
waveguides. This is a further step towards the design of a
realistic piano simulator, which should ultimately also model
effects such as nonlinear hammer-string coupling and sound-
board radiation phenomena.

## II. SECOND-ORDER MODELS OF ONE-DIMENSIONAL WAVE PROPAGATION

### A. A family of PDEs

Consider a general linear second-order (in time) wave
equation, of the form

$$ \frac{\partial^2 y}{\partial t^2} + 2 \sum_{k=0}^{M} q_k \frac{\partial^{2k+1} y}{\partial x^{2k} \partial t} + \sum_{k=0}^{N} r_k \frac{\partial^{2k} y}{\partial x^{2k}} = 0. \quad (1) $$

Here, $y(x,t)$, the solution, is a function of position $x \in \mathbb{R}$ and time $t \ge 0$, and $q_k$ and $r_k$ are real constants; the solution is not uniquely defined until two initial conditions are supplied. (For the moment, we concentrate on the pure initial value problem and ignore boundary conditions—we will return to this subject in Sec. II B.) Because this equation describes a linear and shift-invariant system, it is possible to analyze it through Fourier techniques. Defining the spatial Fourier transform $\hat{y}(\beta, t)$ of $y(x,t)$ by

$$ \hat{y}(\beta, t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} y(x,t) e^{-j\beta x} dx, $$

Eq. (1) can be rewritten as

$$ \frac{\partial^2 \hat{y}}{\partial t^2} + 2 \sum_{k=0}^{M} (j\beta)^{2k} q_k \frac{\partial \hat{y}}{\partial t} + \sum_{k=0}^{N} (j\beta)^{2k} r_k \hat{y} = 0. \quad (2) $$

This second-order linear ordinary differential equation with
real coefficients will have solutions of the form

$$ \hat{y}(\beta, t) = \hat{y}_0(\beta) e^{st} \quad (3) $$

for complex frequencies $s$ which satisfy the characteristic
polynomial equation

$$ s^2 + 2q(\beta)s + r(\beta) = 0 \quad (4) $$

with

$$ q(\beta) = \sum_{k=0}^{M} (j\beta)^{2k} q_k, \quad r(\beta) = \sum_{k=0}^{N} (j\beta)^{2k} r_k. $$

Notice that because only even derivatives appear in the fam-
ily Eq. (1), the functions $q$ and $r$ are real. The characteristic
polynomial equation has roots

$$ s_{\pm} = -q \pm \sqrt{q^2 - r}. $$

The condition that the initial value problem corresponding to Eq. (1) be *well-posed* is that these roots have real parts which are bounded from above as a function of $\beta$; this is in effect saying that solution growth can be no faster than exponential. A more restrictive (and physically relevant) condition is that these roots have nonpositive real part for all $\beta$, so that all exponential solutions are nonincreasing. It is simple to show that this will be true for

$$ q(\beta), r(\beta) \ge 0. \quad (5) $$

For $q$ and $r$ satisfying Eq. (5), the imaginary parts of

<table><caption>TABLE I. Physical model parameters for piano tones C2, C4, and C7.</caption><thead><tr><th></th><th>C2</th><th>C4</th><th>C7</th><th>Units</th></tr></thead><tbody><tr><td>L</td><td>1.23</td><td>0.63</td><td>0.10</td><td>m</td></tr><tr><td>c</td><td>160.9</td><td>329.6</td><td>418.6</td><td>m s<sup>-1</sup></td></tr><tr><td>κ</td><td>0.58</td><td>1.25</td><td>1.24</td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>b<sub>1</sub></td><td>0.25</td><td>1.1</td><td>9.17</td><td>s<sup>-1</sup></td></tr><tr><td>b<sub>2</sub></td><td>7.5 × 10<sup>-5</sup></td><td>2.7 × 10<sup>-4</sup></td><td>2.1 × 10<sup>-3</sup></td><td>m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>F<sub>s</sub></td><td>16 000</td><td>32 000</td><td>96 000</td><td>s<sup>-1</sup></td></tr></tbody></table>

these roots correspond to oscillation frequencies, and the real
parts to loss. Clearly, for real wave numbers $\beta$ such that $q^2$
$\le r$, the real parts of $s_\pm$ are simply $-q$. This case corre-
sponds to normal damped wave propagation; notice in par-
ticular that if $q$ depends on $\beta$, then damping rates will be
wave number (and thus frequency) dependent. If $q^2 > r$, then
both roots are purely real and nonpositive, yielding damped
nontraveling solutions.

Consider a member of the family defined by Eq. (1),

$$ \frac{\partial^2 y}{\partial t^2} = c^2 \frac{\partial^2 y}{\partial x^2} - \kappa^2 \frac{\partial^4 y}{\partial x^4} - 2b_1 \frac{\partial y}{\partial t} + 2b_2 \frac{\partial^3 y}{\partial x^2 \partial t}. \quad (6) $$

The first term on the right-hand side of the equation, in the
absence of the others, gives rise to wavelike motion, with
speed $c$. The second “ideal bar” term¹⁰ introduces disper-
sion, or frequency-dependent wave velocity, and is param-
etrized by a stiffness coefficient $\kappa$. The third and fourth terms
allow for loss, and if $b_2 \ne 0$, decay rates will be frequency
dependent. [A complete model, for a piano string, is obtained
by including a hammer excitation term, $f(x,t)$, possibly ac-
counting for nonlinear effects, on the right-hand side, and by
restricting the spatial domain to a finite interval and supply-
ing a realistic set of boundary conditions.] This model differs
from that of Ruiz⁹ in only the last term (see the Appendix).

The characteristic equation has the form of Eq. (4), with

$$ q(\beta) = b_1 + b_2 \beta^2, \quad r(\beta) = c^2 \beta^2 + \kappa^2 \beta^4. $$

For $b_1, b_2 \ge 0$, condition Eq. (5) is satisfied and this PDE
obviously possesses exponentially decaying solutions, and
what is more, loss increases as a function of wave number.
The PDE of Eq. (6) possesses traveling wave solutions when
$q^2 \le r$, which, for realistic values of the defining parameters
for a piano string, includes the overwhelming part of the
audio spectrum. For instance, for a C2 piano string, de-
scribed by the parameters given in Table I, the lower cutoff
wave number for traveling waves will be $\beta=0.0025$, corre-
sponding to a frequency of 0.080 Hz. There is no upper
cutoff.

In order to relate this PDE model with a digital wave-
guide numerical simulation method, it is useful to write the
expressions for dispersion and loss directly. Taking

$$ s_{\pm} = \sigma \pm j\omega \quad (7) $$

over the range of $\beta$ for which traveling wave solutions exist,
we obtain

$$ \sigma(\beta) = -b_1 - b_2 \beta^2, \quad (8) $$

$$ \omega(\beta) = \sqrt{-(b_1 + b_2 \beta^2)^2 + c^2 \beta^2 + \kappa^2 \beta^4}. \quad (9) $$

We will discuss digital waveguide models in detail in Sec.
IV A.
Finally, we mention that the general form Eq. (1) serves
as a useful point of departure for more accurate models of
loss and dispersion in a stiff string. The more terms are in-
cluded, the better these phenomena can be modeled over the
entire frequency range of interest. Though it is difficult to
associate physical processes directly with the various extra
perturbation terms in the equation, it is at least simple to
ensure, through condition Eq. (5), that the model is well-
posed, an important first step in developing stable numerical
methods.

B. Boundary conditions

In this section, we provide a brief analysis of pinned
boundary conditions, to show that when coupled to our string
model, the same analysis of well-posedness may be applied
(i.e., wave numbers of modal solutions are real). Let us now
restrict the spatial domain for the problem defined by Eq. (6)
to $x \in [0, L]$. As Eq. (6) is of fourth order in the spatial de-
rivatives, we need to supply two boundary conditions at ei-
ther end, i.e., at $x=0$ and $x=L$. Following Chaigne,³ we
apply pinned boundary conditions,

$$
y|_{x=0} = y|_{x=L} = \left. \frac{\partial^2 y}{\partial x^2} \right|_{x=0} = \left. \frac{\partial^2 y}{\partial x^2} \right|_{x=L} = 0. \quad (10)
$$

For a solution of the form $y(x,t)=e^{st+j\beta x}$, from dispersion relation Eq. (4), there are thus four solutions for $\beta$ in terms of $s$,

$$
\beta_{+}^{2}(s) = \frac{-\gamma + \sqrt{\gamma^{2} - 4\kappa^{2}(s^{2} + 2b_{1}s)}}{2\kappa^{2}}, \quad (11a)
$$

$$
\beta_{-}^{2}(s) = \frac{-\gamma - \sqrt{\gamma^{2} - 4\kappa^{2}(s^{2} + 2b_{1}s)}}{2\kappa^{2}} \quad (11b)
$$

with

$$
\gamma = c^2 + 2b_2 s.
$$

At frequency $s$, thus, any linear combination

$$
y(x,t) = e^{st}(a_{+,+}e^{j\beta_{+}x} + a_{+,−}e^{-j\beta_{+}x} + a_{−,+}e^{j\beta_{−}x} \\
\qquad + a_{−,−}e^{-j\beta_{−}x}) \quad (12)
$$

is a solution to Eq. (6). Applying the boundary conditions
Eq. (10) to this solution gives the matrix equation

$$
\underbrace{\begin{bmatrix} 1 & 1 & 1 & 1 \\ e^{j\beta_{+}L} & e^{-j\beta_{+}L} & e^{j\beta_{-}L} & e^{-j\beta_{-}L} \\ -\beta_{+}^2 & -\beta_{+}^2 & -\beta_{-}^2 & -\beta_{-}^2 \\ -\beta_{+}^2 e^{j\beta_{+}L} & -\beta_{+}^2 e^{-j\beta_{+}L} & -\beta_{-}^2 e^{j\beta_{-}L} & -\beta_{-}^2 e^{-j\beta_{-}L} \end{bmatrix}}_{\triangleq A(s)} \times \begin{bmatrix} a_{+,+} \\ a_{+,−} \\ a_{−,+} \\ a_{−,−} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}. \quad (13)
$$

Nontrivial solutions can occur only when det(A)=0, giving
the relation

$$
(\beta_{+}^{2} - \beta_{-}^{2})^{2} \sin(\beta_{+}L) \sin(\beta_{-}L) = 0. \quad (14)
$$

Discounting the case $\beta_+^2 = \beta_-^2$ [which yields an identically zero solution $y(x,t)$], then obviously, solutions are of the form $\beta_+ = n\pi/L$, for integer $n \neq 0$ (similarly for $\beta_-$), and the modal frequencies $s_n$ are, from the solutions Eq. (8) and Eq. (9) of the dispersion relation (4),

$$
s_n = \sigma(n\pi/L) + j\omega(n\pi/L) \quad (15)
$$

over wave numbers for which a traveling solution exists (for
small $b_1$ and $b_2$, this will be true for all $n$).

III. A FINITE DIFFERENCE SCHEME

In order to solve Eq. (6) numerically, we may approxi-
mate its solution over a grid with spacing $X$, and with time
step $T$. Equation (6) can be written as

$$
\delta_t^2 y = c^2 \delta_x^2 y - \kappa^2 \delta_x^2 \delta_y^2 - 2b_1 \delta_{t,0} y + 2b_2 \delta_x^2 \delta_{t,-} y \\
\qquad + O(T, X^2), \tag{16}
$$

where the various difference operators are defined by

$$
\delta_x^2 y(x,t) = \frac{1}{X^2} (y(x+X,t) - 2y(x,t) + y(x-X,t)),
$$

$$
\delta_t^2 y(x,t) = \frac{1}{T^2} (y(x,t+T) - 2y(x,t) + y(x,t-T)),
$$

$$
\delta_{t,0}y(x,t) = \frac{1}{2T}(y(x,t+T) - y(x,t-T)),
$$

$$
\delta_{t,-}y(x,t) = \frac{1}{T}(y(x,t) - y(x,t-T)).
$$

All these operators are “centered” about the point $(x,t)$,
except for the backward difference operator $\delta_{t,-}$, which is
used in order to obtain an explicit algorithm. The approxi-
mation is first-order accurate in the time step $T$, and second-
order accurate in the space step $X$. [It is worth mentioning,
that this is but one among many ways of discretizing (6).]
We may now rewrite Eq. (16) as a difference scheme, oper-
ating on the grid function $y_m^n$, indexed by integer $m$ and $n$,
which will serve as an approximation to $y(x,t)$ at the loca-
tion $x=mX, t=nT$:

$$
y_m^{n+1} = a_{10}y_m^n + a_{11}(y_{m+1}^n + y_{m-1}^n) + a_{12}(y_{m+2}^n + y_{m-2}^n) \\
\qquad + a_{20}y_m^{n-1} + a_{21}(y_{m+1}^{n-1} + y_{m-1}^{n-1}). \tag{17}
$$

Here, the difference scheme coefficients are defined by

$$
\begin{align}
a_{10} &= \frac{(2 - 2\lambda^2 - 6\mu^2 - 4b_2\mu/\kappa)}{(1 + b_1 T)}, \notag \\
a_{11} &= \frac{(\lambda^2 + 4\mu^2 + 2b_2\mu/\kappa)}{(1 + b_1 T)}, \notag \\
a_{12} &= \frac{-\mu^2}{(1 + b_1 T)}, \tag{18}
\end{align}
$$

$$
a_{20} = \frac{(-1 + 4b_2\mu/\kappa + b_1T)}{(1 + b_1T)},
$$

$$
a_{21} = \frac{(-2b_2\mu/\kappa)}{(1 + b_1 T)},
$$

where, for brevity, we have introduced the quantities

$$
\lambda = cT/X, \quad \mu = \kappa T/X^2.
$$

In order to examine the stability of scheme Eq. (17),
which is, like its generating PDE, Eq. (6), linear and shift
invariant, we may apply frequency domain techniques—in the finite difference setting, these techniques are referred to as Von Neumann analysis.⁷,⁸ This analysis proceeds in a fashion exactly analogous to the analysis applied to the continuous time and space systems of the preceding sections. Shortcutting the process somewhat, we may consider a solution of the form $y_m^n = z^n e^{imX\beta}$, where $z = e^{sT}$ (we could equivalently employ a z transform and a discrete-time Fourier transform). We thus obtain the two-step characteristic, or *amplification equation*

$$z^2 + a_1(\beta)z + a_2(\beta) = 0,$$

where the functions $a_1(\beta)$ and $a_2(\beta)$ are defined in terms of the difference scheme parameters of Eq. (18) by

$$a_1(\beta) = -a_{10} - 2a_{11} \cos(\beta X) - 2a_{12} \cos(2\beta X),$$

$$a_2(\beta) = -a_{20} - 2a_{21} \cos(\beta X).$$

The necessary and sufficient stability conditions for Eq. (17) are that the roots of the amplification equation be of magnitude less than or equal to unity, for all wave numbers $\beta$; for this real-valued quadratic, these conditions can be written in terms of $a_1(\beta)$ and $a_2(\beta)$ as

$$|a_1(\beta)| - 1 \le a_2(\beta) \le 1.$$

The right inequality is satisfied for $b_1, b_2 \ge 0$, and after some algebra, it can be shown that the left inequality is equivalent to the condition $\lambda^2 + 4\mu^2 + 4b_2\mu / \kappa \le 1$, further implying that

$$T \le X^2 \frac{-4b_2 + \sqrt{16b_2^2 + 4(c^2X^2 + 4\kappa^2)}}{2(c^2X^2 + 4\kappa^2)}.$$

The relative ease with which an exact bound such as the above may be derived is a direct consequence of the use of a two-step scheme and the relative simplicity of the model itself; for schemes involving more steps of lookback (which results from the discretization of a model with higher time derivatives, such as Ruiz's system), this analysis becomes much more involved, though we do approach it nonetheless in the Appendix. This, in addition to reduced memory requirements, is a further advantage of using a second-order model as a starting point.

Equation (16) is but one of many possible discretizations of Eq. (6)—for instance, replacing $\delta_{t,-}$ by $\delta_{t,0}$ yields an *implicit* algorithm,⁷ and other implicit schemes such as the *θ*-forms discussed in the work of Chaigne⁵ may be of interest in reducing *numerical dispersion*,⁷ and may be of higher formal accuracy (which may be tempered by the stability requirements). To emphasize this point, we have plotted the phase velocities and loss curves for the model system of Eq. (6) versus those of the difference scheme of Eq. (17) in Fig. 1. Notice that this simple difference scheme is a reasonable approximation to the model only for small $\omega$ (i.e., for low frequencies). As we will see later in Sec. V, this deviation from the model PDE will account for differences in simulation results obtained from digital waveguide models, which approximates the phase velocity and loss curves directly.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Phase velocity (top) and loss (bottom) for the model of Eq. (6) (solid line) and for difference scheme Eq. (17) (dashed line) as a function of the frequency. The model parameters are those corresponding to the note C2, as given in Table I.</figcaption></figure>

# IV. THE DIGITAL WAVEGUIDE MODEL

The digital waveguide approach provides computational models for musical instruments primarily in the string, wind, and brass families.¹¹ They have also been developed specifically for piano synthesis.¹¹⁻¹³ This section summarizes the basic ideas of the digital waveguide approach, and relates the parameters of a digital waveguide model to an underlying physical model.

## A. Background

As mentioned in Sec. I, to arrive at a PDE modeling the piano string, it is fruitful to start with the ideal wave equation and add perturbation terms to give more realistic frequency-dependent dispersion and losses. The perturbed PDE is then numerically integrated via a finite-difference scheme (or possibly by another approach, such as finite element methods, etc.). The digital waveguide approach interchanges the order of these operations: the ideal wave equation is integrated first using a trivial finite-difference scheme, and the resulting solutions are perturbed using *digital filters* to add frequency-dependent loss and dispersion. In the case of a strongly dissipative and dispersive string, the modulus of these so-called *loop filters* decreases rapidly with frequency, and phase can become strongly nonlinear. For a frequency-domain implementation, this has no effect on computational complexity, but for a time-domain implementation, a larger filter size may be required in order to match the large variations of the phase response.

It has been known since d'Alembert¹⁴ that the ideal one-dimensional wave equation is solved exactly by arbitrary (sufficiently smooth) wave shapes propagating in both directions. The digital waveguide formulation works directly with these traveling wave components. An isolated traveling wave is trivially simulated in practice using a *delay line*. An ideal vibrating string is then modeled as a pair of delay lines, one for each direction of travel.

For digital implementation, the traveling waves are *sampled* at intervals of $T$ seconds. By Shannon's sampling
<figure><img src="image_4.png" /><figcaption>FIG. 2. Digital waveguide model of a rigidly terminated string.</figcaption></figure>

theorem,¹⁵ the solution remains exact, in principle, at all frequencies up to half the sampling rate. To avoid aliasing, all initial conditions and ongoing excitations must be band-limited to less than half the sampling rate $f_s = 1/T$.

Figure 2 shows the simulation diagram for a digital waveguide model of a rigidly terminated string. The string is excited by the signal E and observed via the signal S. Sampled traveling velocity waves propagate to the right along the upper rail, and to the left along the lower rail. The rigid terminations cause inverting reflections (the two -1 scale factors). The delay lines are denoted $D_i$, $i=1,2,3$, and the $F_i$ blocks are digital filters to be described further below.

Consider a wavelike solution propagating from a point $M_1$ to a point $M_2$ along a string (see Fig. 3, top). The distance $M_1M_2$ will be arbitrarily called $l$ and the propagation time $d$ at the minimal phase velocity. At the observation point $M_2$, the wave will have arrived after having undergone the effects of loss and dispersion. In terms of digital waveguides, the wave will undergo a pure delay [in the frequency domain, a multiplicative phase factor $\exp(-i\omega d)$], times a multiplicative factor $F(\omega)$ representing the loss and the dispersion experienced by the wave during this interval (see Fig. 3, bottom).

Since loss and dispersion are, for this system, linear time-invariant (LTI) phenomena, even when frequency dependent, the perturbations needed for added realism in the digital waveguide string model are LTI digital filters. Since LTI filters commute, we may lump all of the filtering associated with propagation in one direction into a single LTI filter. These filters are denoted $F_i$ in Fig. 2, $i=1,2,3$.

For purposes of computing the output signal S from the input signal E, Fig. 2 may be further simplified.

The two filters labeled $F_1$ can be replaced by a single filter $F_1^2$ (by commuting one of them with the intervening two delay lines $D_1$ and -1 gain). A similar simplification is possible for $F_3$.

In the same way, the two delay lines labeled $D_1$ can be replaced by a single delay line $D_1^2$ (having twice the length of $D_1$), and the two $D_3$ blocks can be replaced by one $D_3^2$ block.

In general, any uniform section of a linear vibrating string which is excited and observed only at its endpoints can be accurately modeled (in one vibrational plane) by a pair of

<figure><img src="image_5.png" /><figcaption>FIG. 3. The physical system and its corresponding delay line/filter.</figcaption></figure>

digital delay lines, each in series with a digital filter. A discussion of more generalized approaches involving nonuniform string sections, and the relationship with finite difference schemes, is provided in a recent dissertation.¹⁶

Since losses and dispersion are relatively weak in vibrating strings and acoustic bores, a low-order filter can approximate very well the distributed filtering (infinite-order in principle) associated with a particular stretch of string or bore. (We should repeat, however, that the approach outlined above is equally applicable in the case of strongly dispersive or lossy systems, though in these cases, higher-order filters may be necessary.) In practice, the desired loss and dispersion filters are normally derived from measurements such as the decay time of overtones in the freely vibrating string.¹⁷⁻¹⁹ In the next section, the filter will be derived from the stiff-string model of Eq. (16). Interestingly, the filter(s) so designed can be mapped back to an equivalent PDE, including many higher-order terms (which may or may not have a physical interpretation). Lumping of traveling-wave filtering in this way can yield computational savings by orders of magnitude relative to more typical finite difference schemes,²⁰⁻²² and this efficiency can be important in applications such as real-time modeling of musical instruments for purposes of automatic sound synthesis or “virtual acoustic instrument” performance.

## B. Relating digital waveguide parameters to the physical model

We address here the problem of relating digital waveguide filter parameters to the loss and dispersion curves from the physical model discussed in Sec. II. For that, we consider the continuous frequency representation of the loop filter and show its relation with the physical parameters. The digital waveguide parameters can then be obtained by discretization. We do not address here the problem of the time domain implementation of the digital waveguide.

According to Eqs. (3) and (7), the transformations of the wave due to propagation along the string segment can be represented in terms of a digital waveguide filter by a multiplicative phase factor $\exp$($sd+j\beta l$). Ideally, the modulus and phase of this expression are related to the filter $F$ by

$$|F(\omega)| = |e^{sd+j\beta l}| = e^{\sigma d} = e^{-(b_1+b_2\beta^2)d}, \quad (19)$$

$$\arg(F(\omega)) = \arg(e^{sd+j\beta l}) = \omega d + \beta l. \quad (20)$$

In order to write the expressions of the modulus and the phase of the loop filter in terms of the frequency $\omega$, it is necessary to express the wave number $\beta$ in terms of $\omega$. From Eq. (9), solving for $\beta$, one gets

$$\beta^2(\omega) = \frac{-\alpha_1 \pm \sqrt{\alpha_1^2 + 4\alpha_2(b_1^2 + \omega^2)}}{2\alpha_2} \quad (21)$$

with

$$\alpha_1 = c^2 - 2b_1b_2, \quad \alpha_2 = \kappa^2 - b_2^2. \quad (22)$$

Since $\beta$ must be real (see Sec. II B), we keep only the solution for which the term inside the root is positive. Then,
$$ \beta(\omega) = \pm \sqrt{\frac{-\alpha_1 + \sqrt{\alpha_1^2 + 4\alpha_2(b_1^2 + \omega^2)}}{2\alpha_2}}. \quad (23) $$

Given that, for realistic piano string modeling, $b_1 \approx 1$, $b_2 \approx 10^{-4}$, $c \approx 200$, $\kappa \approx 1$, and $\omega \approx 400$, we make the simplifying assumptions

$$ b_1 b_2 \ll c^2, \quad b_2^2 \ll \kappa^2 \quad (24) $$

(the second of which also ensures the realness of $\beta$)

$$ b_1^2 \ll \omega^2 \quad (25) $$

which permit the following approximation of $\beta$:

$$ \beta(\omega) \approx \pm \frac{c}{\kappa\sqrt{2}} \sqrt{-1 + \sqrt{1 + 4 \frac{\kappa^2}{c^4} \omega^2}}. \quad (26) $$

In practice, it is helpful to work with more perceptually significant parameters for sound synthesis purposes. For that purpose, we will now suppose that the string is of length $L$, with perfect reflections at the extremities. The delay $D$, which corresponds to the propagation time of the minimum phase-velocity wave along the length $L$ can be expressed as

$$ D = L/c = \pi/\omega_0, \quad (27) $$

where $\omega_0$ is the fundamental frequency (rad/s) of the ideal string and $c$ is the minimum phase velocity given by

$$ c = \omega_0 L / \pi. \quad (28) $$

Thus, $D$ is the propagation time over distance $L$ for a sinusoidal traveling wave tuned to the first resonant mode of the ideal string. This choice of nominal propagation-time $D$ simplifies the frequency-domain computations to follow; however, since phase velocity increases with frequency, the associated propagation filter $F$ will be noncausal in the time domain. This poses no difficulty for frequency-domain implementation.

Next, we may express $\beta$ as

$$ \beta(\omega) \approx \pm \frac{\pi}{L\sqrt{2B}} \sqrt{\xi} \quad (29) $$

with

$$ \xi = -1 + \sqrt{1 + 4B\omega^2/\omega_0^2} \quad (30) $$

in terms of the inharmonicity coefficient$^{23}$ $B$ given by

$$ B = \kappa^2 \omega_0^2 / c^4, \quad (31) $$

where $c$ denotes $c(\omega_0)$ for notational simplicity.

We now have to choose the sign of $\beta$ in the expression for the phase. Since we want the output signal to be *delayed* with respect to the input signal, the loop filter/pure delay combination has to be causal. This means that the phase of the whole transfer function must be negative, i.e.,

$$ -\omega D + \arg(F(\omega)) < 0. \quad (32) $$

This indicates that we choose the negative solution for $\beta$ in the phase expression. Finally, using Eq. (29), we arrive at approximate expressions for the modulus and phase of the filter $F$ as a function of frequency $\omega$,

$$ |F(\omega)| \approx \exp\left(-D\left[b_1 + \frac{b_2 \pi^2 \xi}{2BL^2}\right]\right), \quad (33a) $$

$$ \arg(F(\omega)) \approx \omega D - \pi \sqrt{\frac{\xi}{2B}}. \quad (33b) $$

These expressions serve as the link the PDE model of Eq. (6) to the lumped filters of the digital waveguide. For the sake of simplicity, one can choose an “elementary filter” $\delta F = F^{1/D}$ such that

$$ |\delta F(\omega)| \approx \exp\left(-\left[b_1 + \frac{b_2 \pi^2 \xi}{2BL^2}\right]\right), \quad (34a) $$

$$ \arg(\delta F(\omega)) \approx \omega - \omega_0 \sqrt{\frac{\xi}{2B}}. \quad (34b) $$

The filters $F$ of the digital waveguide, which correspond to propagation over a time duration $D$, can be then easily expressed in terms of $\delta F$ by

$$ F = \delta F^D. \quad (35) $$

The stability of the digital waveguide model is ensured if the modulus of the filter $\delta F$ is less than one. This condition is here always respected, since the expression in the exponential of Eq. (33a) is negative.

## V. NUMERICAL SIMULATIONS

We now address the validity of the analytical approach of the preceding section to determine digital waveguide loop filters. The modulus and phase of the filters corresponding to the digital waveguide model have been directly linked to the parameters of the physical model. For a given set of PDE parameters, thus, we can design a complete digital waveguide model simulating the signal of the vibrating string at a given location, for a predetermined excitation location. For purposes of comparison, we have generated output signals using both the finite difference scheme discussed in Sec. III, and the digital waveguide model for the same set of model parameters. The digital waveguide has been computed in the frequency domain, allowing the use of Eq. (33a). The excitation, a Gaussian function of the form $e^{-(x-x_0)^2/\sigma^2}$, simulates an initial velocity distribution of the string at a distance $x_0=L/8$ from one end of characteristic width $\sigma=K_c L_H$, where $K_c$ is an arbitrary constant and $L_H$ is the hammer width (we do not enter into too much detail here, as hammer modeling is not dealt with in this paper). The signal is observed at the location $9L/10$. We have performed simulations for the notes C2, C4, and C7 using the parameters in Table I (all of which are taken from values provided in the papers by Chaigne and Askenfelt,$^{3,4}$ except for the parameter $b_2$, which comes from the calibration procedure applied subsequently in Sec. VI.

Figure 4 shows the two first periods of the waveforms generated by both approaches. The amplitudes are similar. Nevertheless, there is a slight discrepancy between the two signals, due to the numerical dispersion introduced by the finite difference scheme (see Sec. III).
<figure><img src="image_6.png" /><figcaption>FIG. 4. Velocity signals obtained from the finite difference scheme presented in Sec. III (solid line) and a waveguide model (dotted line), for three different notes and for two periods of sound. The model is excited at distance L/8 from one endpoint, and the output signal is measured at distance 9L/10. Note that the abscissa scales is different for each figure.</figcaption></figure>

This discrepancy can be better seen by comparing the phase velocity of the two signals in Fig. 5. This corresponds to the phase velocity plotted in Fig. 1. The phase velocity of the signal produced by the waveguide is similar to the one of the model.

The long-time behavior of the generated signal is also very similar. In Fig. 6 spectrograms obtained over the whole length of the sound are shown. It is clear that the global damping behavior is similar in the low-frequency range. However, at high frequencies, the finite-difference model suffers from an artificially high propagation gain, as derived in Sec. III. The fundamental frequencies are essentially

<figure><img src="image_7.png" /><figcaption>FIG. 5. Phase velocity for the waveguide model (solid line) and for difference scheme as a function of the frequency for the note C2. The analytically obtained phase is not shown, as it is identical (by definition) to the phase response of the waveguide network.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 6. Spectrograms of the output of the finite difference scheme (at top) and the waveguide model (at bottom) for the note C2.</figcaption></figure>

equal, but the wave dispersion due to string stiffness is different due to the *numerical* dispersion introduced by the difference scheme. (We can see in Fig. 7 a slight tuning difference in the high-frequency partials.) In summary, these figures illustrate the extent to which the waveguide model provides a more accurate digital simulation of stiff, lossy strings with respect to both attenuation and dispersion of wave propagation, when compared with finite difference schemes.

# VI. CALIBRATION OF PHYSICAL MODEL PARAMETERS FOR A GRAND PIANO FROM EXPERIMENTAL DATA

## A. Experimental setup

In order to calibrate the parameters of the physical model, data were collected on an experimental setup consisting of a Yamaha Disklavier C6 grand piano equipped with sensors (see Fig. 8). The string vibrations were measured

<figure><img src="image_9.png" /><figcaption>FIG. 7. Partial frequencies of the output of the finite difference scheme (dotted) and of the waveguide model (plain) as function of the theoretical frequencies.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Experimental setup. The grand piano was isolated in an anechoic room, and both the string vibrations and the acoustical radiated signal were measured. The string vibrations were measured using both an accelerometer located at the bridge level and a laser vibrometer, while the acoustic signal was measured at the head level of the pianist using an artificial headset. Our library of measured data also includes signals corresponding to various hammer velocities [referred to in recent and forthcoming articles concerned with the hammer-string interaction (Ref. 31)], but for this paper, we only need acceleration measurements for each string.</figcaption></figure>

using an accelerometer located at the bridge level. The Disklavier allows the piano to be played under computer control. These measurements were made at the Laboratoire de Mécanique et d'Acoustique, Marseille, France.

In order to ensure that the measurements were taken under similar excitation conditions, we measured the hammer velocity using a photonic sensor (MTI 2000, probe module 2125H) (see Fig. 9).

Since we were interested in exciting a large portion of the frequency spectrum while remaining in the linear domain of vibration, we chose a medium (mezzo-forte) hammer velocity of 2.2 m s⁻¹, which corresponds roughly to a MIDI value of 80. Such a hammer velocity allows the generation of about 140 spectral components for low frequency tones with a reasonable signal-to-noise ratio (see Fig. 10).

For each note of the grand piano, the optical sensor was placed close to the hammer and the accelerometer (B&K

<figure><img src="image_11.png" /><figcaption>FIG. 9. Optical sensor used to measure the hammer velocity. The velocity is obtained through the duration corresponding to the travel time of the hammer between two reflectors placed on it.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 10. Spectrum of the note A0.</figcaption></figure>

4374) at bridge level. For notes corresponding to double or triple sets of strings, the accelerometer was placed as close as possible to the strings (see Fig. 11). Due to the imprecision of MIDI coding, several measurements were made, until a target value of the hammer speed was obtained. We have deemed an uncertainty of ±0.1 m s⁻¹ for the hammer velocity to be acceptable, as the estimation of modal frequency and decay rates is relatively insensitive to such an error. The acceleration measured at the bridge level was then digitally recorded at 16 bits at a sampling rate of 44.1 kHz, before being entered in the database.

## B. Estimation of parameters

Because this model is intended for use in the context of musical sound synthesis, we here discuss the calibration of $b_1$ and $b_2$, and the determination of the stiffness parameter for a given string.

<figure><img src="image_13.png" /><figcaption>FIG. 11. Accelerometer at the bridge level.</figcaption></figure>
To estimate the damping factor associated with each modal component of the signal, we used a signal processing technique based on the theory of analytic signals.²⁴,²⁵ The analytic signal representation provides an easy way of extracting both the instantaneous frequency and the damping of each modal component, through band-pass filtering. To isolate each component in frequency, we used a truncated Gaussian window, the frequency bandwidth of which was chosen so as to minimize smoothing effects over the attack duration and to avoid overlapping two successive frequency components. The Gaussian was employed since its time-bandwidth product is minimized. As a consequence, it optimizes the exponential damping support after convolution with a causal single component for a given band pass filtering.²⁶

The analytic signal associated with each individual component facilitates the estimation of both the instantaneous frequency and the amplitude modulation law of the component. The frequency dependent damping factor is directly related to the amplitude modulation law of each partial. According to the physics of a single string vibrating in one plane, the amplitude modulation of each component is expected to be exponential. This makes the damping factor easy to estimate by the measurement of the slope of the logarithmic representation of the amplitude of the analytic signal.¹⁹ This technique has advantages over other parametric methods such as Prony's method,²⁷ mainly due to its ability to extract a coherent mean damping factor when multiple components are present. In fact, the hammer usually strikes not one, but two or three strings simultaneously. The coupling gives rise to perceptually significant phenomena such as beating and two-stage decay;²⁸ these effects are not accounted for in our model Eq. (6). For these multistring notes, the calculated damping coefficients σ can be thought of describing the global perceived decay of the sound.

For the same reason, the spectral representation of each partial is the result of the summation of several contributions due to each string and each polarization of vibration. The phase of the analytic signal allows an accurate estimation of the mean frequency of the partials. Actually, it permits the calculation of the instantaneous frequency which is a time-dependent function oscillating around a mean value. This mean value coincides with the spectral centroid of the partial²⁹ and consequently with the more likely perceived frequency. From the mean frequency values estimated this way for each partial, B may be deduced for each note, and consequently ω₀, the fundamental frequency of the corresponding ideal string. The inharmonicity factor B is plotted as a function of the frequency in Fig. 12; B is an increasing function of the note number, except over the bass range, where the strings are double-wrapped (this behavior has also been measured by Conklin³⁰).

Combining Eqs. (8) and (29), one obtains

$$ \sigma(\omega) = -b_1 - b_2 \left( \frac{\pi^2}{2BL^2} \left[ -1 + \sqrt{1 + 4B(\omega/\omega_0)^2} \right] \right). \quad (36) $$

Then, from the value of σ obtained for each partial, b₁ and b₂ may be estimated for a given tone. The evolution of these

<figure><img src="image_14.png" /><figcaption>FIG. 12. The measured inharmonicity factor B.</figcaption></figure>

parameters as a function of the fundamental frequency is shown in Fig. 13.

We see that b₁ and b₂ are both increasing functions of MIDI note number, indicating increasing loss as one approaches the treble range. In Fig. 13, we have also fit extremely simple curves to the loss parameter data. The fits are linear as a function of the fundamental frequency, and are given by

$$ b_1 = 4.4 \times 10^{-3} f_0 - 4 \times 10^{-2}, \quad (37a) $$

$$ b_2 = 1.0 \times 10^{-6} f_0 + 1 \times 10^{-5}. \quad (37b) $$

These simple empirical descriptions of b₁ and b₂ allow the reproduction of piano tones whose damping will be close to that of the perceived acoustic note. A multistring waveguide

<figure><img src="image_15.png" /><figcaption>FIG. 13. Values of “equivalent” b₁ and b₂ fitted from measured data as a function of the fundamental frequency.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 14. Normalized modulus and phase of F for selected tones.</figcaption></figure>

model has also been designed,²⁶,³¹ allowing for beating and two-stage decay. Its use for synthesis purposes is discussed in Sec. VII. The deviations in the curves (Fig. 13) from the linear fits must be attributed to the impedance ratio between the strings and the soundboard, which varies over the length of the bridge. As a result, the curves should not be interpreted as impedance curves as they do not represent measurements at a single point and data are taken only at string modal frequencies. For a detailed discussion of soundboard impedance and its measurement, we refer to the work of Giordano.³²,³³

## VII. WAVEGUIDES AND SOUND SYNTHESIS

The determination of $b_1$, $b_2$, B, and $\omega_0$ for each note allows for an explicit expression of the behavior of the filter F as a function of note number. In order to represent the evolution of the elementary loop filter in terms of note value, we show in Fig. 14 the modulus and phase of the elementary filter $\delta F$, normalized with respect to the time delay D,

$$ \delta F = F^{1/D}. \qquad (38) $$

In order to understand the general behavior of the modulus and the phase of the loop filter with respect to the note played, we expand their expressions [Eq. (34a)] for $4B(\omega/\omega_0)$ to third order near zero. We obtain

$$ |\delta F(\omega)| \approx \exp\left(-\left[b_1 + b_2 \frac{\omega^2}{c^2}\right]\right), \qquad (39a) $$

$$ \arg(\delta F) \approx \frac{\omega^3 B}{2 \omega_0^2}. \qquad (39b) $$

The modulus (which also accounts for the losses at the endpoints), is decreasing with note number as shown in Fig. 14. This is mainly due to the increase in $b_1$. But the parameter $b_2$, which allows for frequency dependent loss is also

<figure><img src="image_17.png" /><figcaption>FIG. 15. Amplitude of partials one to six for the note E1 as a function of time and frequency.</figcaption></figure>

increasing, leading to a modification of the slope of the modulus versus frequency for different note numbers. We note that this behavior is slightly different for the wrapped strings (A0, A1) than for the other strings (A2, A3). We also note that although B is mainly an increasing function of note number, the phase of the filters $\delta F$ grows less rapidly for the bass tones than for higher tones. This is due to the fact that the phase of the filter depends not only on B, but also on the fundamental frequency, as evidenced by Eq. (39b). Though this expression is meaningful only for the first few partials, it is clear that phase dispersion decreases as a function of fundamental frequency.

In the case of the piano, the string is struck at a distance of approximately one-eighth to one-sixteenth of its length from the bridge, depending on the note. We are only interested in the vibration generated at the bridge termination, as this is the mechanism by which energy is transmitted to the soundboard. This situation corresponds to a digital waveguide structure identical to the one presented in Fig. 2, except that the loop filters $F_2$ and $F_3$ are combined. Using the parameters estimated by experiment, one can reproduce the vibration generated by a single string. Figure 15 shows the evolution of the amplitude of the first six partials of the vibration velocity for the note E1, respectively, measured on the piano and generated by the digital waveguide model. From a perceptual point of view, the two sounds are identical.

If the tones are produced by two or three strings struck simultaneously, the basic digital waveguide model still generates a signal having the same amplitude and damping of the modes. It does not, however, account for the modulations and double decays due to the coupling of the strings at the bridge.²⁸ Using two or three coupled digital waveguides, and thus allowing for energy transfer between the strings, one can easily reproduce this phenomenon.²⁰,²⁶,³¹ Figure 16 shows the time evolution of the amplitude of the first six partials for the tone C2 (two strings), using the coupled-
<figure><img src="image_18.png" /><figcaption>FIG. 16. Amplitude of partials one to six for the note C2 as a function of time and frequency.</figcaption></figure>

strings model described in a recent publication.²⁶ The modulations are essentially perfectly reproduced. Moreover, this model follows directly from the physics of coupled strings. In fact, the loop filters are again related to the parameters of the physical model, and numerous sound transformations are conceivable. One could, for instance, extend the use of such a waveguide model to practically unrealizable situations involving, for example, widely mistuned strings or coupled strings with differing material properties.

## VIII. CONCLUSIONS

We have presented a model of transverse vibrations on a string which includes the effects of stiffness and frequency-dependent loss. This model possesses several advantages over those proposed previously, in particular that it can be framed as a well-posed initial-boundary value problem (leading to stable numerical methods), and also that it can be easily related to digital waveguides. The source of these good properties is the fact that this model allows only two traveling-wave type solutions; the nonphysical third unstable term in the model of Ruiz (which can lead to difficulties both analytically and numerically) is thus eliminated in favor of higher-order spatial terms. It is also simple to write expressions for dispersion and loss as a function of frequency in terms of the model system parameters—such information is critical for the design of the terminating filters in a digital waveguide implementation.

For the sake of comparison, we have performed numerical simulations of the model system, with pianolike parameters, using both finite differences and a digital waveguide; the most notable distinction is the complete lack of numerical dispersion (which leads to mode mistuning) in the waveguide implementation. On the other hand, the waveguide allows the computation of a solution (“sound”) only at preselected points on the string, whereas a finite difference scheme computes the entire string state in (sampled) physical

form. This is not a drawback for sound synthesis applications because, only the behavior of the string at the bridge is of interest in most stringed instruments. Moreover, physically accurate outputs from additional points along the string are easily added to a digital-waveguide simulation at the price of one small digital filter each.

A set of experimental data measured from a grand piano was used in order to calibrate the PDE model parameters over the entire keyboard range. String vibration was measured at the bridge through the use of an accelerometer, for each note on the piano, and for an average hammer velocity. The piano employed was equipped with sensors to provide hammer velocity data; from these measurements, all the parameters relating to the relevant PDE model were estimated. Given that the model itself is not completely physical—that is, the various loss mechanisms, interstring coupling, as well as energy transfer to the soundboard are modelled, for simplicity, as internal to the string itself—these parameters must be considered as those describing an “equivalent” string, under fixed termination. The equivalent parameters, are, however, sufficient for the resynthesis of piano tones to a high degree of fidelity, when a digital waveguide is employed. The digital waveguide model was also extended in order to directly take into account the effects of interstring coupling, through the use of two or three coupled waveguides.

The modeling of the excitation mechanism for the piano string (i.e., the hammer) is also of great importance, and must be carried out with some care; we have not addressed this issue here. As has been shown in the work of Chaigne and Askenfelt, it is possible to design a nonlinear hammer, which, when applied to a stiff string with frequency-dependent loss, produces signals quite similar to those measured on a real piano. The problem of extracting hammer parameters from measured data is also worthy of future research.

## ACKNOWLEDGMENTS

This work was supported in part by the DSP (France) under Contract No. 016060 and by a fellowship from Region PACA (France).

## APPENDIX: THE PIANO STRING MODEL OF CHAIGNE AND ASKENFELT

The results in this section have appeared, in a similar form, in the thesis of Ruiz.⁹ We have added various comments regarding well-posedness and numerical stability.

The stiff string model in the thesis of Ruiz⁹ and in the papers by Chaigne and Askenfelt³,⁵ is described by the following equation:

$$ \frac{\partial^2 y}{\partial t^2} = c^2 \frac{\partial^2 y}{\partial x^2} - \kappa^2 \frac{\partial^4 y}{\partial x^4} - 2b_1 \frac{\partial y}{\partial t} + 2b_3 \frac{\partial^3 y}{\partial t^3}. \quad (A1) $$

This model differs from Eq. (6) only by the replacement of the term $2b_2(\partial^3 y/\partial x^2 \partial t)$ by $2b_3(\partial^3 y/\partial t^3)$; this model also allows for frequency-dependent loss, but the system itself is of a quite different character, due to the increased degree of the equation with respect to the time variable. We spend a
little time here explaining the significance of the difference, which has a radical effect on the analysis of the system as a whole.

We can examine the well posedness of the system by inserting a solution of the form $y(x,t)=e^{st+j\beta x}$ into Eq. (A1), in order to obtain a dispersion relation,

$$-2b_3s^3+s^2+2b_1s+c^2\beta^2+\kappa^2\beta^4=0. \quad (A2)$$

This is a cubic in the variable $s$ [the quantity on the left-hand side is referred to as the symbol⁷ of Eq. (A1)], and again, as discussed in Sec. II, a necessary condition for the system of Eq. (A1) to be well posed is that the roots of this equation have real part bounded from above as a function of spatial frequency $\beta$. It is simple to see that the real part of at least one root of Eq. (A2) will be positive and unbounded as a function of wave number. As this is a third-degree polynomial equation with real coefficients, one root will always be real, and the two others occur as a complex conjugate pair (or perhaps as two other real roots). Consider Eq. (A2) as $|\beta|$ becomes large. In this case, the three roots will behave as

$$-2b_3s^3+\kappa^2\beta^4\approx 0$$

and will be evenly distributed over a circle of radius $(\kappa^2\beta^4/(2b_3))^{1/3}$. If $b_3>0$ (as suggested in the numerical experiments in the papers by Chaigne and Askenfelt⁴), then there will be one positive real root of the magnitude mentioned above, clearly unbounded as a function of wave number $\beta$. (If $b_3<0$, there will be two roots in the right half plane, of this same magnitude, at an angle approaching $\pm 60$ degrees with respect to the positive real axis.) We have thus shown that the initial value problem corresponding to the system of Eq. (A1) is, formally speaking, ill-posed.

We can extract some more detailed information by asking under what conditions the roots of Eq. (A2) have positive real part. A straightforward application of the Routh–Hurwitz stability criterion³⁴ to Eq. (A2) shows that in fact, if $b_3>0$, there is always exactly one real positive root, regardless of the value of the wave number $\beta$.

The following question then arises: How can we explain the apparently stable behavior exhibited by simulations⁴ of these equations? Indeed, for realistic piano string parameters, the numerical integration routine provided in the paper by Chaigne and Askenfelt⁴ is stable, and produces piano sounds of excellent quality. A first guess might be that the above analysis is incomplete due to the neglect of boundary conditions. Using the boundary conditions supplied by Chaigne, however, leads to an analysis identical to that performed in Sec. II B—the modal frequencies for the string system defined by Eq. (A1) will be given by solutions of the dispersion relation Eq. (A2) under the replacement of $\beta$ by $n\pi/L$ for integer $n$. For any $n$, there will be exactly one modal frequency $s_n$ with positive real part. Thus the instability persists even in the presence of boundary conditions.

We must then conclude that discretization has a stabilizing effect on system of Eq. (A1). To explore this idea in more detail, consider the discretization,⁴ which can be written as

$$\begin{aligned} y_m^{n+1} = & a_{10}y_m^n + a_{20}y_m^{n-1} + a_{11}(y_{m+1}^n + y_{m-1}^n) + a_{12}(y_{m+2}^n \\ & + y_{m-2}^n) + a_{21}(y_{m+1}^{n-1} + y_{m-1}^{n-1}) + a_{30}y_m^{n-2}. \end{aligned} \quad (A3)$$

This difference scheme involves three steps of lookback, reflecting the degree of the model system of Eq. (A1). Here, the difference scheme coefficients are defined by

$$\begin{aligned} a_{10} &= (2 - 2\lambda^2 - 6\mu^2 + b_3/T)/D, \\ a_{20} &= (-1 + b_1T + 2b_3/T)/D, \quad a_{11} = (\lambda^2 + 4\mu^2)/D, \\ a_{12} &= (b_3/T - \mu^2)/D, \quad a_{21} = a_{30} = (-b_3/T)/D, \end{aligned} \quad (A4)$$

where, for brevity, we have again used

$$\lambda = cT/X, \quad \mu = \kappa T/X^2,$$

as well as

$$D = 1 + b_1T + 2b_3/T.$$

Let us now examine the characteristic polynomial, which can be written as

$$z^3 + a_1(\beta)z^2 + a_2(\beta)z + a_3(\beta) = 0 \quad (A5)$$

with

$$\begin{aligned} a_1(\beta) &= -a_{10} - 2a_{11}\cos(\beta X) - 2a_{12}\cos(2\beta X), \\ a_2(\beta) &= -a_{20} - 2a_{21}\cos(\beta X), \quad a_3(\beta) = -a_{30}. \end{aligned}$$

The solution to the recursion will be bounded and decay if the solutions to this equation are confined to the interior of the unit circle for all $\beta \in [-\pi/X, \pi/X]$. It is simple to show that this is in fact true, for any of the choices of parameters given in the papers by Chaigne and Askenfelt.⁴ This does not mean, however, that the difference scheme can be considered to be numerically stable in the Von Neumann sense.⁷ This is a rather subtle point, and is worth elaborating.

According to the Lax–Richtmeyer equivalence theorem,⁷ if the initial-boundary value problem is well-posed, the solution to a finite difference scheme will converge to the solution of the model problem if it is consistent and stable. In this case, though, the model system is not well posed, and thus no finite difference can possibly converge to a stable solution in some limit as the time step $T$ and the grid spacing $X$ approach zero. The difference scheme Eq. (A3) is indeed consistent with (A1) to first-order accuracy (and we note that if one does wish to use this ill-posed model system, it is in fact possible to design second-order accurate explicit methods), but it is possible to show (as we expect) that it cannot be stable in the limit as $T$ becomes small. Because the recursion is of third order, the analysis is somewhat involved, and requires the application of the Schur–Cohn recursive procedure⁷ (the discrete time analog of the Routh–Hurwitz stability test, again allowing us to check the stability of a polynomial without explicitly calculating the roots). Nevertheless, it is possible to show in this way that a necessary condition that the roots of the polynomial of Eq. (A5) be inside the unit circle is that

$$b_3/T \le \lambda^2 + 4\mu^2 \le 1.$$

(The second inequality is exactly the necessary stability condition given in the paper by Chaigne and Askenfelt.³) Al-
though we have not provided all the details, we note that it
suffices to perform the Schur–Cohn test at the wave number
$\beta = \pi/X$ in order to arrive at these conditions. Clearly, these
two conditions cannot be satisfied if

$$
T \le b_3 \Rightarrow \text{Sample rate} \ge 1/b_3
$$

and thus for a small enough time step, the system poles must
cross to the exterior of the unit circle, *regardless of the grid
spacing X*.

It should be said, however, that because $b_3$ is in general
quite small for realistic piano string models (on the order of
$10^{-9}$), for any reasonable sample rate in the audio range, the
recursion does not exhibit this unbounded growth. On the
other hand, as we have shown in this paper, it is simple
enough to dispense with the nonphysical solution and all the
concomitant analysis by making use of a simpler second-
order model.

<footnote n="1">¹L. Hiller and P. Ruiz, “Synthesizing musical sounds by solving the wave equation for vibrating objects,” J. Audio Eng. Soc. **19**, 462–551 (1971).</footnote>

<footnote n="2">²R. Bacon and J. Bowsher, "A discrete model of a struck string," Acoustica **41**, 21–27 (1978).</footnote>

<footnote n="3">³A. Chaigne and A. Askenfelt, “Numerical simulations of struck strings. I. A physical model for a struck string using finite difference methods,” J. Acoust. Soc. Am. **95**, 1112–1118 (1994).</footnote>

<footnote n="4">⁴A. Chaigne and A. Askenfelt, “Numerical simulations of struck strings. II. Comparisons with measurements and systematic exploration of some hammer-string parameters,” J. Acoust. Soc. Am. **95**, 1631–1640 (1994).</footnote>

<footnote n="5">⁵A. Chaigne, “On the use of finite differences for musical synthesis. Application to plucked stringed instruments,” J. Acoust. **5**, 181–211 (1992).</footnote>

<footnote n="6">⁶N. Fletcher and T. Rossing, *The Physics of Musical Instruments* (Springer-Verlag, New York, 1991).</footnote>

<footnote n="7">⁷J. Strikwerda, *Finite Difference Schemes and Partial Differential Equations* (Wadsworth and Brooks/Cole Advanced Books and Software, Pacific Grove, CA, 1989).</footnote>

<footnote n="8">⁸B. Gustaffson, H.-O. Kreiss, and J. Oliger, *Time Dependent Problems and Difference Methods* (Wiley, New York, 1995).</footnote>

<footnote n="9">⁹P. M. Ruiz, "A technique for simulating the vibrations of strings with a digital computer," Ph.D. thesis, University of Illinois, 1970.</footnote>

<footnote n="10">¹⁰K. Graff, *Wave Motion in Elastic Solids*, 3rd ed. (Prentice-Hall, Englewood Cliffs, NJ, 1974).</footnote>

<footnote n="11">¹¹J. O. Smith III, *Digital Waveguide Modeling of Musical Instruments*, http://www.ccrma.stanford.edu/~jos/waveguide/http://www-ccrma.stanford.edu/~jos/waveguide/, 2002 (printed version forthcoming).</footnote>

<footnote n="12">¹²J. O. Smith III and S. A. Van Duyne, “Committed piano synthesis,” *Proceedings of the 1995 International Computer Music Conference, Banff (Computer Music Association, 1995)*, pp. 319–326.</footnote>

<footnote n="13">¹³S. A. Van Duyne and J. O. Smith III, “Developments for the commuted piano,” *Proceedings of the 1995 International Computer Music Conference, Banff (Computer Music Association, 1995)*, pp. 335–343.</footnote>

<footnote n="14">¹⁴J. le Rond d'Alembert, "Investigation of the curve formed by a vibrating string, 1747," in *Acoustics: Historical and Philosophical Development*, edited by R. Bruce Lindsay (Dowden, Hutchinson & Ross, Stroudsburg, 1973), pp. 119–123.</footnote>

<footnote n="15">¹⁵J. O. Smith III and P. Gossett, “A flexible sampling-rate conversion method,” *Proceedings of the International Conference on Acoustics, Speech, and Signal Processing, San Diego* (IEEE Press, New York, 1984), Vol. 2, pp. 19.4.1–19.4.2. An expanded tutorial based on this paper and associated free software are available online at The Digital Audio Resampling Home Page: http://www.ccrma.stanford.edu/~jos/resample/http://www-ccrma.stanford.edu/~jos/resample/</footnote>

<footnote n="16">¹⁶S. Bilbao, Wave and scattering methods for the numerical integration of</footnote>

<footnote>partial differential equations, Ph.D. thesis, Stanford University, 2001.</footnote>

<footnote n="17">¹⁷J. O. Smith III, “Techniques for digital filter design and system identification with application to the violin,” Ph.D. thesis, Stanford University (CCRMA), 1983. Available as CCRMA Technical Report STAN-M-14. Portions available online at http://www.ccrma.stanford.edu/~jos/http://www-ccrma.stanford.edu/~jos/</footnote>

<footnote n="18">¹⁸J. Laroche and J. L. Meillier, “Multichannel excitation/filter modeling of percussive sounds with application to the piano,” IEEE Trans. Speech Audio Process. **2**, 329–344 (1994).</footnote>

<footnote n="19">¹⁹V. Välimäki, J. Huopaniemi, M. Karjalainen, and Zoltan Jánosy, “Physical modeling of plucked string instruments with application to real-time sound synthesis,” J. Audio Eng. Soc. **44**, 331–353 (1996).</footnote>

<footnote n="20">²⁰J. O. Smith III, “Efficient synthesis of stringed musical instruments,” *Proceedings of the 1993 International Computer Music Conference, Tokyo (Computer Music Association, 1993)*, pp. 64–71. Available online at http://www.ccrma.stanford.edu/~jos/cs/</footnote>

<footnote n="21">²¹M. Karjalainen, V. Välimäki, and Z. Jánosy, “Towards high-quality sound synthesis of the guitar and string instruments,” *Proceedings of the 1993 International Computer Music Conference, Tokyo (Computer Music Association, 1993)*, pp. 56–63. Available online at http://www.acoustics.hut.fi/~vpv/publications/icmc93-guitar.htm</footnote>

<footnote n="22">²²T. Hikichi and N. Osaka, "An approach to sound morphing based on physical modeling," *Proceedings of the International Computer Music Conference* (Beijing, Computer Music Association, 1999), pp. 108–111.</footnote>

<footnote n="23">²³H. Fletcher, E. D. Blackham, and R. Stratton, "Quality of piano tones," J. Acoust. Soc. Am. **34**, 749–761 (1961).</footnote>

<footnote n="24">²⁴D. Gabor, “Theory of communication,” J. Inst. Electr. Eng., Part 1 **93**, 429–457 (1946).</footnote>

<footnote n="25">²⁵R. Kronland-Martinet, Ph. Guillemain, and S. Ystad, “Modelling of natural sounds using time-frequency and wavelet representations,” Organised Sound **2**, 179–191 (1997).</footnote>

<footnote n="26">²⁶M. Aramaki, J. Bensa, L. Daudet, P. Guillemain, and R. Kronland-Martinet, “Resynthesis of coupled piano string vibrations based on physical modeling,” J. New Music Res. **30**, 213–226 (2002).</footnote>

<footnote n="27">²⁷J. D. Markel and A. H. Gray, *Linear Prediction of Speech* (Springer-Verlag, New York, 1976).</footnote>

<footnote n="28">²⁸G. Weinreich, “Coupled piano strings,” J. Acoust. Soc. Am. **62**, 1474–1484 (1977). Also contained in Ref. 36. See also Sci. Am. **240**, 94 (1979).</footnote>

<footnote n="29">²⁹Ph. Guillemain and R. Kronland-Martinet, “Characterization of acoustic signals through continuous linear time-frequency representations,” IEEE Special Issue on Time-Frequency and Time-Scale Analysis, **84**, 561–585 (1996).</footnote>

<footnote n="30">³⁰H. A. Conklin, Jr., "Design and tone in the mechanoacoustic piano. III. Piano strings and scale design," J. Acoust. Soc. Am. **100**, 1286–1298 (1996).</footnote>

<footnote n="31">³¹J. Bensa, F. Gibaudan, K. Jensen, and R. Kronland-Martinet, “Note and hammer velocity dependence of a piano string model based on coupled digital waveguides,” *Proceedings of the International Computer Music Conference, La Havana, Cuba* (Computer Music Association, 2001), pp. 95–98.</footnote>

<footnote n="32">³²N. Giordano, “Mechanical impedance of a piano soundboard,” J. Acoust. Soc. Am. **103**, 2128–2133 (1998).</footnote>

<footnote n="33">³³N. Giordano, “Sound production by a vibrating piano soundboard: Experiment,” J. Acoust. Soc. Am. **104**, 1648–1653 (1998).</footnote>

<footnote n="34">³⁴M. Van Valkenburg, *Network Analysis* (Dover, New York, 1975).</footnote>

<footnote n="35">³⁵J. O. Smith III, Music applications of digital waveguides. Technical Report STAN-M-39, CCRMA, Music Dept., Stanford University, 1987. A compendium containing four related papers and presentation overheads on digital waveguide reverberation, synthesis, and filtering. CCRMA technical reports can be ordered by calling (650)723-4971 or by sending electronic mail requests to info@ccrma.stanford.edu</footnote>

<footnote n="36">³⁶*Five Lectures on the Acoustics of the Piano*, edited by A. Askenfelt (Royal Swedish Academy of Music, Stockholm, 1990). Lectures by H. A. Conklin, A. Askenfelt, E. Jansson, D. E. Hall, G. Weinreich, and K. Wogram. Sound example CD included. Publication number 64. Available online at http://www.speech.kth.se/music/5_lectures/</footnote>