1. Introduction

Physical sound synthesis uses mathematical models based on the physics of sound production to synthesize sound. In other words, physical sound synthesis uses a model to simulate the sound producing object, rather than the sound produced by an object. In this paper we focus on mass-spring systems: networks of masses, springs, and dampers. The mathematical model of mass-spring systems is based on differential equations. To approximate these differential equations on a digital computer, numerical methods are used. An important question to ask when using a numerical method is how well does this method approximate the differential equations used in the system? Numerical methods can become unstable. This means that the numerical solution can deviate arbitrarily far from the exact solution. In many cases the error can grow without bound, making the results of the numerical method meaningless.

As well, we want to be able to quantify the accuracy of our approximation. Most musical sounds are composed of a number of different frequencies. The lowest of these frequencies is called the fundamental frequency [1]. The fundamental frequency determines the perceived pitch of the sound, and an error in the fundamental frequency will cause the sound to be out of tune. The higher frequency components influence the timbre, or tone colour, of the sound [1] and errors in these components give the sound a different timbre than it should. An error in frequency caused by a numerical method is known as frequency warping.

The decay rate determines how quickly the amplitude (or volume) of the sound decreases. For example, a note on a piano can be heard for 20 or 30 seconds after it is struck, while on a banjo it becomes imperceptible after only 3 or 4 seconds because the decay rate of a banjo is much larger than that of a piano. Numerical methods may add extraneous damping to vibrating systems that are undamped. This is known as numerical damping.

There have been several sound synthesis systems built using mass-spring systems and described in the literature [2–6]. Most of these systems have used a numerical method called the symplectic Euler method. The previous literature on mass-spring systems used in sound synthesis describe how these systems work (i.e., the equations used and the finite difference equations used to approximate them), but have not addressed the issues of the stability and accuracy of the numerical methods. This has been an important part of the criticism of mass-spring systems in the physical synthesis literature. Mass-spring system have been criticized as being computationally expensive [7], lacking an analysis of stability [7], and having an unknown accuracy [8].

The symplectic Euler method has been studied by researchers outside the sound synthesis community. The book by Hairer et al. [9] presents a thorough analysis of symplectic numerical methods including the symplectic Euler method. Using different methods—they do not use the -transform—they arrive at the stability condition for the undamped mass-spring system as which is the same as our results when the damping is zero. They do not analyze the symplectic Euler method in terms of frequency warping or its effect on damping. The thesis by Beck [10] contains a proof that symplectic Euler method is symplectic, which implies that it has no numerical damping for an undamped mass-spring system, which agrees with our conclusions. Beck's thesis looks at using the symplectic Euler on the Lotka-Volterra (predator-prey) equations. It also does not analyze its frequency warping or its effect on damping. In a previous paper [11] we have given an analysis of the symplectic Euler method when used to simulate undamped mass-spring systems. In this paper, we extend that analysis to include damped mass-spring systems.

The contributions of this paper are the equations for the stability, the frequency warping, and numerical damping of damped mass-spring systems presented in Section 3. This paper is not presenting a new method of doing sound synthesis or suggesting improvements to existing methods. The question it proposes to answer is, if given a specification of a mass-spring system, that is, the values for the mass, spring stiffness, and viscous damping constants, and the connections between them—will the system be stable, and if so, what sound will it produce? The three main questions addressed in this paper are the following.

Section 2 introduces the mass-spring system and explains why the symplectic Euler method is often used to discretize the differential equations of a mass-spring system. Section 3 presents the analysis of the symplectic Euler method. We begin by using the symplectic Euler method to discretize a mass-spring system containing only one mass. We use the -transform to find the symplectic Euler method's effect on the frequency and decay rate of the system, and find the conditions for stability of the system. This section contains the main contributions of the paper. We end this section by demonstrating the consistency of our theoretical results with the results of a computer simulation of a mass-spring system. In Section 4 we show how the results in the previous section can be extended to mass-spring systems with more than one mass. Section 5 concludes with a summary of our results.

2. The Mass-Spring Model

The mass-spring model builds complex musical instruments from simple components: masses, springs, and dampers. Each element is discretized using finite difference methods. The behavior of the system depends solely on the network and the physical equations of each of the components. No other physical equations are used.

Figure 1 shows a simple mass-spring model, where and are masses, , and are springs and is a damper.

Figure 1: Simple mass-spring system.

2.1. Choosing a Numerical Method for Mass-Spring Systems

To simulate the vibrations of a mass-spring system on computer, we need to use a numerical method to discretize the differential equations of the system. There are many numerical methods to choose from. What should we look for when choosing a numerical method for sound synthesis? Humans can hear sounds that have frequencies from 20 to 20 000 Hz. Notes played on typical musical instruments such as a piano, a guitar, a trumpet, and so forth. may last for several seconds. This means that a note may contain thousands or tens of thousands of cycles. The energy of a mass-spring system is the sum of its potential and kinetic energies, which depend on the amplitudes of the vibrations. It is therefore important that a numerical method used in sound synthesis can be able to conserve energy for thousands of cycles. If the numerical method causes the energy to increase over time, the simulation will become unstable. Conversely, if the numerical method causes the energy to decrease, the sound will decay more rapidly than it should. Numerical methods that do not conserve energy have proved to be a problem in fields such as molecular [12] and planetary simulation [13]. There has been an interest in recent years in numerical methods that can accurately simulate the qualitative aspects of physical systems. Symplectic numerical algorithms, among other properties, conserve energy over long periods of time [14].

Figure 2 shows an undamped mass-spring system containing one mass and one spring.

Figure 2: Undamped mass-spring system.

If we regard the equilibrium position of the spring to be , the force of the spring, according to Hooke's Law, is , where is the spring stiffness coefficient and is the position of the mass at time . We can then write the differential equation for the system, using for the acceleration and for the mass, as This is a second-order differential equation. The general solution is [15]

Setting initial condition (the initial position) and (the initial velocity) the particular solution is where is the radial frequency of the system. So the solution of the system is simply a cosine wave of frequency . Since there is no damping the amplitude of the cosine wave should not change over time (i.e., the system should conserve energy). Next, we examine how well this system is simulated by three first-order numerical methods: the forward Euler, the backward Euler, and the symplectic Euler. In these simulations we set at radians per second and the sample rate at 1000 samples per second (i.e., the frequency is the sample rate).

The forward Euler method is defined as where is the first derivative of with respect to time and is the length of each time step. We use the subscript notation to represent numerical approximations (e.g., denotes the numerical approximation of at time step ). The forward Euler method gains energy over time, causing it to be unstable for undamped or lightly damped systems. Figure 3 shows the result of simulating the above mass-spring system using the forward Euler method with equal to the sampling frequency. We can see that this results in an unstable system.

Figure 3: Forward Euler approximation of the undamped mass-spring system ().

The backward Euler method is defined as The backward Euler is called an implicit method since it uses the derivative at the new point which has not yet been determined. The backward Euler method loses energy over time. Figure 4 shows the result of simulating the above mass-spring system using the backward Euler method with equal to the sampling frequency. We can see that this results in cosine wave that is quickly damped.

Figure 4: Backward Euler approximation of the undamped mass-spring system ().

The symplectic Euler is defined as where is the displacement, is the velocity, and is the acceleration. We first calculate the new velocity , since it can be calculated using the known acceleration and the known velocity . We can then use to calculate the new position, . This makes the entire system explicit. This method is therefore sometimes referred to as the explicit version of the symplectic Euler [10], since there is also an implicit version [9, 10]. Since the symplectic Euler method, as its name implies, has been shown to be a symplectic numerical method [9, 10] unlike the forward and backward Euler methods, it should conserve energy. Figure 5 shows the result of simulating the above mass-spring system using the symplectic Euler method with equal to the sampling frequency. We can see that the results of this simulation are much nearer to the analytic solution than either the forward or backward Euler and that the amplitude of the vibration appears to be constant. The ability of the Symplectic Euler Method to conserve energy is probably the main reason why it has been used in several of the mass-spring systems built for sound synthesis, such as the CORDIS system [4–6] and the TAO system [2, 16]. The symplectic Euler is a first order method, meaning that its global truncation error (the cumulative error) is proportional to the time step . In a previous paper [11] we explored the possibility of using higher order symplectic methods. We concluded that in cases where the mass-spring system is being used to simulate a continuous system such as a vibrating string or two dimensional membrane, the extra accuracy of the higher order method is not worth the increase in computational cost. When a small number of masses are used to model a continuous system, the resulting mathematical model is not very accurate, and increasing the accuracy of the numerical method does not noticeably improve the sound of the simulation. If a large number of masses are used to model the continuous system, the model is much more accurate, but the sound produced will contain high-frequency components. In order to keep the system stable and to avoid aliasing, the sample rate has to be increased for both the symplectic Euler and the higher order methods. At a high sample rate the symplectic Euler method is quite accurate for the low and medium frequency components of the sound, and so there is still no appreciable difference in the sound of the simulation between the first-order symplectic Euler and higher order symplectic methods. This makes the symplectic Euler method a good choice for most sound synthesis applications using mass-spring systems. But it is still important to know the accuracy and stability limits of the method in order to set the sample rate of the simulation and resolve problems in situations where the sound produced is not as expected.

Figure 5: Symplectic Euler approximation of the undamped mass-spring system ().

2.2. Mass-Spring Discretization

We now use the symplectic Euler method to discretize the mass-spring model.

Mass Element
We can derive the behavior of a mass from Newton's 2nd law: where is the force acting on the mass at time is the mass, and is the acceleration of the mass at time . Since acceleration is the derivative of the velocity we can write (8) as We then have a system of two first-order differential equations: We then use the symplectic Euler method to discretize these equations. Substituting in (7) gives us where is the displacement, is the velocity and is the numerical approximation of the force acting on the mass at time step .

Note that if was used instead of in (7) we would have the backward Euler approximation. This is the way this numerical method has been described in the sound synthesis literature [2–4]; as the backwards Euler method with the forces delayed by one time step, rather than the symplectic Euler method.

Spring Element
The equations for the spring are derived using Hooke's Law , where is a constant denoting the spring stiffness. We write them as Here we let , and represent the distance from the equilibrium position of mass and mass at either end of the spring. We use to denote the force acting on mass at one end of the spring at time step . The force, , acting on mass at the other end of the spring is, according to Newton's third law, equal and opposite to .

Damper Element
The damper element is used to represent viscous friction. This is the object's resistance to motion and is assumed to be proportional to the velocity. The formula for the damper is , where is a constant denoting the coefficient of viscosity, the force, and the relative velocity of the two ends of the damper. This can be written as where and represent the forces acting on the masses at the ends of the damper, and and are the corresponding velocities.

2.3. Synthesizing Sound Using a Mass-Spring System

The mass-spring system works by discretizing the physical equations of each of the elements—the masses, springs, and dampers—of a mass-spring system such as the one shown in Figure 1. At each time step the sums of the forces acting on each of the masses are calculated. These sums consist of the forces of the springs and dampers directly connected to the mass and any external forces acting on the mass. The external forces are used to simulate various physical interactions with the instrument, such as plucking, hitting, and bowing. The forces are used to calculate the new positions and velocities of the masses. Once the new positions and velocities have been calculated, they are used to calculate the new forces acting on each mass. This cycle then repeats for the duration of the simulation. The position of one of the masses at each time step in the simulation is written to a sound file to represent the sound produced at this point in the simulated instrument. Alternatively, the vibrations of several masses can be summed together and written to the sound file.

Since each element interacts only with the elements connected to it, the number of calculations depends on the number of masses and average number of connections of each mass. For simple systems these calculations can be done in real time [17]; more complex systems must be run offline. More detailed accounts of the implementation of mass-spring systems can be found in Cadoz et al. [4], and Pearson [2].

3. Analysis of a Damped Mass-Spring System Systems with a Single Mass

In this section we look at the stability and accuracy of the symplectic Euler method when used to simulate a damped mass-spring system containing a single mass, a single spring, and a single damper. We start by finding the analytical solution of this system. We find the -transform of this equation when it is approximated by the symplectic Euler method. We then find the damping and the frequency of the discrete mass-spring system represented by transformed equation. We also find the conditions under which this system is stable.

3.1. The Analytical Solution of the Single Mass Damped Mass-Spring System

Figure 6 shows a damped mass-spring containing only one mass, , one spring, , and one damper, .

Figure 6: Damped mass-spring system.

The equation for this system is where is the distance of the mass from its equilibrium position, is the mass, the spring stiffness, and the viscous damping coefficient.

We can solve this equation by finding the roots of the characteristic equation We use the substitutions and : The roots of this characteristic equation are For the system to vibrate we require that , so The condition dividing vibrating systems from those that do not vibrate occurs when This value is known as critical damping.

The general solution of (15), when , can be shown to be [15] where and are constants depending on the initial conditions. This can be written as [15] where This shows that the damped mass-spring system has a starting amplitude of . This amplitude is being decreased by the term . The frequency (actually the quasi-frequency since the system is not strictly periodic) is . As the damping approaches zero this equation becomes .

3.2. The Damped Mass-Spring System Using the Symplectic Euler Method

We next examine the damping and frequency of the damped mass-spring system, when approximated by the symplectic Euler method, by using the -transform.

The equation of the damped mass-spring system shown in Figure 6 is where is an external force acting on the mass. Using the substitutions from the previous section we can write (24) as We now discretize this equation using the symplectic Euler method by substituting (25) in (7): We should note that (26) can also be derived by substituting the equations for the spring (12) and the damper (14) together with an external force in the equation for the mass (11). This is how mass-spring systems actually work: by calculating the new values of each mass, spring, and damper using the equations from Section 2.2 at each time step. We can write (26) as a scalar equation by substituting the second line in the first line: Since, from the first line of (26), , we can write as . Using this substitution in (27) gives us which simplifies to We can shift the time step back by one to get

3.3. The -Transform

The -transform takes signals from the time domain and transforms them to signals on the -plane. The -plane represents the signals in terms of amplitude growth (the growth or decay rate of the signal) and frequency, which are particularly useful in the analysis of musical signals. The -transform is defined [18]: where . The -transform of a system can be represented on the -plane, which is shown in Figure 7. The -plane uses polar coordinates with being the distance from the origin and the angle. The amplitude growth or decay of the signal is represented by and the frequency by . The frequency varies from on the positive real axis, to radians per sample on the negative real axis. This frequency, also known as the Nyquist limit, is the maximum frequency a discrete system can have. If the frequency goes above the Nyquist limit, it becomes indistinguishable from a frequency less than the Nyquist limit and the resulting frequency will be perceived as the lower of the two frequencies. This is known as aliasing, and causes inaccuracies in synthesized sound. The transfer function of a system is defined as the -transform of its output divided by the -transform of its input. The poles of a system are defined as the roots of the denominator of its transfer function. The poles are the normal modes or the natural frequencies of the system. A system is called stable if, when the input is absolutely summable, the output is absolutely summable. The system is stable on the -plane if all its poles lie inside the unit circle [18] and marginally stable if it has a pole on the unit circle, but no pole outside the unit circle. A marginally stable system has a bounded output in some conditions, such as when the system has no input; but oscillations in a marginally stable system do not die away but persist indefinitely [19]. If any poles are outside the unit circle, the system is unstable. An important feature of the -transform is that multiplying the -transform of a signal by is the same as delaying the signal by time steps [20], that is,

Figure 7: The -plane.

3.4. Using the -Transform to Find the Poles of the System

Using the property of (32), the -transform of (30) is and the transfer function is The poles of the transfer function occur when The roots of (35) are

3.5. Using the Poles to Analyze the Mass-Spring System

The poles represent the modes or natural frequencies of the system. If the poles are complex, the discrete system will vibrate at a frequency in the range radians per sample. In this section we first determine what the conditions are for complex poles. We then find an equation for the damping of the discrete system when it has complex poles. The damping is determined by the distance of the pole from the origin on the -plane. We then determine the frequency of the discrete system by finding the angle of the poles on the -plane.

If the poles are real, then the frequency of the discrete system is either zero, if the larger magnitude pole is on the positive real axis, or radians per sample, if the larger magnitude pole is on the negative real axis. We find the conditions determining whether the larger pole is positive or negative. We also find conditions for a pole being outside the unit circle, which will make the discrete system unstable.

We consider 2 cases in (37): is complex and is real.

3.5.1. Case Is Complex

We look at the case in which contains an imaginary component, that is, when We first calculate the range of values for of this case. If by multiplying both sides by and rearranging. The roots of this equation are This simplifies to If is complex, falls between these two roots.

We now look at the damping for case ( is complex). In this case (37) becomes We can calculate the length of using the imaginary and real parts: where denotes . This simplifies to We see from (44) that, for this case, the length of does not depend on the frequency.

Figure 8 shows an example of the poles on the -plane of the symplectic Euler approximation of the damped mass-spring equation. The value for used for is . Root 1 and root 2 are the two roots of (37): root 1 has the surd added, while in root 2 it is subtracted. The arrows show how these roots—the poles of the system—move as varies from zero to above radians per sample. The poles start on the positive real axis where the frequency is zero. While is complex, each pole traces a semicircle, root 1—where the surd is added—is the upper semi-circle and root 2—where the surd is subtracted—is the lower semi-circle. If becomes large enough, the solutions of (37) become real again—but now they are on the negative real axis. Root 1 moves toward the origin, while root 2 becomes increasingly negative. When root 2 moves outside the unit circle, the system becomes unstable. Figure 9 shows the same plot with at . The radius of the circular region, where is complex, is now very small.

Figure 8: Damped mass-spring system using symplectic Euler-poles on -plane with

Figure 9: Damped mass-spring system using symplectic Euler-poles on -plane with

The radius of the circle containing the poles is , which depends on both and the time step . If the damping or the time step increase, the radius of the circle will become smaller. If the system is undamped (i.e., ), the circle will have a radius of one. This means that for the undamped system there is no numerical damping. This is consistent with the fact that the symplectic Euler method conserves energy.

We next look at the frequency warping when is complex. From (42), we can calculate the frequency using the real and imaginary components. Using to denote the actual frequency obtained using the symplectic Euler,

Figure 10 shows the frequency warping for the undamped mass-spring system. Note that for the undamped system the low frequencies are very accurate, and the range of , where is complex, goes from to radians per sample. As increases, the digital frequency (the actual frequency produced by the symplectic Euler approximation) becomes increasingly warped upward. When has reached radians per sample, the digital frequency is radians per sample—the Nyquist limit.

Figure 10: Frequency warping for undamped system.

Figure 11 shows an example of the frequency warping of the symplectic Euler method when the damping is quite high: . The analog frequency is calculated as , and is slightly lower than . Note that the digital frequency has reached the Nyquist limit of radians per sample at around radians per sample. This is consistent with (41), which gives as the lower limit and as the upper limit for being complex. As increases, the range where is complex decreases and the frequency warping becomes more pronounced.

Figure 11: Frequency for damped mass-spring system, .

3.5.2. Case Is Real

Equation (41) gives us 2 conditions for to be real. Equation (36) yields 2 roots, which are the two poles of the system. The pole with the larger magnitude will dominate the system. If the pole with the larger magnitude is real and positive, the frequency of the system is zero. If it is real and negative, the system will vibrate at radians per sample.

We show that, if is greater than the root of (41) in which the surd is added, then (36) will have its largest magnitude root less than zero. Since the surd being added or subtracted in (36) is the same, what we want to show is simply that

We also assume that is less than . The case where is considered later. We can write the left side of (46) as From these facts we can derive that Since the pole with the larger magnitude is on the negative real axis of the -plane, the frequency is radians per sample.

Similarly, if is less than the root of (41) in which the surd is subtracted, that is, if then So in this case, since the pole with the larger magnitude is on the positive real axis of the -plane, the frequency is zero.

As the leftmost pole approaches , the damping becomes smaller. The system will become unstable when becomes less than . This happens when from (37) Squaring both sides and simplifying results in For a fixed time step, , we can solve for : So for stability we require that

As the value of increases, the circular region on the -plane where is complex becomes smaller, as can be seen in Figure 9 where . When the circular region disappears. For , there are only two possible digital frequencies: zero and radians per sample. If the pole with the larger magnitude is negative, the digital frequency will be radians per sample, otherwise it will be zero. From (37), the pole with the larger magnitude is negative when so, when , the curve dividing systems that do not vibrate and those that vibrate at radians per sample is

3.6. Regions of the -Plane

We are using a mass-spring system as a way of mathematically modelling a vibrating physical system. We then use numerical methods to approximate this mathematical model and implement a simulation of the system on a computer. There are 2 sources of error.

In this paper, our concern is with the second error. We want to know, given the mathematical model, how accurate is the computer simulation? The mathematical model, which is a continuous model, can be analyzed on the -plane. Analogous to the -plane, the -plane represents the poles of the analog (continuous) system where . The -plane uses rectangular coordinates with the horizontal axis—the real axis—representing the decay rate and the vertical axis—the imaginary axis—representing the frequency. The decay rate is denoted by and the frequency by . The mathematical model, which we refer to as the analog system, is stable if all its poles are on the left hand side of the -plane. Ideally, all stable analog systems would result in stable computer simulations. However, this is not necessarily the case, since numerical methods are approximations. In this section, we show which systems having stable mathematical models will be stable when simulated using the symplectic Euler method. We do this by showing graphically which parts of the left hand side of the -plane will be mapped to stable systems on the -plane, that is, which poles on the -plane will, when the system is simulated using the symplectic Euler method, result in poles within the unit circle on the -plane.

We can find the region on the -plane that maps to stable poles on the -plane by using (54). In this case, however, we need it in terms of and . First we find in terms of and : We then substitute (57) in (53): The system is then stable if We can also solve (41) and (56) in terms of and . Equation (41) divides poles that are complex from those that are real. Substituting (57) into (41) we get Similarly, solving (56)—the curve dividing vibrating poles from nonvibrating poles—for in terms of results in Using these equations, we can divide the left half on the -plane into sections. These regions represent the qualitative properties that any pole on the -plane within the region will have when the systems is approximated with the symplectic Euler method. Figure 12 shows the regions of the -plane for positive frequencies. The negative frequencies are mirror images of the positive ones. If any of the poles of the analog system lie within the unstable region, the discrete system resulting from the symplectic Euler method will be unstable.

Figure 12: Regions of the -plane.

3.7. The Accuracy of the Damping

In this section we determine the accuracy of the damping of a mass-spring system discretized by the symplectic Euler method. We determine the equation for the digital damping and examine some of its points of interest. We find that it has a sharp change of direction when it moves from one region (from Figure 12) to another. It turns out that in some places, counter-intuitively, increasing the damping coefficient actually decreases the digital damping. We also find that the equation for digital damping has a singular point where the damping becomes infinite, and we determine exactly where that point is.

From (22), the damping of the analog system is . On the -plane this is . The damping of the discrete system is where . How do these two values compare? If we sample the continuous damping at each time step, we have as the samples of the analog damping, where is the sample number and the length of the time step. The discrete damping is where denotes the damping of the discrete system. So, Using the value of from (37) we can compare the analog to the digital damping. Figure 13 shows the damping of the symplectic Euler method when . When the analog damping is small (i.e., is near zero), it is quite accurate. Note that, since the accuracy of the damping depends on both and , we can increase the accuracy of the damping by decreasing the time step , because, for a fixed value of , decreasing will move the digital damping of the system toward the right side of the figure where the value of is very close to . There are two points when the function has a sharp corner: the first is when leaves the region where it is complex and enters the region where the frequency is zero; the second is where enters the region where the system vibrates at radians per sample. Once is in the region where the digital frequency is radians per sample, the digital damping decreases (i.e., approaches ) and the system becomes unstable when the digital damping is greater than zero. The positions of the vertical lines marking where the digital frequency becomes zero and radians per sample are found by solving (60) and (61), respectively, for in terms of , and finding the value of , given the known value of . For (61) this works out to For (60), solving for results in a very long fourth degree equation.

Figure 13: Digital damping versus analog damping—.

Solving for when results in the value for the line at which leaves the region where it is complex and enters the region where the digital frequency is zero, and for the line dividing systems that have digital frequencies of zero from those having digital frequencies of radians per sample. We can also see these values by drawing a horizontal line, for this example at , in Figure 12. We see that the points at which the horizontal line enters the regions and match those of Figure 13. Figure 14 shows the digital damping when . Solving for when results in the value for the line at which leaves the region where it is complex and enters the region where the digital frequency is zero, and for the line dividing systems that have digital frequencies of zero from those having digital frequencies of radians per sample.

Figure 14: Digital damping versus analog damping—.

Figure 15 shows the digital damping when . Solving for when results in the value for the line at which leaves the region where it is complex and enters the region where the digital frequency is radians per sample.

Figure 15: Digital damping versus analog damping— .

From (63) we see that when approaches zero, the digital damping approaches . This happens when from (36) If and , then