Physical Modeling of Membranes for Percussion Instruments. Federico Fontana, Davide Rocchesso. Centro di Sonologia Computazionale, Dipartimento di ...
ACUSTICA·
acta acustica
Vol. 84 (1998) 529-542
©
S. Hirzel Verlag· EAA
529
Physical Modeling of Membranes for Percussion Instruments
Federico Fontana, Davide Rocchesso Centro di Sonologia Computazionale, Dipartimento di Elettronica e Informatica, Universita degli Studi di Padova, Via Gradenigo 6/A, 35131 Padova, Italy
Summary Recent research on Physical Modeling has led to 2-D discrete-time structures based on the Digital Waveguides. These structures are well suited for efficient yet accurate simulation of wave propagation in an ideal membrane. Nevertheless, real membranes exhibit a different behaviour, due to the environmental conditions and to the material they are made of. In this work we consider some aspects, crucial for the audio signal, of the physical phenomena concerning real membranes, and we will develop a 2-D waveguide model encompassing the effects of these aspects. In order to excite the simulated membrane, we will consider a hammer model previously developed for piano strings, and here adapted to fit the hammer-membrane interaction. PACS no. 43.75.Hi, 43.75.Wx, 43.40.Dx
1. Introduction Physical Modeling is becoming one of the main methods for sound synthesis, due to the naturalness of sound dynamics under different playing conditions. Remarkable results have been obtained especially in the synthesis of stringed and wind-instrument sounds, where wave propagation is essentially one-dimensional [I]. Efficient techniques, such as Digital Waveguide Networks [2], have been devised to provide an efficient yet accurate model of resonators (e.g., the string of the violin or the bore of the clarinet) and their interactions with exciting elements (e.g., the bow or the reed). The simulation of two- and three-dimensional wave propagation media poses important problems due to the fact that any waveform can be decomposed into an infinite sum of plane waves, and any feasible model has to reduce this sum to a small number of terms. The digital waveguide paradigm has been adapted to the multidimensional case by introduction of waveguide meshes, which can be thought of as approximations of continuous surfaces or volumes by grids of wires [3,4]. In this paper we focus on two-dimensional elastic media, namely circular membranes, which are often found as resonators in percussion instruments. Prior work on this subject includes the rectangular and hexagonal waveguide meshes as efficient structures for membrane modeling [3,5,4]. We propose to use a triangular tessellation of the surface since this mesh structure shows better properties in minimizing spatial aliasing and numerical dispersion error1. The resulting structure is called the Triangular Digital Waveguide Mesh. Several improvements can be made to the basic model in order to achieve realistic sounds. For example, nonlin-
ear filters can be inserted at the rim to boost high frequencies and obtain sounds similar to cymbals [7]. Even simple membrane-based instruments show a behavior remarkably different from the ideal elastic membrane, due to several reasons. One important fact is that the air load strongly affects the position of resonances, and in many practical cases it lines up some of the partials along a harmonic series [8], thus giving pitch to the sound. By using loaded waveguide junctions we can take this effect into account, and insert attenuation due to internal losses and radiation, so that the alignment of partials and the decay time can be tuned based on the sound requirements. Another important component for sound quality is a correct excitation. To this end, the results of studies in hammer-string interaction, as found in the piano literature [9, 10, II, 12, 13, 7], can be adapted to the membrane. However, most of the models of hammer developed so far suffer from instability problems when run at the sampling rate requested by real-time applications. Moreover, it is not always straightforward how to tune the hammer parameters in order to get the right hardness and weight. Indeed, a versatile hammer model has been recently developed [14,9], and we succeeded in the task of coupling it with the triangular waveguide mesh. As a result of our investigations and experiments, we propose a complete system which is a highly parallel structure •• made of simple elements (second order digital filters) connected by delay units. This suggests best implementation of it using two-dimensional pipelines within custom chips. The project has been developed within the context of Physical Modeling, even though some complex phenomena of real membranes, such as the air load, are taken into account by simulating their effects rather than their causes.
Received 13 August 1996, accepted 7 July 1997. 1 While revising this article we have been informed of new interesting developments of the rectangular waveguide mesh [6]. Namely, the simulation is restated as a finite-difference scheme and an interpolation technique is used for reducing the numerical artifacts due to discretization.
2. The Triangular Digital Waveguide Mesh In this section, a discrete-time structure modeling an ideal membrane is given. First, we summarize the background
ACUSTICA
530
Fontana and Rocchesso: Physical modeling of membranes
theory and the tools needed to define the model. Then we will show that the same equations and structures can be obtained by means of a Finite Difference Scheme.
. acta acustica Vol. 84 (1998)
For a lossless (series) junction of N ideal strings of impedance respectively ZO,l, ... ,ZO,N, we can write at time t (due to losslessness) the two equations
2.1. One and Two-dimensional Wave Propagation The one-dimensional wave equation for the transverse velocity v (t, x) at time t and longitudinal position x of an ideal vibrating string is (1)
where Vi and Ii are respectively the transverse velocity and the force at the junction of the i-th string, and Vj is the velocity of the junction. Using equations (2), (4), (5) and (6) we find out the relation [16]
Vi-(t) where c is the traveling speed of the wave. Solving equation (1) without any particular condition we obtain the D' Alembert solution [15]
V(t, x)
= v+(x - ct) + v_(x + ct),
N
=
L ZO,kVk+(t) 2 k=l N - Vi+(t),
(7)
LZO,k k=l
(2)
meaning that the transverse velocity is the superposition of two velocity waves, v+ and v_, traveling into opposite directions. The two-dimensional wave equation for the transverse velocity v (t, x, y) at time t and position (x, y) of an ideal vibrating membrane is
= Vj(t) - Vi+(t)
which gives the wave signal Vi-, outgoing to the i-th string, in terms of the incoming waves vi+' If every string has the same impedance, equation (7) simplifies in (8)
(3) where c is, again, the traveling speed of the wave. Solving equation (3) without any particular condition we obtain [15]
(27r
v(t,x,y)
=
J
o
ve(xcosO+ysinO-ct)dO,
meaning that the transverse velocity is the superposition of an infinite number of velocity plane waves, Ve , traveling in all directions. The same equations apply to the force I. In particular, for the ideal string we have
I(t, x)
=
f+(x - ct)
+ I-(x + ct),
(4)
where 1+ and I-are force waves. The wave speed for the string is c = JT /1], where T is the tension of the string and 1] its linear density. The wave speed for the membrane is c = JTm/a, where Tm is the tension per unit length of the membrane and a its superficial density.
2.2. Junction of Ideal Strings The impedance Zo = v'T7i for an ideal string gives a relation between velocity waves and force waves:
1+
=
f- =
Zov+, -Zov_.
(5)
2.3. The Digital Waveguide and The Lossless Scattering Junction As models for the ideal string and for the junction of ideal strings we use respectively the Digital Waveguide and the Lossless Scattering Junction. For our purpose, we just give a simple definition of them, suggesting [2, 16] for a thorough treatment. To define the Digital Waveguide we rewrite equation (2) in the discrete-time, with time sample T, obtaining the equation
v(nT, x) = v+(x - cnT)
+ v_ (x + cnT),
where x is now defined in the discrete set of positions X = {xo ± cnT , n = 0, 1, ... }, Xo being an arbitrary position. If the velocity wave is known for every element in X at a time nT, velocity at the following time sample for a fixed position x can be found out summing the wave samples located in the adjacent positions, x ± cT :
v(nT
+ T,
x) = v+(x - c(nT + T)) + v_(x + c(nT = v+((x - cT) - cnT) + v_((x + cT)
+ T)) + cnT),
again for every element x EX. This suggests considering a structure composed of a pair of delay lines: in each one of them a discrete-time signal is present, and they travel independently in opposite directions (see Figure 1). The Digital Waveguide does not exhibit aliasing error for signals of bandwidth less than Fs /2, where is Fs = lJT.
ACUSTICA·
acta acustica
Fontanaand Rocchesso:Physicalmodelingof membranes
Vol. 84 (1998)
531
2.4. Digital Waveguide Meshes
Figure I. Digital Waveguidesconnectedthrougha Lossless Scattering Junction.
-
•
~
•
e that the mesh implements a second-order difference scheme for the differential wave equation (3) of the ideal membrane;
•
e that signal components travel along the mesh at a speed depending on their frequency or, in other words, that there is a dispersion error on the mesh. Moreover, this dispersion error is dependent on direction, and therefore the simulated wave propagation is anisotropic.
-
--
-
-
-
-
-
More recently, Van Duyne and Smith have studied different topologies, such as the hexagonal mesh for 2-D membranes and the Tetrahedral Digital Waveguide Mesh for 3-D spaces [4]. After a careful examination ofthe features of the previously-analyzed topologies, we decided to investigate the properties of the Triangular Digital Waveguide Mesh, in the hope that it can provide better accuracy than the rectangular mesh at a reasonable complexity. As a partial justification of such a choice we can say that, from a signal-processing viewpoint, an hexagonal sampling (which corresponds to a triangular tiling) of a circularly-bandlimited 2-D signal minimizes the aliasing artifacts for a given sampling density [17]. This consideration is applicable to (circularly bandlimited) static pictures of waves generated by a lumped excitation of an isotropic membrane. Moreover, as we will show, even the dynamic behavior of the triangular mesh is much less sensitive to direction. The triangular mesh is shown in Figure 3, where the main coordinate axes, x and y, are indicated together with two other axes, 1 and m, which follow the natural alignment of the waveguide branches.
••
0
4
Figure2. Geometryof the 2-D Digital WaveguideMesh.
y,t
x
Figure 3. Geometryof the TriangularDigital WaveguideMesh.
To define the Lossless Scattering Junction we just rewrite equation (8) in the discrete time 2
vi_(nT)
= N
N
L vk+(nT)
We now consider 2-D structures made of Digital Waveguides of one delay unit and Lossless Scattering Junctions. Van Duyne and Smith have found how to relate these structures to the ideal membrane, and how to measure their performance. In particular, they have considered the rectangular 2-D Digital Waveguide Mesh [3] [5] shown in Figure 2, where we indicate Digital Waveguides with 'I' and Lossless Scattering Junctions with 'e', and they have shown
- vi+(nT),
(9)
A simulation of the triangular mesh structure gives Figure 4, where we can see respectively the transverse velocity signal 6 and 9 time samples after the excitation of the central junction (with a hammer model described below), and Figure 5, where we can see the interference of wavefronts coming from the excitation of two different junctions2 • Following the theoretical approach outlined above, we are now going to prove that the triangular mesh implements a second-order finite difference scheme for the differential wave equation of the ideal membrane. In this framework, f.'i\\.\.~-d\.\\~'i:~'i\C~'i,C\\~'i:'i\~'i, ~\ ~\ m:~ \l'i,~\\l\ \)~ca\l'i,~ a'i,\.a'i\dm:d
k=l
and notice that, through this relation, we can find the outgoing discrete-time velocity waves from the incoming velocity waves for a connection of N Digital Waveguides. Figure I shows two Digital Waveguides connected by a Lossless Scattering Junction.
2 For drawing the pictures, we considered a stretched version of the mesh as part of a grid made of square cells; then we filled the cells not occupied by junctions with values obtained by linear interpolation,and finally we stretchedback the grid to recover the correct size of the mesh
ACUSTlCA·
532
Fontanaand Rocchesso:Physicalmodelingof membranes
acta acustica Vol. 84 (1998)
The inversion of this system gives the two equations
av { av ay
av
av
-=-+ax 8l am
av
1
=
J3 8l
-
1
av '
J3 am
which give a rule to calculate the partial derivatives along the main axes as functions of the partial derivatives along the axes I and m. Hence, deriving again and summing the new equations on the second partial derivatives, each one with a proper weighting factor, we have
Finally, using equation (3), we obtain (10)
Figure 4. Evolutionof the wavefronton the triangularmesh.
which is a form of equation (3) in the three not-independent axes x, I and m. Considering now the Z-transform (with T = 1) of the velocity signal for an arbitrary junction j, we have from equation (9) the two relations
(11) and, summing the six outgoing wave signals, 2
N
Vj(z) = N L
(12)
Vi-(z).
i=l
Letting Vj,i(t) be the velocity signal at the i-th junction adjacent to the junction j, we can rewrite equation (12) in a form involving the wave signals ingoing to the adjacent junctions at the following time sample: 2 N Figure 5. Interferenceof two wavefrontsof velocityon the triangular mesh.
method exists for calculating the dispersion error3. The proof considers velocity signals and waves. The partial derivatives computed along the axes I and m are in relation with the partial derivatives along the main orthogonal axes through the equations
Vj(z) = N L
zVj,i+(z),
(13)
i=l
where Vj,i+ (t) is the wave signal going to the i-th adjacent junction, and coming from junction j. Similarly, we can rewrite equation (11) involving again the signals at the adjacent junctions, with obvious meaning of Vj,i- (t) :
(14)
av = ~ av + J3 av 8l 2 ax 2 ay av 1 av J3 av . { am - 2 ax - 2 ay 3 The same results on the analysis of dispersion in a triangular mesh have been obtained independently,with a similar procedure, by VanDuyne and Smith [19]
Multiplying equation (13) by Z~2 and adding to it both members of equation (14), we obtain the equation 2
(1+z~2)Vj(z)=
N
NLz-1Vj,i(Z). i=l
. acta acustica
ACUSTICA
Fontanaand Rocchesso:Physicalmodelingof membranes
Vol. 84 (1998)
The inverse Z-transform of it, after a multiplication by z, is
+ 1) + vj(n -
vj(n
2 1) = N
N
2: vj,;(n),
(15)
;=1
2
2 N
2
3" Co
Using these definitions, we can rewrite equation (15) at time sample nT for the junction at the position (x, l, m) in the form
+ T,x,l,m)
vj(nT
where we recognize the left member to be the second-order finite time difference, and the right member to be the secondorder finite space difference along the axes x, land m, with their central terms omitted. The comparison between equation (15) and equation (10) proves the equivalence of the mesh algorithm with a finitedifference scheme. We also notice that signal travels at a nominal speed Co = 1/ y2 spatial samples per time sample (or, equivalently, coFs waveguide units per second). In fact, equation (15) is a finite-difference version of equation (10) where the speed of propagation on the membrane equals the speed of propagation on the mesh, such that 1 3
-vj(nT, x + cxT, l, m) + -vj(nT, x - cxT, l, m) 3 3 1
V2cT/~m:S
~l
=
(16)
1
= ~m =
V2cT
thus ensuring that the conditions (16) are satisfied in the equality. The benefits of this fact are twofold: on the one hand the discrete system is stable, and on the other hand the numerical dissipation and dispersion turn out to be minimized [20]. Now we compute the dispersion error of the triangular mesh by measuring the dispersion factor k(f) = c/(coFs). First, we can see the velocity signal corresponding to the mesh junctions as a function vj (nT, x, l, m), defined in the discrete-time and in a non-orthogonal discrete-space. For the linear dependency of x, land m, the velocity signal on a junction is represented by many values of this function, anyway a complete knowledge of it allows the definition of a new, non-ambiguous function4 Vj (nT, X, fj) in the independent coordinate axes (x, fj) through the relation _
vj(nT, x
1
1
J3
2
2
2
+ -l + -m, -l
J3
- -m) 2
1
In
E Z(dx)
x Z(dy)}
, M
1
+ -vj(nT,x,l,m-cmT),
+ -vj(nT,x,l,m+cmT) 3
3
or, using the function vj, in the form
+ T, x, fj) + vj(nT - T, x, fj) =
vj(nT 1
1 x + cxT,fj) + -vj(nT, x3 _ J3
-vj(nT,
3 1_ _ 1 + -v 3) ·(nT , x + -c-T 2x, ,x
+ -v(nT 3)
-
1_
_
1
x-
+ -v 3) ·(nT , x
y
+ -c-T) 2 y
J3
-c-T 2x,
1
y- -
1
_ J3 y - -c-T) 2 y
1
J3 y- + -c-T) 2Y'
+ -c-T 2x, -c-T x 2,
cxT, fj)
-c-T) 2 y
where c( *) are components of speed along the direction ofthe axis (*) such that C( *) T equals the waveguide length d( *). The last equation may be Fourier transformed over the space domain, so that it gets mapped into a function defined in the discrete-time and in the spatial frequencies (6;, ~jj). After the application of Euler's formulas, this function becomes:
Vj(nT
1/2 I = 10 .../3/2 ' 1
suchthat we can definefor it the discreteFouriertransform [21].
+ T,~x,~jj) + Vj(nT - T,~x,~jj)
(17)
= bVj(nT,~x,~jj), with b=
2
3" cos (27r~xcxT)
+ J3 c-c-) T} 2 ero
xxoe
I
60
4.6. Tuning the Filter Solving equation (27) in the term
J
ere) mn
1+
tJ. (f(id») '7
200
400
600
Hz
800
1000
1200
1400
jj;;~) results in Figure 16. Shifting of the ideal modes (dashed line) into the modes (solid line) obtained using the loaded junction. '0' and' x' denote the calculated mode positions respectively of the unloaded and loaded membranes.
(id)2 j mn
=
L
o
,
mn
'7
and, substituting equation (29), we have -10
(id)2
(re) mn
J
I I
j mn
= 1 + ko
(id»)2 f ( fmn
nd)
-20
+ mo
In order to define a harmonic series, we can set ko = -0.8 and mo = 0.15. With this choice, the last equation results in Table III. The first six axial modes define a harmonic series, with adjacent modes separated by a distance of 1/4 of the fundamental mode 11. If central modes are forced to decay faster, as indicated in Section 4.1, the model, if excited near the rim, produces a pitched sound after an inharmonic attack. If excited near the center, the model produces a short inharmonic sound, since the axial modes are not excited. This is what physically happens in many percussion instruments, particularly in kettledrums, even if the harmonic series is different, as it can be seen comparing Tables III and II.
In Figure 16, the plot in solid line is the frequency response of the model having its Loaded Junctions set as seen above. The calculated positions of the modes are marked with' x '. The frequency response of the ideal, unloaded model is also given, in dashed line, with the calculated positions of the ideal modes being marked with '0'. The correspondence between the calculated and the effective shifted mode positions is good. An accurate look reveals that this correspondence is best for the less shifted modes, i.e. the modes affected by low (positive or negative) air load, ranging from the second to the tenth mode. Position error is caused by the bilinear transform, which introduces a phase distortion when turning the load filter into the discrete-time formulation. Tests on the filter have shown that this error is approximatively proportional to the air load. It is always possible to adjust the warping of the frequency
-30, I I
-40 [D "0
I
-501
-70-80
o
200
400
600
Hz
800
1000
1200
1400
Figure 17. Frequency distribution of the modes of the loaded and attenuated membrane. 'x' denotes the calculated mode positions.
axis introduced by the bilinear transform by changing the factor 2Fs to some other value. In our application we use the standard 2Fs factor, corresponding to integration by trapezoidal rule [24], since the overall performance is satisfactory for our purposes. Providing the system with the attenuation model described in Section 4.1, and observing it after the central modes have decayed, the harmonic series may be clearly observed, as the plot of Figure 17 shows. The calculated positions of the axial modes are marked with 'x'.
5. Excitation of the Membrane The way to excite a vibrating membrane-for example using a drumstick rather than a felt mallet-influences greatly the
ACUSTICA·
acta acustica
Fontanaand Rocchesso:Physicalmodelingof membranes
Vol. 84 (1998)
TableIII. Modeshiftingwith k
= -0.8
and m
541
= 0.15.
mn
01
11
21
02
31
12
41
22
03
51
32
61
f::;~)/ fi~e)
0.97
1
1.24
1.31
1.49
1.62
1.75
1.92
1.98
2.01
2.22
2.26
150
140
130
120 aJ ""0
110
Figure 18. Modelof hammer-membraneinteraction. 100f
sound obtained by it. The same influence has the exciter used in the membrane model. Physical models of exciters are usually discrete-time filters, modeling simple mechanical systems made of a mass and a spring, coupled with the resonator through waveguide links, as shown in Figure 18. Many studies have been devoted to the simulation of percussive exciters, especially within the framework of piano modeling [10,25, II, 12, 13,7]. These systems are commonly named hammers, since they were first designed to model the piano hammer. The interaction between hammer and string involves non-linear phenomena, which are summarized in the model using a nonlinear spring. The behavior of the system is mainly dependent on how the non-linearity has been handled. Most of the models cited above suffer from some problems, especially in real-time applications. Namely, they can be unstable if the sampling frequency is too low, or their parameters are difficult to tune. A Hammer-String Interaction model which overcomes these difficulties has been recently proposed by Borin and De Poli [14, 9]. We coupled this model to the triangular mesh, thus having a stable exciter whose parameters are easy to tune. The model is based on a non-linear spring-mass model defined by the relations
(31)
where f is the force applied on the spring, Xh is the displacement of the hammer, and ~x the spring compression. The hammer action can be controlled by means of the mass mh, the spring constant kh, and p, which defines a non-linear map since it is usually larger than one. The system of equations (31) is coupled to the waveguide model ofthe membrane by equating the sum of the forces at the branches of a junction to the force governing the hammer. A discrete-time solution is found by using the bilinear-transform approximation for time derivatives, and by solving the implicit nonlinear func-
90 1
80
o
l_
o
200
400
600
00
800
0
0
0
OQ)
00
1000 1200 1400 1600 1800 2000
Hz
Figure 19. Frequency spectrumof the transient attack using a stiff hammer. '0' denotes the calculatedmode positions.
tion of f (n) [26]. This latter step can be carried out in closed form if we use p = 2, a value found to be good in practice. The following steps result: D,
• w(n) =
~x(n - 1) T
1
+6"
+ vh(n
- 1)
6
L (vi+(n)
+ vi+(n
- 1)) - f(n - 1),
i=l
- 1 + 2abw(n) - Jl + 4abw(n) ( ( )) • f( n ) 2ab2 1w n , T • vh(n) = vh(n - 1) - (f(n) + f(n - 1)), 2mh where
Vh
is the hammer velocity,
a ~ kTP,
b ~ ~ (2~h
-
6~J'
and 1(.) is the unit step function. It is possible to show that the contact condition for the hammer with the string can be tested on the sign of the pseudo-compression w(n) [9]. Figure 19 shows the spectrum of the transient attack signal, exciting the standard mesh with a stiff hammer similar to a drumstick (mh = 40 g and kh = 107 Nfm). Figure 20 shows the spectrum of the attack signal using a soft hammer similar to a felt mallet (mh = 70 g and kh = 3.106 Nfm). The number of hammer models to be applied to the waveguide mesh depends on the performances requested from the simulated membrane. Two hammers, one located at the center and one located near the rim of the mesh, may be enough for a test of the model. Hammers add computations. In particular, at least two multiplications per time sample are always needed for each
ACUSTICA·
542
Fontana and Rocchesso: Physical modeling of membranes
150
140"
130f 1
120f
rg
1
1101 I
100I
90 80-·~~-
o
200
-- -"400
600
800
~ ---
-~---
~--
1000 1200 1400 1600 1800 2000
Hz Figure 20. Frequency spectrum of the transient attack using a soft hammer. '0' denotes the calculated mode positions.
hammer. In addiction, one access to a look-up table is needed at each time sample for any active hammer (in contact with the waveguide mesh).
6. Results and Ongoing Research Simulations of the membrane model have been conducted under different playing conditions. The sound samples obtained are certainly encouraging, although even better results might be obtained by modeling a complete percussion instrument. Samples can be found via the WWW page of one of the authors: http://www.dei.unipd.it/ricerca/csc/people/roc. Possible topics for further research include the improvement of the dissipation-filter design procedure, in order to obtain a precise attenuation for the central modes. Higher order transfer functions of load filters may allow to model more complex behaviors of real membranes. Finally, a single filter modeling all the attenuation effects, the mode shifting effects and the excitation, might be designed.
References [I] J. O. Smith: Physical modeling synthesis update. Computer Music Journal 20 (1996) 44-56. [2] J. O. Smith: Physical modeling using digital waveguides. Computer Music Joumal16 (1992) 74-91. [3] S. A. Van Duyne, J. O. Smith: Physical modeling with the 2D digital waveguide mesh. Proc. of the International Computer Music Conference, Tokyo, Japan, ] 993.
acta acustica Vol. 84 (1998)
[4] S. A. VanDuyne, J. O. Smith: The tetrahedral digital waveguide mesh. Proc. of the IEEE ASSP Workshop on Applications of Signal Processing to Audio and Acoustics, Mohonk, New York, 1995. [5] S. A. Van Duyne,J. O. Smith.: The 2D digital waveguide mesh. Proc. of the IEEE ASSP Workshop on Applications of Signal Processing to Audio and Acoustics, Mohonk, New York, 1993. [6] L. Savioja, V. Viilimiiki: Improved discrete-time modeling of multi-dimensional wave propagation using the interpolated digital waveguide mesh. Proc. Int. Cont. on Acoustics, Speech and Signal Processing, Munich, Germany, 1997. [7] S. A. Van Duyne, J. R. Pierce, J. O. Smith: Traveling wave implementation of a lossless mode-coupling filter and the wave digital hammer. Proc. of the International Computer Music Conference, Aarhus, Denmark, 1994. 411-418. [8] N. H. Fletcher, T. D. Rossing: The physics of musical instruments. Springer-Verlag, New York, 1991. [9] G. Borin, G. De Po\i: A hysteretic hammer-string interaction model for physical model synthesis. Proc. Nordic Acoustical Meeting, Helsinki, Finland, June 1996. 399-406. [10] A. Chaigne, A. Askenfelt: Numerical simulations of piano strings I. A physical model for a struck string using finite difference method. J. Acoust. Soc. Am. 95 (1994) 1112-1118. [II] D. E. Hall: Piano string excitation. I: nonlinear modeling. J. Acoust. Soc. Am. 92 (1992) 95-105. [12] A. Stulov: Hysteretic model of the grand piano hammer felt. J. Acoust. Soc. Am. 97 (1995) 2577-2585. [13] H. Suzuki: Model analysis of hammer-string interaction. J. Acoust. Soc. Am. 82 (1987) 1145-1151. [14] G. Borin, G. D. Poli: A hammer-string interaction model for physical model synthesis. Proc. of the XI Colloquium on Musical Informatics, Bologna, Italy, 1995. 89-92. [15] P. M. Morse: Vibration and sound. American Institute of Physics for the Acoustical Society of America, New York, 199 I, 1st ed. 1936, 2nd ed. 1948. [] 6] J. O. Smith: Music applications of digital waveguides.CCRMA - Stanford University. Report STAN-M-39, 1987. [17] D. E. Dudgeon, R. M. Mersereau: Multidimensional digital signal processing. Prentice Hall, Inc., Englewood Cliffs, NJ, ]984. [18] V. Comincioli: Analisi numerica: Metodi, modelli, applicazioni. McGraw-Hill Libri Italia, Milano, 1995. [19] S. A. Van Duyne, J. O. Smith. Personal communications, Fall ]995. [20] A. Chaigne: On the use of finite differences for musical synthesis. Application to plucked stringed instruments,. J. Acoustique 5 (1992) 181-2] I. [21] G. Cariolaro: La teoria unificata dei segnali. UTET, Bologna, 1991. [22] T. I. Laakso, V. Viilimiiki, M. Karjalainen, U. K. Laine: Splitting the unit delay - Tools for fractional delay filter design. IEEE Signal Processing Magazine 13 (1996) 95-105. [23] S. A. Van Duyne, J. O. Smith: A simplified approach to modeling dispersion caused by stiffness in strings and plates. Proc. of the International Computer Music Conference, Aarhus, Denmark, 1994. 407-410. [24] A. V. Oppenheim, R. W. Schafer: Discrete-time signal processing. Prentice-Hall, Englewood Cliffs, NJ, 1989. [25] G. Borin, G. De Poli, A. Sarti: Sound synthesis by dynamic systems interaction. - In: Readings in Computer-Generated Music. D. Baggi (ed.). IEEE Computer Society Press, 1992, 139-160. [26] F. Fontana: Modelli di membrane vibranti realizzati con reticoli waveguide. Dissertation. Universita di Padova, Italy, 1996.
