# Parameter fitting for piano sound synthesizer based on physical and perceptual attributes

Julien Bensa, Olivier Gipouloux, Richard Kronland-Martinet

► To cite this version:

Julien Bensa, Olivier Gipouloux, Richard Kronland-Martinet. Parameter fitting for piano sound synthesizer based on physical and perceptual attributes. Journal of the Acoustical Society of America, 2005, 118 (1), pp.495-504. 10.1121/1.1929230. hal-00105865

HAL Id: hal-00105865

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

Submitted on 4 Sep 2020

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Parameter fitting for piano sound synthesis by physical modeling

Julien Bensa,ᵃ) Olivier Gipouloux,ᵇ) and Richard Kronland-Martinet

Laboratoire de Mécanique et d'Acoustique, 31 Chemin Joseph Aiguier, 13402 Marseille, France

A difficult issue in the synthesis of piano tones by physical models is to choose the values of the parameters governing the hammer-string model. In fact, these parameters are hard to estimate from static measurements, causing the synthesis sounds to be unrealistic. An original approach that estimates the parameters of a piano model, from the measurement of the string vibration, by minimizing a perceptual criterion is proposed. The minimization process that was used is a combination of a gradient method and a simulated annealing algorithm, in order to avoid convergence problems in case of multiple local minima. The criterion, based on the tristimulus concept, takes into account the spectral energy density in three bands, each allowing particular parameters to be estimated. The optimization process has been run on signals measured on an experimental setup. The parameters thus estimated provided a better sound quality than the one obtained using a global energetic criterion. Both the sound's attack and its brightness were better preserved. This quality gain was obtained for parameter values very close to the initial ones, showing that only slight deviations are necessary to make synthetic sounds closer to the real ones.

## I. SCOPE OF THE PROBLEM

Even though the synthesis of piano tones by physical modeling has been widely addressed, musicians call for further improvements. Actually, the specificity and the complexity of the instrument lead to a unique timbre dynamically linked to the playing. As a consequence, listeners have developed cognitive knowledge of piano tones that allows them to evaluate sound quality by taking into account various aspects and subtleties of the sound.¹ Among numerous perceptually important contributions to piano tones, one can cite for example the inharmonicity of the partials, the frequency-dependent damping, the beating due to the coupling of several strings (possibly attached to different notes), and the relations between the hammer velocity and the sound brightness. All these “sound ingredients” are linked to the physics of the instrument, making a physical approach more adapted to the synthesis if the physical parameters are precisely known. However, the perceptual influence of these ingredients is far from obvious, making it even more difficult to adjust these parameters when one seeks to minimize the perceptual difference between a synthesis sound and an original piano sound. In this paper, we estimate the parameters of a physical model that describes the phenomenon occurring during and after the hammer-string interaction. This estimation is based on the minimization of a perceptual criterion measuring the difference between a recorded vibration of the string and a synthesized one. This approach allows one to use a single signal to estimate the parameters governing the hammer-string interaction, which are usually measured in

static.² As already mentioned, an accurate estimation of those particular parameters is extremely useful (and difficult under dynamic conditions) since the hammer-string interaction determines a perceptually important part of the piano sound (especially the dynamic of the sound). Hence, this study constitutes a crucial step toward the design of realistic synthetic sounds. The estimation method is schematized in Fig. 1 and is an improved version of the method presented in Ref. 3. Starting from initial parameters, the optimization process computes a new set of parameters through the model algorithm and the perceptual criterion. This criterion measures perceptual differences between the synthetic signal given by the model and the original signal recorded on an experimental setup. The “optimized” parameters obtained provide a synthetic signal perceptually closer to the original.

In Sec. II, we describe the physical model, exhibiting the parameters that have to be estimated. A model was chosen that takes into account the hammer-string interaction and the string vibration. It does not constitute a full piano model since there is no soundboard and coupling between strings, for instance. However, a piano model taking into account all the elements in the sound production would be extremely complex and would exhibit many parameters to be taken into account in the optimization process. Our estimation method is presented for a simple case, which may be, in the future, applied to models that are more complex. In Sec. III, we address the direct problem, which consists in the study of the string behavior with respect to the model parameters. This aspect is mainly based on a numerical calculus of the solutions provided by the model. Some critical aspects of the physical behavior are also described, such as the hammer force behavior and the felt compression's influence. Section IV is devoted to the inverse problem consisting in estimating the parameters of the mechanical model by using reference

<footnote>ᵃ) Electronic mail: bensa@lma.cnrs-mrs.fr</footnote>

<footnote>ᵇ) Also at: Modélisation et Calcul Scientifique, Université Lyon 1, Bâtiment ISTIL, 15 Boulevard Latarjet 69622 Villeurbanne Cedex, France.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Principal stages of the calibration process.</figcaption></figure>

signals recorded on the experimental setup. This setup measures vibrations from one isolated string struck by a piano hammer, and of the corresponding hammer velocities. The setup is suitable for our model since there is no soundboard and string coupling. The use of real piano sounds would require using another model, since the damping and the frequencies of the partials may be modified by the finite impedance at the bridge and by the coupling between several strings. The collected data are further computed using an optimization technique that combines gradient and simulated annealing methods minimizing perceptual criteria. These criteria are based on the spectral energy in three frequency bands, analogous to the so-called perceptual tristimulus decomposition.⁴ In Sec. V, the estimated parameters are used to validate the method through sound synthesis.

## II. STRING AND HAMMER MODEL

Several models of piano-string and hammer-string interaction have been proposed. They generally consist of a set of PDEs governing the string motion⁵⁻⁷ and the hammer displacement that are coupled by a power law describing the compression and losses of the felt as a function of the hammer and string displacement.²,⁸⁻¹⁰ The model of wave propagation on the string we used was proposed in Ref. 7 and is given by

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

with $y_s$ the string displacement, $c$ the wave speed, $\kappa$ the stiffness coefficient, $b_1$ and $b_2$ the loss coefficients (a symbols table may be found in Table I). The first term on the right-hand side gives rise to wave-like motion with speed $c$. The second term introduces dispersion, and the last two terms introduce frequency-dependent loss via mixed time-space derivative terms. This model differs from that of Hiller and Ruiz⁵,¹¹ only by the replacement of the term $2b_3(\partial^3 y/\partial t^3)$ in Ruiz's equation by $2b_2(\partial^3 y/\partial x^2 \partial t)$. This replacement avoids a nonphysical third unstable solution (which can lead to difficulties both analytically and numerically, see Ref. 7 for more details).

To take into account the modification of the tone with the dynamic, we use a nonlinear hammer model. Ghosh⁸ was one of the first to propose a nonlinear spring model of the hammer felt,² obeying the power law:

<table><caption>TABLE I. Symbols used.</caption><thead><tr><th>Symbol</th><th>Definition</th></tr></thead><tbody><tr><td>y<sub>s</sub></td><td>String displacement</td></tr><tr><td>y<sub>x0</sub></td><td>Initial string displacement</td></tr><tr><td>u<sub>s</sub></td><td>String velocity</td></tr><tr><td>u<sub>x0</sub></td><td>Initial string velocity</td></tr><tr><td>c</td><td>Wave speed</td></tr><tr><td>κ</td><td>String stiffness coefficient</td></tr><tr><td>L</td><td>String length</td></tr><tr><td>b<sub>1</sub></td><td>String loss parameter 1</td></tr><tr><td>b<sub>2</sub></td><td>String loss parameter 2</td></tr><tr><td>yH</td><td>Hammer displacement</td></tr><tr><td>u<sub>H</sub></td><td>Hammer velocity</td></tr><tr><td>u<sub>H0</sub></td><td>Initial hammer velocity</td></tr><tr><td>F<sub>H</sub></td><td>String force acting on the hammer</td></tr><tr><td>f<sub>H</sub></td><td>Hammer force acting on the string</td></tr><tr><td>K<sub>H</sub></td><td>Felt stiffness coefficient</td></tr><tr><td>p</td><td>Stiffness exponent</td></tr><tr><td>μ</td><td>Felt loss coefficient</td></tr><tr><td>g</td><td>Window function</td></tr><tr><td>x<sub>0</sub></td><td>Contact location on the string</td></tr></tbody></table>

$$ f_H(x) = K_H x^p, \quad (2) $$

where $x$ refers to the compression of the felt, $K_H$ is the stiffness coefficient, and $p$ the stiffness exponent. Using Ghosh's model and Ruiz's string wave equation, Chaigne and Askenfelt⁶,¹² have produced realistic simulation of the hammer-string interaction. For a more realistic model, some papers suggest a hysteretic law for the hammer felt. The felt is actually compressed and extended several times during the contact and, because the relaxation is not instantaneous, the hardness of the felt increases. Boutillon⁹ modeled this phenomenon using several values of $p$ for the increasing and decreasing part of the spring stiffness characteristic. But with this model the felt deformation tends to zero with the unload- ing of the force although measurements¹³ show that the felt is still deformed after the force is removed. Stulov² proposed a mathematical model to account for the hysteretic feature of the felt, starting from a simple model of material with memory (obtained by replacing constant elastic parameters by time-dependent operators). Giordano and Millis¹⁴ have shown that this model reproduces well the force exerted by a hammer on a sensor, in the static case. They have also mounted a very small accelerometer on the string, opposite to the hammer-string contact, and shown that this model is not well adapted for simulating the force characteristic in “dynamic” conditions.

The model we used was originally developed by Hunt and Crossley¹⁰ to describe the interaction between two colliding objects (and adapted by Rochesso and Avanzini in Ref. 15 to model the hammer-string contact hysteresis). The losses are introduced by adding a first time derivative to Ghosh's model. Rochesso and Avanzini show that the force-compression diagrams are close to the ones of a real hammer. The force depends on the felt compression $x$ and of the compression speed $\partial x/\partial t$,

$$ f_H(x) = K_H x^p \left( 1 + \mu \frac{\partial x}{\partial t} \right), \quad (3) $$

where $\mu$ is the felt loss coefficient ($\mu > 0$). The system of PDEs describing the hammer-string contact is then
$$ \frac{\partial^2 y_s}{\partial t^2} = c^2 \frac{\partial^2 y_s}{\partial x^2} - \kappa^2 \frac{\partial^4 y_s}{\partial x^4} - 2b_1 \frac{\partial y_s}{\partial t} + 2b_2 \frac{\partial^3 y_s}{\partial x^2 \partial t} \\ + g(x, x_0) f_H(y_s - y_H), \qquad (4) $$

$$ M_H \frac{d^2 y_H}{dt^2} = -F_H(y_s - y_H) \qquad (5) $$

with the boundary and initial conditions:

$$ y_s(x, 0) = y_{s0}(x), \quad \frac{\partial y_s}{\partial t}(x, 0) = v_{s0}, \quad y_0(0, t) = y_s(L, t) = 0, \\ \frac{\partial^2 y_s}{\partial x^2}(0, t) = \frac{\partial^2 y_s}{\partial x^2}(L, t) = 0, \qquad (6) $$

$$ y_H(0) = y_s(x_0, 0), \quad \frac{\partial y_H}{\partial t}(0) = v_{H0} \qquad (7) $$

where

$$ F_H(u) = \begin{cases} K_H u^p \left( 1 + \mu \frac{\partial u}{\partial t} \right) & \text{if } u > 0, \\ 0 & \text{if } u \le 0 \end{cases} \qquad (8) \\ f_H(u) = F_H(u) \left\{ M_S / L \int_L g(x, x_0) dx \right\}^{-1} $$

with $y_H$ the hammer displacement, $y_{s0}$ and $v_{s0}$ are the initial string displacement and velocity, $v_{H0}$ is the initial hammer velocity, and $x_0$ is the contact location. Finally, $M_S$ (respectively, $M_H$) denote the string (respectively, hammer) mass; $L$ is the length of the string, and $g$ describes the area of the string in contact with the hammer during the interaction. [Please refer to Ref. 6 for more details about the definition of system (8).]

This nonlinear system does not admit any analytical solution. We show in Sec. III that the numerical solution of this system reproduces well the main phenomena in the hammer-string interaction, namely the nonlinear behavior with respect to the hammer velocity and the hysteresis of the hammer force.

### III. DIRECT PROBLEM: NUMERICAL SOLUTION

There is no explicit solution of the system (4)–(7). To compute the solution, we chose finite differences method in the space and time variables (one may use finite volume or finite elements methods, but in the case of one-dimensional space variables, these approaches are very similar). At first we reduce from two to one the order of the time derivative in the system (4) and (5) by introducing two new variables: the string velocity $v_s = \partial y_s / \partial t$ and the hammer velocity $v_H = \partial y_H / \partial t$. This change of variable allows one to work with the displacements and velocities of the string and the hammer. Thereby, the initial conditions may be applied exactly on each of the unknowns. Also, the order of the time derivative is reduced, simplifying the discretization.

Using these variables, one may write system (4)–(7) as a system of first-order derivatives in time:

$$ \frac{\partial y_s}{\partial t} = v_s, \qquad (9) $$

$$ \frac{\partial v_s}{\partial t} = c^2 \frac{\partial^2 y_s}{\partial x^2} - \kappa^2 \frac{\partial^4 y_s}{\partial x^4} - 2b_1 v_s + 2b_2 \frac{\partial^2 v_s}{\partial x^2} \\ + g(x, x_0) f_H(y_s - y_H, v_s - v_H), \qquad (10) $$

$$ \frac{\partial y_H}{\partial t} = v_H, \qquad (11) $$

$$ \frac{\partial v_H}{\partial t} = -\frac{1}{M_H} F_H(y_s - y_H, v_s - v_H) \qquad (12) $$

and the initial and boundary conditions are rewritten as

$$ y_s(x, 0) = y_{s0}(x), \quad v_s(x, 0) = v_{s0}, \quad y_s(0, t) = y_s(L, t) = 0, \\ \frac{\partial^2 y_s}{\partial x^2}(0, t) = \frac{\partial^2 y_s}{\partial x^2}(L, t) = 0, \qquad (13) $$

$$ y_H(0) = y_s(x_0, 0), \quad v_H(0) = v_{H0}. \qquad (14) $$

Then, $F_H$ and $f_H$ are now defined by

$$ F_H(u, v) = \begin{cases} K_H u^p (1 + \mu v) & \text{if } u > 0 \\ 0 & \text{if } u \le 0, \end{cases} \qquad f_H(u, v) = F_H(u, v) \left\{ M_S / L \int_L g(x, x_0) dx \right\}^{-1} $$

It should be noted that $F_H$ and $f_H$ are now linear with respect to the velocity variable $v$, simplifying the numerical problem.

#### A. Time discretization and initial conditions

The system above (9)–(12) may be written

$$ \frac{\partial X}{\partial t} = \mathcal{A}(X), \qquad (15) $$

where $X = (y_s, v_s, y_H, v_H)^T$ ($M$) denotes the transposed matrix of $M$ and $\mathcal{A}$ contain both the spatial differential operator and the nonlinear function corresponding to the right-hand side of Eqs. (9)–(12). Equation (15) is discretized in time by a Crank–Nicholson scheme:¹⁶

$$ \frac{X^{n+1} - X^n}{\delta t} = \frac{1}{2} (\mathcal{A}(X^{n+1}) + \mathcal{A}(X^n)). \qquad (16) $$

This scheme is used because it has good stability properties (unconditionally stable). Here, $X^n$ denotes the discrete value of $X(t)$ at the time $t=n\delta t$, where $\delta t$ is the time step.

We need to compute initial conditions on both the velocity and the displacement of the hammer and the string. When a classical formulation is used (only in displacement variable),⁶ one needs to discretize the velocity in terms of displacement and to use a pseudo-iteration time at $t=0-\delta t$ to adjust the initial velocity value. Using our formulation, one may directly apply the initial values on the velocity unknowns.
<figure><img src="image_3.png" /><figcaption>FIG. 2. Hammer force curves as a function of time, for a note C4 and three hammer velocities. The hammer parameters and the striking positions are as in Ref. 12.</figcaption></figure>

## B. Space discretization and boundary conditions

The discretization of the space derivatives in the operator A is done using a classical centered finite differences scheme. Consider the string ]0, L[ discretized by N+2 points {$x_i=ih$}$_{i=0,...,N+1}$ where $h=L/(N+1)$. Then one may approximate the space derivatives:

$$ \frac{\partial^2 X}{\partial x^2}(x_i) \approx h^{-2}(X_{i-1} - 2X_i + X_{i+1}), \quad (17) $$

$$ \frac{\partial^4 X}{\partial x^4}(x_i) \approx h^{-4}(X_{i-2} - 4X_{i-1} + 6X_i - 4X_{i+1} + X_{i+2}), \quad (18) $$

where $X_i$ denotes $X(x_i)$.

We take into account the boundary conditions coming from Eqs. (13) and (14) to apply the discretization equations (17) and (18) on the displacement at the extremities $x=0$ and $x=L$. Actually, to apply those equations at $x=h$ (respectively, $x=L-h$), the value of the displacement needs to be introduced at “virtual” nodes ($x=0, x=-h$) (respectively, $x=L, x=L+h$). Consider the case at $x=0$. The boundary condition $y_s(x=0)=0$ gives $y_{s,i=0}=0$. The boundary condition on the

string curvature $\partial^2 y_s / \partial x^2 = 0$ may be discretized at $x=0$ by $y_{s,i=-1}-2y_{s,i=0}+y_{s,i=1}=0$, which allows $y_{s,i=-1}$ to be eliminated using $y_{s,i=-1}=-y_{s,i=1}$.

## C. Numerical results

The discretization in both the time and the space domains using Eqs. (16)–(18) leads to a system that is numerically solved at each time iteration. During the hammer-string contact, this system is nonlinear: a Newton-Raphson algorithm and a Choleski matrix decomposition algorithm are then used. After the contact, the system is linear: a Choleski decomposition is performed one time for all, and each time iteration needs only two triangular system resolutions.

As shown in Figs. 2–4, the system behaves in accordance with what is foreseen by the theory. The hammer-string contact lasts about 3 ms and decreases with the increase in the initial hammer velocity. Moreover, the reflections of the waves propagating between the hammer and the agraffe can be seen on the hammer force shape (especially for 2 and 4 m s⁻¹). Figure 3 shows the force as a function of the compression of the felt. We clearly see a

<figure><img src="image_4.png" /><figcaption>FIG. 3. Hammer force curve as a function of the felt compression, for a note C4 and a hammer velocity of 4 m s⁻¹. The hammer parameters and the striking positions are $\mu=0.1$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. The spectra of velocity signal for three hammer velocities for a note C4.</figcaption></figure>

hysteretic behavior [due to the term $K_H\mu(v_s-v_H)$ of the non-linear force law] reflecting losses in the felt. In Fig. 4 the spectral modulus of the velocity signal is plotted for different initial hammer velocities $v_{H0}$. We see a modification of the spectral shape, namely increased amplitude of the high frequency partials with the hammer velocity (corresponding to increased brightness of the sound). Similar results may be found namely in Refs. 12 and 17–19.

## IV. INVERSE PROBLEM: FITTING THE PARAMETERS OF THE MODEL WITH EXPERIMENTAL DATA

We deal here with the problem of identifying the parameters of the model so that the solution is perceptually close to the measured signal. To do so, we consider the solutions $(y_s, v_s, y_H, v_H)$ of the system (9)–(14) as a function of all the physical parameters characterizing the vibration $(c, \kappa, b_1, b_2, K_H, \mu, M_H, p, v_{H0})$. We do not, however, consider the parameters $x_0$ and $M_S$ in this identification process. $x_0$ is supposed to be known accurately and since only the ratio $M_S/M_H$ plays a role in Eqs. (9)–(14), we consider only $M_H$.

Some of the parameters of the PDE model described in Sec. II can be accurately measured experimentally (the initial hammer velocity $v_{H0}$, the wave speed $c$, etc.). In contrast, parameters like the hammer stiffness $K_H$ and the felt loss coefficient $\mu$ are very difficult to measure, especially under dynamic conditions. To get a realistic and perceptively accurate estimation of those parameters, we propose to look for the set of model parameters allowing a “good” match between the synthetic and original signal. Since the dependency of the computed signal as a function of the physical parameters is not explicit, we used nonlinear optimization techniques to find this set of parameters. First, Sec. IV A describes the experimental setup. Section IV B deals with the selection of the optimization criterion accounting for the perceptual feature of the sound. Finally, Sec. IV C describes the optimization method and its numerical implementation.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Experimental setup.</figcaption></figure>

### A. Experimental setup

We have designed an experimental setup (see Fig. 5) to measure the vibration of a string struck by a hammer at different velocities. This setup has been described in Ref. 20. It consists of a massive concrete support with, on the extremes, a piece of a bridge and an agraffe taken from a real piano. The string is then tightened between the bridge and the agraffe and tuned manually. The string length is 1 m. It is clear that the string is not totally uncoupled from its support, but this experiment was satisfactory for recording signals of struck strings in order to validate our calibration method. The string is struck with a hammer linked to an electronically piloted key. By imposing different voltages on the system, the hammer velocity can be controlled in a reproducible way. The precise velocity is measured immediately after escapement through an *optic sensor* (MTI 2000, probe module 2125H) pointing at the side of the head of the hammer. The vibration on the string (velocity signal) is measured by a *laser vibrometer* (Polytec OFV 050). The signals are directly recorded on Digital Audio Tape. We collected signals corresponding to hammer velocities between 0.8 and 2.3 m s⁻¹.

### B. Choice of a perceptual criterion

The optimization criterion has to be chosen with care to take into account perceptual phenomena. For that, one has to be aware that the model is a simplified representation of reality. It does not take into account some physical phenomena leading to small perceptual contributions. This is for example the case for the beats due to the coupling between the two polarizations of the string<sup>21,22</sup> and for the finite impedance caused by the bridge.<sup>23</sup> Nevertheless, we expect the synthetic sound to be perceptually close to the original one, even though the wave forms corresponding to the signals may not be exactly identical. This latter notion implies that an error criterion based on the temporal signal behavior of the real and synthetic sounds is not appropriate for the optimization process.
On the other hand, the spectral repartition of the energy plays an important role, ensuring for example the inharmonicity of the partials and the frequency dependent damping. Since slight phase modifications in the frequency domain do not alter the perception of the sound, we based our optimization criterion on the power spectral energy density. Actually, even though the phase contains relevant information about the physical parameters, this information is also included in the spectral energy density. For example, the parameter $κ$ related to the dispersion of waves acts on the group delay (which corresponds to the derivative of the phase of the Fourier transform), but also on the frequencies of the spectrum's partials (since the wave dispersion is responsible for the inharmonicity). One may notice that, even though we are interested in the behavior of a nonstationary signal, it is not necessary to base our optimization process upon a time-frequency representation. The energy in the system is entirely determined by the initial condition, given by the hammer-string contact. After the first milliseconds, the system behaves linearly and the energy decreases in the same way until extinction. The way the energy decreases is then fully determined by the shape of the frequency components of the signal spectrum.

The optimization process consists in minimizing the quadratic error $C$ between the energy of the Fourier transform $Y_S^c$ of the string velocity $v_S$ obtained by Eqs. (4)–(7) and the energy of the Fourier transform $S^m$ of the measured signal $s$:

$$C = \frac{1}{2} \int (|Y_S^c|^2 - |S^m|^2)^2 d\omega. \quad (19)$$

Nevertheless, the direct use of such a criterion does not sufficiently take into account the hearing processes. Actually, when we first used it in the optimization process, we obtained a synthetic sound for which the attack was too smooth compared to the attack of the original sound. Moreover, the brightness of the sound was reduced. This reduction can be explained by the fact that the most important part of the energy is localized in the low frequency range of the spectrum. Consequently, this “global” criterion privileges low frequency components at the expense of high frequency ones, which are essential for the brightness and the sharpness of the attack.

To overcome this high frequency loss, we adapted the criterion to take perceptual features into account, by-following the so-called tristimulus concept.⁴ This approach consists in separating the spectral range into three spectral bands. The first band $B_1$ is centered around the fundamental frequency of the signal, whereas the second $B_2$ and the third band $B_3$ in our approach, respectively, contain the partials two to eight and all the components above the eighth partial. This approach has been successfully used for estimation of sound synthesis models by Ystad et al.²⁴ This spectral split has several advantages.

$B_1$. Tracking of the fundamental ensures good estimation of the geometric and static parameters. This fundamental component is of importance for the sound since it is highly correlated to the pitch. We used it to estimate the wave speed $c$ and the loss coefficient $b_1$.

$B_2$. Components two to eight contain a large part of the energy. These components bring most of the perceptual information related to the damping behavior of the sound, and their frequencies are responsible for the sensation of inharmonicity.²⁵ Hence, this frequency band is used to estimate the loss coefficient $b_2$ as well as the string stiffness coefficient $κ$.

$B_3$. High frequency components are also of great importance since they can dramatically change the timbre of the sound, even though their energy is low. These components are mainly correlated to the velocity and the stiffness of the hammer felt. They are characteristic of the dynamics of the playing. We used this frequency band to estimate the remaining coefficients ($K_H, \mu, M_H, p, v_{H0}$) governing the hammer behavior.

Such a decomposition of the “global” spectral power density criterion into three criteria acting on three frequency bands accounts for the limitations of the model and for the characteristics of our perception of sounds. As we discuss later, the decomposition permitted us to estimate parameters providing better synthetic sounds than the ones obtained using classical criteria.

### C. The optimization process

We describe here the optimization applied to our problem. Two kinds of optimization algorithms were used. For the bands $B_1$ and $B_2$, we used a simple gradient method since only one minimum was usually found and the criterion was regular enough to be minimized by such a method. But the optimization problem in the last band $B_3$ (concerning the hammer parameters) is more complex. Numerical experiments show the existence of local minima, and a classical gradient method does not permit convergence to the global minimum. We then considered an adapted simulated annealing method [(ASA), see Ref. 26], an algorithm known to give good results in such a situation.

For both methods, the gradient of the solution needs to be computed in regard to the physical parameters (see in the Appendix how to compute this gradient). The problem may be stated as the following. Let us denote:

$$(z_j)_{j=1,...,9} = (c, \kappa, b_1, b_2, K_H, \mu, M_H, p, v_{H0}). \quad (20)$$

A set of physical parameters $(z_j)_{j=1,...,9}$ is selected so that the system (9)–(14) gives a solution that sounds as close as possible to the experimental signal. More precisely, considering for each band the quadratic criterion $C$ as defined in Eq. (19), one denotes $\tilde{C}(z_1, ..., z_9)$ the value of $C$ computed from the solution of the problem (9)–(14) using $(z_1, ..., z_9)$ as physical parameters. Then the general optimization problem is

$$\text{Find } (z_j)_{j=1,...,9} \text{ such that } \tilde{C}(z_1, ..., z_9) \text{ is minimum.} \quad (21)$$

By using the tristimuli formalism explained in Sec. IV B, we obtain the three-step process: Consider $(z_1^0, ..., z_9^0)$ (the superscript denotes the internal loop index needed to solve the nonlinear equation) given by the measured parameters on the experimental setup.
(1) Find $(z_1^1, z_3^1)$ such that $\tilde{C}(z_1^1, z_2^1, z_3^1, z_4^0, z_5^0, z_6^0, z_7^0, z_8^0, z_9^0)$ is minimum on the band $B_1$.

(2) Find $(z_1^2, z_4^1)$ such that $\tilde{C}(z_1^1, z_2^1, z_3^1, z_4^1, z_5^0, z_6^0, z_7^0, z_8^0, z_9^0)$ is minimum on the band $B_2$.

(3) Find $(z_5^1, z_6^1, z_7^1, z_8^1, z_9^1)$ such that $\tilde{C}(z_1^1, \dots, z_9^1)$ is minimum on the band $B_3$.

To solve each of those problems, we then used the two optimization algorithms described below.

## 1. Gradient algorithm

Concerning the minimization of the criterion in the bands $B_1$ and $B_2$, the situation is regular enough (unique minimum) to use a classical gradient algorithm. Starting from the expression of the criterion (19), one can derive with respect to $(z_j^0)_{j=1,\dots,9}$:

$$ \nabla_z \tilde{C} = 2 \int (|Y_S^c|^2 - |S^m|^2) |Y_S^c| \nabla_z Y_S^c d\omega. $$

$\nabla_z$ denotes the gradient with respect to the physical parameters. Because the Fourier transform is linear, one can deduce $\nabla_z Y_S^c$ from $\nabla_z v_S$ (see the Appendix for how to compute the gradient $\nabla_z v_S$). To converge to the local minimum we then constructed the following algorithm:

* $\epsilon$ given; $(z)_{i=1,n}=(z^0)_{i=1,n}$ given

* while $\tilde{C} > \epsilon$ do

Compute $(y_s, v_s, y_H, v_H)$ solution of (9)–(14) with parameter set $(z)_{i=1,n}$

* for $j=1,n$ do

Compute $(y_s^{(j)}, v_s^{(j)}, y_H^{(j)}, v_H^{(j)})$ as explained in the Appendix

Compute $\tilde{C}^{(j)}$ using $(y_s^{(j)}, v_s^{(j)}, y_H^{(j)}, v_H^{(j)})$ ($\tilde{C}^{(j)}, j^{th}$ component of $\nabla_z \tilde{C}$)

* end for

chose $h$ such that $\tilde{C} + \nabla_c \tilde{C} \cdot h \le \tilde{C}$

* for $j=1,n$ do

$z_j \leftarrow z_j + h_j$

* end for

* end while

This algorithm was successfully used to find the minimum of the criterion $C$ in bands $B_1$ and $B_2$.

## 2. Simulated annealing algorithm

For the last frequency band $B_3$, one needs an algorithm that is able to escape local minima in order to find the global minimum. Let us then consider an optimization method more adapted to this situation ASA. We give here only the basic principles of this method. A detailed description can be found in Ref. 26. As its name implies, the simulated annealing exploits an analogy between the way in which a metal cools and freezes into a minimum energy crystalline structure (the annealing process) and the search for a minimum in a more general system. As cooling proceeds, the system becomes more ordered and approaches a “frozen” ground state at “Temperature” $T=0$. Hence the process can be thought of as an adiabatic approach to the lowest energy state. If the initial temperature of the system is too low or cooling is done

<table><caption>TABLE II. Estimated and optimized physical coefficients for the string.</caption><thead><tr><th rowspan="2"></th><th colspan="4">String parameters</th></tr><tr><th>c</th><th>κ</th><th>b<sub>1</sub></th><th>b<sub>2</sub></th></tr></thead><tbody><tr><td>Initial</td><td>310.6</td><td></td><td>0.500</td><td>1.004E-10</td></tr><tr><td>Final</td><td>304.101</td><td>1.737</td><td>0.568</td><td>1.628E-06</td></tr></tbody></table>

insufficiently slowly the system may become quenched forming defects or freezing out in metastable states (i.e., trapped in a local minimum energy state). Applied to our situation, this algorithm may be summarized as the following: considering an initial point $(z^0)_{i=1,n}$, a new point $(z)_{i=1,n}$ is generated using a given probability distribution function $g(z)$ and the criterion $\tilde{C}$, evaluated at this point. The new point is accepted or refused, following an acceptance function $h(z^0, z)$ that depends both on the parameter called Temperature, and on the difference between values of the criterion $\tilde{C}$ on the two points $(z^0, z)$. If the new point is accepted, it becomes the old one, the Temperature is decreased, and the process is iterated.

## V. RESULTS

We applied the optimization process on the sounds recorded by the experimental setup. As initial parameters, measured and published values were used. For the string, the parameters were extracted from measurements on the experimental setup. The parameters corresponding to the hammer characterization were collected in the literature.¹² Our results showed that the optimization process leads to a new set of physical parameters. The string parameters obtained are given in Table II. Table III shows the hammer parameters. This set of “optimized” values is the result of the optimization process and depends on the choice of the convergence criterion. The three-band criterion we used made it possible to estimate new parameters ensuring a better synthesis from a perceptual point of view. It is not sure, however, that they correspond to the actual physical values. Nevertheless, as shown in Tables II and III, these values are physically consistent and do not deviate too much from the starting point. Moreover, those optimized parameters were very similar from one measurement to another (for the same hammer velocity). The wave speed $c$ was reduced, meaning that the string tension was weaker than expected. The loss coefficients, $b_1$ and $b_2$ were greater than those estimated on the experimental setup. For the hammer, note that the felt loss coefficient increased, as did the initial hammer velocity. For the stiffness exponent $p$, the value was almost unchanged.

This set of new parameters was then used to compute a synthetic signal [using Eqs. (4)–(7)]. In Fig. 6 are plotted the spectral modulus of the original and the synthesized signals

<table><caption>TABLE III. Estimated and optimized physical coefficients for the hammer.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Hammer parameters</th></tr><tr><th>K<sub>H</sub></th><th>μ</th><th>p</th><th>M<sub>H</sub></th><th>ν<sub>0</sub></th></tr></thead><tbody><tr><td>Initial</td><td>1.516 47E+09</td><td>0.1000</td><td>2.55</td><td>0.006 790</td><td>0.740 78</td></tr><tr><td>Final</td><td>1.099 49E+10</td><td>0.1998</td><td>2.56</td><td>0.005 687</td><td>0.799 93</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>FIG. 6. Spectrum modulus (velocity signal) for a hammer velocity of 0.84 m s⁻¹. Top: original measured signal, bottom: synthesized optimized signal.</figcaption></figure>

with optimized parameters. Even though they are not identi-
cal, they look very similar. The amplitudes of the partials are
in good agreement for the first seven components, but differ
for higher frequencies. Figure 7 represents the wave shape of
the real and the synthesized sounds, before and after the
optimization. The similarity between the signals has been
clearly improved by the optimization of the parameters. Nev-
ertheless, the synthesized signal contains slightly more en-

ergy in the high frequency range. As discussed in Sec. IV B,
the tristimuli criterion, rather than the “global” criterion, was
chosen for its better ability to take into account the energy of
high frequency components. We believe here that the main
part of the differences is due not only to the optimization
method but also to the model approximation. Obviously, the
optimization method is constrained by the model limitation;

<figure><img src="image_8.png" /><figcaption>FIG. 7. Comparison of the measured and computed velocities of the string, for a hammer velocity of 1.1 m s⁻¹. Dashed line: measured signal, dash-dotted line: synthesis before optimization, plain synthesis after optimization.</figcaption></figure>
better results would probably be obtained if first the model itself were reconsidered. This aspect is not in the scope of this article.

Another important point validating our approach is that with the optimized parameters the model still behaves as the measurements with the dynamic. This means that, if we compute the model for another hammer velocity, the signal obtained is still close to the signal measured for the same hammer velocity.

Informal listening test has shown that the quality of the resulting sounds increases using the optimized parameters (sound examples can be found at Ref. 27). The attack is more realistic and the decay closer to the original even though the tone still lacks realism. We believe that the explanation for this lack resides in the simplicity of the hammer model, which is not able to accurately simulate the nonlinear behavior of the hammer-string interaction. Actually, previous synthesis works have proved the accuracy of the string model itself.⁷

## VI. CONCLUSION

One of the main issues when using physical modeling for synthesis purposes is the fitting of the model's parameters. These parameters can sometimes be extracted from data obtained by experiments, but in some cases the dynamic behavior of the models reduces the accuracy of static measurements. This is the case for piano sound models, in which it is difficult to estimate parameters like the stiffness coefficient and exponent characterizing the hammer.

We have addressed here the problem of parameter fitting for piano modeling, aiming at estimating the optimal parameters for synthesis purposes. This means that rather than focusing on the precision of the physical parameters, we wanted the synthetic sound to be as close as possible to the real one. This led us to a technique associating optimization processes and perceptual concepts. To take into account our perception of sounds, a minimization criterion (three-band criterion) was designed following the tristimulus concept. The frequency range was separated into three bands, respectively, containing the fundamental, partial 2 to 8, and the remaining high frequency components. Having in mind that real sounds are much more complex than the ones generated by the physical model, we compared signal characteristics that are as little as possible sensitive to small variations. This eliminated an estimation based on the time signal to the benefit of its spectral energy density. This choice does not minimize the influence of the phase information, and consequently the time behavior of the signal, since most of it is also contained in the frequency distribution of the energy. We showed how each band can be used to estimate the model's parameters thanks to an optimization process that combines a gradient method and a simulated annealing algorithm aimed at avoiding convergence problems in case of multiple local minima. The model's parameters make it possible to generate piano tones that are perceptually better than those obtained using a “global” criterion. Fine piano characteristics that are generally poorly synthesized are thus reproduced. This is the case for the attack time and the brightness. Nevertheless,

slight differences between the real and the synthetic sounds can still be heard. We believe these differences are mainly due to the weakness of the model, especially for the string-hammer interaction. Also note that even though the parameters obtained by our technique are optimized for synthesis purposes, their values correspond to what is expected from a physical point of view. The generality of the method proposed in this paper makes it easy to be applied to other models of music instruments and sound generators.

## APPENDIX: COMPUTATION OF THE PHYSICAL PARAMETERS DERIVATIVES

Here we explain how to get the gradient of the solution of system [Eqs. (9)–(14)] with respect to the physical parameters ($c, \kappa, b_1, b_2, K_H, \mu, M_H, p, v_{H0}$). We do not mention $v_{s0}$, $v_{s0}$, $x_0$, since we will not calculate the derivatives for those parameters. We suppose that $v_{s0}=0$, $v_{s0}=0$ and that $x_0$ can be precisely measured on the experimental setup. We will consider the first derivative with respect to one of those parameters: the celerity $c$. One may easily generalize this calculation for all the other parameters and for higher order derivatives.

We first write that $(y_s(c), v_s(c), y_H(c), v_H(c))$ are solutions of the system (9)–(14) for a value $c$ of the celerity and that $(y_s(c+h), v_s(c+h), y_H(c+h), v_H(c+h))$ are solutions of the system (9)–(14) for a value $c+h$ of the celerity. By calculating term by term the difference between the two corresponding systems, dividing by $h$, and passing to the limit when $h$ tends to zero, we obtain a linear system of equations, governing the derivatives $(y_s^{(1)}, v_s^{(1)}, y_H^{(1)}, v_H^{(1)})$ of $(y_s, v_s, y_H, v_H)$ in regard of $c$:

$$ \frac{\partial y_s^{(1)}}{\partial t} = v_s^{(1)}, \qquad (A1) $$

$$ \begin{aligned} \frac{\partial v_s^{(1)}}{\partial t} = & c^2 \frac{\partial^2 y_s^{(1)}}{\partial x^2} + 2c \frac{\partial^2 y_s}{\partial x^2} - \kappa^2 \frac{\partial^4 y_s^{(1)}}{\partial x^4} \\ & - 2b_1 v_s^{(1)} + 2b_2 \frac{\partial^2 v_s^{(1)}}{\partial x^2} + g(x, x_0) \\ & \times \nabla f_H(y_s - y_H, v_s - v_H) \cdot {}^t(y_s^{(1)}, v_s^{(1)}, y_H^{(1)}, v_H^{(1)}), \end{aligned} \qquad (A2) $$

$$ \frac{\partial y_H^{(1)}}{\partial t} = v_H^{(1)}, \qquad (A3) $$

$$ \frac{dv_H^{(1)}}{dt} = - \frac{1}{M_H} \nabla F_H(y_s - y_H, v_s - v_H) \cdot {}^t(y_s^{(1)}, v_s^{(1)}, y_H^{(1)}, v_H^{(1)}) \qquad (A4) $$

and the initial and boundary conditions are given by

$$ \begin{aligned} y_s^{(1)}(x, 0) = 0, \quad v_s^{(1)}(x, 0) = 0, \quad y_s^{(1)}(0, t) = y_s^{(1)}(L, t) = 0, \\ \frac{\partial^2 y_s^{(1)}}{\partial x^2}(0, t) = \frac{\partial^2 y_s^{(1)}}{\partial x^2}(L, t) = 0, \end{aligned} \qquad (A5) $$
$$y_H^{(1)}(0) = y_s^{(1)}(x_0, 0), \quad v_H^{(1)}(0) = 0, \qquad (A6)$$

where $\nabla f(y_s - y_H, v_s - v_H)$ denote the gradient of $f(y_s - y_H, v_s - v_H)$ in regard to $(y_s, v_s, y_H, v_H)$ computed at the point $(y_s, v_s, y_H, v_H)$. The derivatives $(y_s^{(j)}, v_s^{(j)}, y_H^{(j)}, v_H^{(j)})_{j=1,...,9}$ in regard to the physical parameters are computed following the same scheme (16)–(18) as for the string and hammer displacements and velocities. More precisely, at each time iteration, once the values $(y_s^n, v_s^n, y_H^n, v_H^n)$ computed, the values $(y_s^{(j,n)}, v_s^{(j,n)}, y_H^{(j,n)}, v_H^{(j,n)})$ are obtained just by the solution of one linear system,

$$\frac{\partial y_s^{(j)}}{\partial t} = v_s^{(j)}, \qquad (A7)$$

$$\begin{aligned} \frac{\partial v_s^{(j)}}{\partial t} = & c^2 \frac{\partial^2 y_s^{(j)}}{\partial x^2} - \kappa^2 \frac{\partial^4 y_s^{(j)}}{\partial x^4} - 2b_1 v_s^{(j)} + 2b_2 \frac{\partial^2 v_s^{(j)}}{\partial x^2} \\ & + (c^2)^{(j)} \frac{\partial^2 y_s^{(j)}}{\partial x^2} - (\kappa^2)^{(j)} \frac{\partial^4 y_s^{(j)}}{\partial x^4} \\ & - 2(b_1)^{(j)} v_s + 2(b^2)^{(j)} \frac{\partial^2 v_s^{(j)}}{\partial x^2} + g(x, x_0) \\ & \times \nabla f_H(y_s - y_H, v_s - v_H) \cdot {}^t(y_s^{(j)}, v_s^{(j)}, y_H^{(j)}, v_H^{(j)}) \\ & + g(x, x_0) f_H^{(j)}(y_s - y_H, v_s - v_H), \end{aligned} \qquad (A8)$$

$$\frac{\partial y_H^{(j)}}{\partial t} = v_H^{(j)}, \qquad (A9)$$

$$\begin{aligned} \frac{\partial v_H^{(j)}}{\partial t} = & -\frac{1}{M_H} \nabla F_H(y_s - y_H, v_s - v_H) \cdot {}^t(y_s^{(j)}, v_s^{(j)}, y_H^{(j)}, v_H^{(j)}) \\ & - \left(\frac{1}{M_H}\right)^{(j)} F_H(y_s - y_H, v_s - v_H) \\ & - \frac{1}{M_H} F_H^{(j)}(y_s - y_H, v_s - v_H) \end{aligned} \qquad (A10)$$

and the initial and boundary conditions are given by

$$y_s^{(j)}(x, 0) = y_{s0}^{(j)}, \quad v_s^{(j)}(x, 0) = v_{s0}^{(j)}, \qquad (A11)$$

$$y_s^{(j)}(0, t) = y_s^{(j)}(L, t) = 0, \quad \frac{\partial^2 y_s^{(j)}}{\partial x^2}(0, t) = \frac{\partial^2 y_s^{(j)}}{\partial x^2}(L, t) = 0, \qquad (A12)$$

$$y_H^{(j)}(0) = y_s^{(j)}(x_0, 0), \quad v_H^{(j)}(0) = v_{H0}^{(j)}. \qquad (A13)$$

<footnote n="1">¹J. Bensa, D. Dubois, R. Kronland-Martinet, and S. Ystad, in *Perceptive and Cognitive Evaluation of a Piano Synthesis Model*, LNCS 3310 of CMMR 2004, Lecture Note in Computer Science Series Vol. 1, edited by U. Will (Springer, Berlin, 2004), pp. 232–245.</footnote>

<footnote n="2">²A. Stulov, “Hysteretic model of the grand piano hammer felt,” J. Acoust. Soc. Am. 97, 2577–2585 (1995).</footnote>

<footnote n="3">³J. Bensa, A. Gipouloux, and R. Kronland-Martinet, “Understanding the influence of physical parameters on piano sound by the use of an optimization method,” in Proceedings of the International Symposium of Musical Acoustics (ISMA2001), The Musical and Architectural Acoustics Laboratory Perugia, Italy, 2001.</footnote>

<footnote n="4">⁴H. Pollard and E. Jansson, “A tristimulus method for the specification of musical timbre,” Acustica 51, 162–171 (1982).</footnote>

<footnote n="5">⁵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="6">⁶A. Chaigne and A. Askenfelt, “Numerical simulation of piano strings. I. a physical model for a struck string using finite difference methods,” J. Acoust. Soc. Am. 95, 1112–1118 (1994).</footnote>

<footnote n="7">⁷J. Bensa, S. Bilbao, R. Kronland-Martinet, and J. O. Smith III, “The simulation of piano string vibration: From physical model to finite difference schemes and digital waveguides,” J. Acoust. Soc. Am. 114, 1095–1107 (2003).</footnote>

<footnote n="8">⁸M. Gosh, Indian J. Phys. 7, 365–382 (1932).</footnote>

<footnote n="9">⁹X. Boutillon, “Model for piano hammers: Experimental determination and digital simulation,” J. Acoust. Soc. Am. 83, 746–754 (1988).</footnote>

<footnote n="10">¹⁰K. Hunt and F. Crossley, “Coefficient of restitution interpreted as damping in vibroimpact,” ASME J. Appl. Mech. 42, 440–445 (1975).</footnote>

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

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

<footnote n="13">¹³T. Yanagisana and K. Nakanura, “Dynamic compression characteristics of piano felt,” J. Acoust. Soc. Jpn. 40, 725–729 (1984).</footnote>

<footnote n="14">¹⁴N. Giordano and J. P. Millis, “Hysteretic behavior of piano hammers,” in Proceedings of the International Symposium of Musical Acoustics, The Musical and Architectural Acoustics Laboratory, Perugia, Italy, 2001, pp. 237–40.</footnote>

<footnote n="15">¹⁵F. Avanzini and D. Rochesso, “Non linear contact force,” in Proceedings of the COST G6 Conference on Digital Audio Effects, Department of Computer Science and Information Systems, University of Limerick, Limerick, Ireland, 2001.</footnote>

<footnote n="16">¹⁶J. C. Strikwerda, *Finite Difference Schemes and Partial Differential Equations* (Wadsworth and Brooks, Pacific Grove, CA, 1989).</footnote>

<footnote n="17">¹⁷D. E. Hall, “Piano string excitation. VI. Nonlinear modeling,” J. Acoust. Soc. Am. 92, 95–105 (1992).</footnote>

<footnote n="18">¹⁸D. E. Hall and A. Askenfelt, “Piano string excitation. V. Spectra for real hammers and strings,” J. Acoust. Soc. Am. 83, 1627–1638 (1988).</footnote>

<footnote n="19">¹⁹F. Avanzini, B. Bank, G. Borin, G. De Poli, and D. Rocchesso, “Musical instrument modeling: The case of the piano,” in Proceedings of the Workshop on Current Research Directions in Computer Music, MOSART Research training network, Barcelona, 2001, pp. 124–125.</footnote>

<footnote n="20">²⁰J. Bensa, K. Jensen, and R. Kronland-Martinet, “A hybrid resynthesis model for hammer-strings interaction of piano tones,” EURASIP Journal on Applied Signal Processing 2004, 1021–1035 (2004).</footnote>

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

<footnote n="22">²²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="23">²³N. Giordano, “Mechanical impedance of a piano soundboard,” J. Acoust. Soc. Am. 103, 2128–2133 (1998).</footnote>

<footnote n="24">²⁴S. Ystad, “Sound modeling applied to flute sounds,” J. Audio Eng. Soc. 48, 810–825 (2000).</footnote>

<footnote n="25">²⁵H. Järvelainen, V. Välimäki, and M. Karjalainen, “Audibility of inharmonicity in strings instruments sounds, and implications to digital sound synthesis,” in Proceedings of the International Computer Music Conference, Computer Music Association, Beijing, 1999, pp. 359–362.</footnote>

<footnote n="26">²⁶L. Ingber, Adaptive Simulated Annealing, http://www.ingber.com/#ASA.</footnote>

<footnote n="27">²⁷http://www.lma.cnrs-mrs.fr/~kronland/JASA/sounds.html.</footnote>

<footnote n="28">²⁸*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 No. 64. Available online at http://www.speech.kth.se/music/5_lectures/.</footnote>