Parameter fitting for piano sound synthesis by physical modeling

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

Parameter fitting for piano sound synthesis by physical modeling Julien Bensa, a͒ Olivier Gipouloux, b͒ 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. 1 Among numerous perceptually important contributions to piano tones, one can cite for example the inharmonicity of the partials, the frequencydependent 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. 2 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 a͒ 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. 4 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 5-7 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. 2, [8][9][10] The model of wave propagation on the string we used was proposed in Ref. 7 and is given by ‫ץ‬ 2 y s ‫ץ‬t 2 = c 2 ‫ץ‬ 2 y s ‫ץ‬x 2 − 2 ‫ץ‬ 4 y s ‫ץ‬x 4 − 2b 1 with y s the string displacement, c the wave speed, 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 timespace derivative terms. This model differs from that of Hiller and Ruiz 5,11 only by the replacement of the term 2b 3 ‫ץ͑‬ 3 y / ‫ץ‬t 3 ͒ in Ruiz's equation by 2b 2 ‫ץ͑‬ 3 y / ‫ץ‬x 2 ‫ץ‬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 8 was one of the first to propose a nonlinear spring model of the hammer felt, 2 obeying the power law: 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 6,12 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 9 modeled this phenomena 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 unloading of the force although measurements 13 show that the felt is still deformed after the force is removed. Stulov 2 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 14 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 10 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 forcecompression diagrams are close to the ones of a real hammer. The force depends on the felt compression x and of the compression speed ‫ץ‬x / ‫ץ‬t, where is the felt loss coefficient ͑ Ͼ 0͒. The system of PDEs describing the hammer-string contact is then  Contact location on the string ‫ץ‬ 2 y s ‫ץ‬t 2 = c 2 ‫ץ‬ 2 y s ‫ץ‬x 2 − 2 ‫ץ‬ 4 y s ‫ץ‬x 4 − 2b 1 ‫ץ‬y s ‫ץ‬t + 2b 2 ‫ץ‬ 3 y s ‫ץ‬x 2 ‫ץ‬t with the boundary and initial conditions: where 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 hammerstring 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 = ‫ץ‬y s / ‫ץ‬t and the hammer velocity v H = ‫ץ‬y H / ‫ץ‬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: and the initial and boundary conditions are rewritten as Then, F H and f H are now defined by 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 where X = t ͑y s , v s , y H , v H ͒ ͓ t ͑M͒ denotes the transposed matrix of M͔ and 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: 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␦t, where ␦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͒, 6 one needs to discretize the velocity in terms of displacement and to use a pseudo-iteration time at t =0−␦t to adjust the initial velocity value. Using our formulation, one may directly apply the initial values on the velocity unknowns.

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: 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 ‫ץ‬ 2 y s / ‫ץ‬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 hammerstring 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 −1 ͒. Figure 3 shows the force as a function of the compression of the felt. We clearly see a hysteretic behavior ͓due to the term K H ͑v s − v H ͒ of the nonlinear 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 , , b 1 , b 2 , K H , , 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 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.

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 extremities, 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 −1 .

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 21,22 and for the finite impedance caused by the bridge. 23 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:

͑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, byfollowing the so-called tristimulus concept. 4 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. 24 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. 25 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 , , 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: 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 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 Find ͑z j ͒ j=1,. . .,9 such that C͑z 1 , ... ,z 9 ͒ is minimum.

͑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.
To solve each of those problems, we then used the two optimization algorithms described below.

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,. . ., 9 : ٌ z denotes the gradient with respect to the physical parameters. Because the Fourier transform is linear, one can deduce ٌ z Y S c from ٌ z v s ͑see the Appendix for how to compute the gradient ٌ z v s ͒. To converge to the local minimum we then constructed the following algorithm: This algorithm was successfully used to find the minimum of the criterion C in bands B 1 and B 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 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 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 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. 12 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 ve-locity͒. 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 with optimized parameters. Even though they are not identical, 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. Nevertheless, 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; 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. 7

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 stringhammer 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 , , b 1 , b 2 , K H , , M H , p , v H0 ͒. We do not mention y s0 , v s0 , x 0 , since we will not calculate the derivatives for those parameters. We suppose that y 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