A Physical Model for Spectral Models Tutorial
(Modeling Sinusoids Numerically for Additive Sound Synthesis)
By Stephen McGovern, - -2007
email all questions and comments to:
draft version. Last recorded update: 04-25-2007

Preface and Table of Contents

The article is broken into the following sections:

1. Introduction 2. Simple Sinusoidal Sounds 3. Decaying Sinusoidal Sounds 4. Decaying Sinusoidal Sounds with Harmonics 5. Setting Initial Conditions 6. Reducing Computational Error 7. Conclusions 8. Bibliography

Most of the sections contain snippets of C++ code near the end. Below each snippet are links to complete and ready to compile C++ and Java code. It is suggested that some readers might prefer to review the code first and the related equations second.

1 - Introduction

There are many models that can be used to synthesize sound. Of these, physical and spectral models are common place. In a spectral model one has information about the frequency content of a sound. This information could convey how the frequency content changes over time. In any case this information is used synthesize sound.

In a physical model one had information about the physical process by which a sound is produced. This information is then used to run a simulation and synthesize sound.

This paper will present a synthesis method that takes spectral information about a sound and then produces a physical model that fits that information. This physical model is then run as a physics simulation. The method will, in essence, fit a physical model to a spectral one and thereby allow many spectral models to be implemented as physical ones.

2 - Simple Sinusoidal Sounds

Virtually all musical notes have a pitch. This pitch has a fundamental frequency that forms a simple sinusoidal waveform. Such a wave form is pictured below in Fig. 1.

To describe the above waveform in mathematical terms, we can use the following formula.

(1)

t is time, A is the peak amplitude, and is the angular frequency of the note. is given by

(2)

The variable f_note is the fundamental frequency of the note we want to synthesize. This frequency happens to be relatively easy to find. The standard reference pitch is A4 which has a frequency of 440 Hz ^1,7. What do we do if we want a frequency of a different note? Well, the frequency of each key on a piano, whether it be black or white, is related to the frequency of the previous one by the relation f_p+1 = f_p 2^1/12. The variable p is an integer index to all the keys on a piano. If f₀ = f_A4, then we can be generalized the relation with the equation

(3)

If our desired note is C2 we are 33 keys below A4. To find the frequency of C2, we solve Eq. 3 as follows.

(4)

We use a negative number because we are going from a higher frequency note to a lower one.

We now return to our discussion of our simple sinusoid. It happens that the equation we used to describe our waveform also the motion of the spring mass system pictured below.

This system is also described by the differential equation

(5)

and equivalently by Hooke's law

(6)

Where F = ma, or force = mass x acceleration. k is a property of the spring called the spring constant. It is given by

(7)

x is the position of the mass. If x = 0 the spring is neither stretched nor compressed. At this position the force it exerts on the mass is equal to zero. This position, however, is also the one for which the mass has its highest velocity. Conversely, when the mass is at the peak of it's oscillation, the force from the spring is at its highest, and the velocity of the mass is zero. Substituting F = ma in Eq. 6, we get

(8)

For simplicity, this article will assume that the mass m is always equal to 1 kg. Dividing both sides by m and neglecting units, we get

(9)

Where ever our mass is, we can find it's acceleration simply by plugging in a value for x. Now, let's say our mass starts off at time t = 0 with the initial conditions

(10a)

(10b)

(10c)

If we choose a very small time interval t = we can approximate the the values of x(t+), v(t+), and a(t+). It is easiest to take one time step for each audio sample. If our sample rate is given by f_{sample rate}, our time step is then given by

(11)

We can now find the position of the block at time t = by multiplying by our initial velocity and adding the result to our initial position. Thus the value of x at time t = is found to be

(12)

We can see from Eq. 9 and Eq. 10b that our acceleration was zero at x(0) = 0. Since x(t) has changed, we now have a new acceleration. Substituting in x() we get

(13)

Notice that the acceleration is dependent on x and not on time. This new acceleration in turn leads us to the following new velocity

(14)

The steps taken in Eq. 12, 13, and 14 can be generalized as follows

(15a)

(15b)

(15c)

Because we have found x(),v(), and a(x()), we can advance our time to t = 2 and use Eq. 15 to find x(2),v(2), and a(x(2)). By repeating these steps, we can approximate the values of x(t),v(t), and a(x(t)) at any point in the future.

Starting with the initial conditions from Eqs. 10 we can find future positions, velocities and accelerations by using the following simple C++ code. Note: is spelled out as "tau."

int main(){	
	double t60=2.0, w, k, x=0.0, a, v=565.5, tau, *Out, f_note=100.0;
	long length, sample_rate=44100, i;
	
	length = (long)( t60*(double)sample_rate );
	Out = new double[length];
	tau = 1.0/(double)sample_rate;
	w = 2.0*PI*f_note;
	k = w*w;
	
	for( i = 0; i < length; i++){
		x += v*tau;
		a = -k*x;
		v += a*tau;
		*(Out+i) = x;
	}
	write_sound(Out, sample_rate, length);	
	delete Out;

return 0;
}

Figure 3.   Numerical solution for a harmonic oscillator. Our time step is set to = 1/fsample rate. Where is spelled out as "tau." The spring constant, k, is given in terms of the frequency fnote.

Download full source code: Sine.cpp Sine.java

3 - Decaying Sinusoidal Sounds

Few sounds in the real world are described well by a simple sinusoid. Musical notes usually loose energy and die out in a relatively short amount of time. If we were to compare this to a mass-spring system again, this energy could be lost to friction. We could produce a friction force simply by submerging the mass in a fluid as pictured below.

This would result in a waveform that looked like the following

To describe this scenario mathematically, we can use by using Eq. 16.

(16)

A is a constant and b is a decay factor that relates to the note length. b is given by

(17)

where t₆₀ is the decay time of our note. It is reasonable to say that this decay time is the time length of the note we want to synthesize. It is referred to as the t₆₀ because the term e^-bt/2 in Eq. 16 will decay by -60 dB when t = t₆₀.

In section two the force on our mass was described by Hooke's law:

(18)

For a mass submerged in a fluid, we will again have spring force -kx. We will also have a drag force opposing the motion of the mass. Adding this force to Eq. 18 we get

(19)

Drag force is often related to velocity. Setting F_drag = bv as shown below should be sufficient.

(20)

The above force equation is equivalent to the differential equation

(21)

Recalling that out mass is equal to 1 kg and dividing it out of Eq. 20, we get

(22)

With the decay term our harmonic oscillator becomes a damped harmonic oscillator. A damped harmonic oscillator has a frequency lower than that of an undamped one. Because of this we have to reset the value of k to stay in tune. Our new value of k is given by Eq. 23.

(23)

Our generalized equations of motion from Eq. 15 can now be written with our new acceleration.

(24a)

(24b)

(24c)

Using Eqs. 24 we can replace the C++ code from Fig. 3 with the following.

int main(){	
	double t60=2.0, w, b, k, x=0.0, a, v=287.7, tau, *Out, f_note=50.0;
	long length, sample_rate=44100, i;
	
	length = (long)( t60*(double)sample_rate );
	Out = new double[length];
	tau = 1.0/(double)sample_rate;
	w = 2.0*PI*f_note;
	b = -2.0*log(0.001)/t60;
	k = w*w + b*b/4.0;
	
	for( i = 0; i < length; i++){
		x += v*tau;
		a = -k*x -b*v;
		v += a*tau;
		*(Out+i) = x;
	}
	write_sound(Out, sample_rate, length);	
	delete Out;

return 0;
}

Download full source code: DecayingSine.cpp DecayingSine.java

In the above code sample, the intial conditions have been set to mimic the sound of a boomy Roland TR-808 kick drum.

4 - Decaying Sinusoids with Harmonics

What we have done so far is good. However, musicians need a lot more than just single decaying sinusoids to make music. Real world sounds contain many frequencies. Because of this, the logical next step is for us to use a summation of decaying sinusoids. This is shown below in Eq. 25 and Fig. 7.

(25)

Expanding Eq. 25, we get

(26)

Each term in the above series represents a separate frequency and can be treated as an independent damped oscillator. This is illustrated below as a group of spring-mass systems for a harmonic sound.

If our sound is harmonic, we can find for the n^th frequency simply by using the following equation.

(27)

In all cases, b and k can be set for the n^th frequency by using the following sequences.

		(28)
		(29)

In the above, n is an integer given by n = 1, 2, 3, 4,....

Using Eq. 27, Eq. 28, and Eq. 29, we can now modify our C++ code and synthesize the fundamental frequency along with harmonics.

int main(){	
	double t60=2.0, w, b, k, x=0.0, a, v=1800.0, tau, *Out, f_note=400.0;
	double t60_2=2.0, w2, b2, k2, x2=0.0, a2, v2=1800.0;
	long length, sample_rate=44100, i;
	
	length = (long)( t60*(double)sample_rate );
	Out = new double[length];
	tau = 1.0/(double)sample_rate;
	w = 2.0*PI*f_note;
	b = -2.0*log(0.001)/t60;
	k = w*w + b*b/4.0;
	w2 = 2.0*w;
	b2 = -2.0*log(0.001)/t60_2;
	k2 = w2*w2 + b2*b2/4.0;
	
	for( i = 0; i < length; i++){
		x += v*tau;
		a = -k*x -b*v;
		v += a*tau;
		x2 += v2 * tau;
		a2 = -k2 * x2 - b2 * v2;
		v2 += a2 * tau;
		*(Out+i) = x+x2;
	}
	write_sound(Out, sample_rate, length);	
	delete Out;

return 0;
}

Download full source code: Harmonics.cpp Harmonics.java

5 - Setting the Initial Velocity

To make the synthesis cleaner and more elegant, we should set our initial velocity so that the highest maximum of our decaying sine wave has the value we want it to have. Our equations are apt to get quite messy if we try to do this for the summation of all our harmonics. We can, however, find the peak of each individual harmonic and then find highest possible maximum for the summation. This is done by solving a fairly basic boundary value problem. To set the value of the highest maximum for the n^th harmonic we use the following equation.

(30)

_n is the angular frequency of the n^th given by Eq. 27 and A_n is the constant that sets the maximum peak of x_n. The value of A_n is found by using the following equation.

(31)

The variable h_n is equal to the highest maximum of the n^th harmonic. The variable t_n(max) is the time taken for the oscillator to go from x_n(0) = 0 to x_n(t) = h_n. It is found by using Eq. 32.

(32)

To find the highest possible maximum for the summation of harmonics we merely sum h_n from n to N as shown below.

(33)

This is far from an exact number, but it should at least get us in the ball park.

If we wanted to use the code from Fig. 6 but have the wave peak at x_n(t_n(max)) = h_n we could use the following code as a substitute.

int main(){	
	double t60=2.0, w, b, k, x=0.0, a, v, tau, f_note=100.0, *Out, h, t, A;
	long length, sample_rate=44100, i;
	
	length = (long)( t60*(double)sample_rate );
	Out = new double[length];
	tau = 1.0/(double)sample_rate;
	w = 2.0*PI*f_note;
	b = -2.0*log(0.001)/t60;
	k = w*w + b*b/4.0;
	t = atan(2*w/b) / w;
	h = 0.5;  
	A = h * pow( E, b*t/2 ) / sin(w*t);
	v = A * w;
	
	for( i = 0; i < length; i++){
		x += v*tau;
		a = -k*x -b*v;
		v += a*tau;
		*(Out+i) = x;
	}
	write_sound(Out, sample_rate, length);	
	delete Out;

return 0;
}

Download full source code: Velocity.cpp Velocity.java

Because we are solving the problem numerically we will have a little bit of round off error. Our wave will not peak exactly at h. At low frequencies our error in the peak is very low. At frequencies below 150 Hz it can be as low as 00.001%. As one approaches the Nyquist frequency, however, our output x_i will be overwhelmed by error.

6 - Reducing Computational Error

All numerical methods suffer from error. There are several ways to quantify the error produced by this synthesis method. One is the peak magnitude of the wave. It can be too high or two low. Another is that the pitch is slightly out of tune. The pitch error can occur not only in the fundamental, but also in the harmonics. This can can cause the harmonics to drift out of phase with one another. Still another error is that the wave form does not decay exactly as we want it to. These problems can vary in severity. All of them tend to worsen at higher frequencies. Below is an image of a 1000 Hz wave where the third harmonic is drifting out of phase with the fundamental.

Figure 10a.   Plot of a wave with the second and third harmonic drifting out of phase with the fundamental frequency.	Figure 10b.   Plot of a wave where one step has been taken on the fundamental frequency, two on the second harmonic, and three on the third harmonic.

The phase drift can be fixed easily and effectively. The drift occurs because there is error in the frequency. This error occurs at a rate that is proportional to the frequency we are trying to synthesize. The frequency is dependent only upon k and b so at any point in time, the frequency error in a given harmonic should be roughly the same. Because of this, we can synchronize the error in the second and fundamental harmonic simply by making twice as many time steps This means that the error in the second harmonic occurs at a rate twice that of the fundamental. Likewise, the error in the third harmonic occurs at a rate three times that of the fundamental. If we make our the time steps in our harmonics smaller so that

(34)

The frequency error in each harmonic, will be the same as the frequency error in the fundamental. Substituting Eq. 27 into Eq. 34 we find,

(35)

To implement this, we make two loops for the second harmonic, and three for the third harmonic every time we make one calculation for the fundamental harmonic. We will be throwing away a lot of numbers, but our harmonics will line up very nicely. This is shown below in Fig. 11.

This technique can be generalized for inharmonic sounds by making the same number of iterations per period for each frequency, and then interpolating the data to obtain the chosen sample rate.

int main(){	
	double t60=1.0, w, b, k, x=0.0, a, v=3500.0, tau, *Out, f_note=900.0;
	double t60_2=1.0, w2, b2, k2, x2=0.0, a2, v2=3500.0, tau2;
	double t60_3=1.0, w3, b3, k3, x3=0.0, a3, v3=3500.0, tau3;
	long length, sample_rate=44100, i, i2, i3;
	
	length = (long)( t60*(double)sample_rate );
	Out = new double[length];
	tau = 1.0/(double)sample_rate;
	w = 2.0*PI*f_note;
	b = -2.0*log(0.001)/t60;
	k = w*w + b*b/4.0;
	tau2 = tau/2.0;
	w2 = 2.0*w;
	b2 = -2.0*log(0.001)/t60_2;
	k2 = w2*w2 + b2*b2/4.0;
	tau3 = tau/3.0;
	w3 = 3.0*w;
	b3 = -2.0*log(0.001)/t60_3;
	k3 = w3*w3 + b3*b3/4.0;
	
	for( i = 0; i < length; i++){
		x += v*tau;
		a = -k*x -b*v;
		v += a*tau;		
		for(i2=0;i2<2;i2++){
			x2 += v2*tau2;
			a2 = -k2*x2 -b2*v2;
			v2 += a2*tau2;
		}				
		for(i3=0;i3<3;i3++){
			x3 += v3*tau3;
			a3 = -k3*x3 -b3*v3;
			v3 += a3*tau3;
		}
		*(Out+i) = x+x2+x3;
	}
	write_sound(Out, sample_rate, length);	
	delete Out;

return 0;
}

Download full source code: PhaseAlign.cpp PhaseAlign.java

We still have a frequency error in our fundamental frequency. The easiest way to fix this is to make smaller and take several loops over all the code in the for-loop in Fig. 11. is now given by

(36)

Where u is an integer constant giving the number of extra global loops for each output sample. Before we can determine what u should be set to, we must first state what we will consider to be an acceptable amount of error. In music there is something know as the just noticeable difference. This term refers to the smallest perceivable pitch interval. It is often stated to be around 5 cents (there are 100 cents is each half-step). A half-step, also known as a semitone, is the pitch difference between two successive keys on a piano. The frequency relation for each successive key on a piano is

(37)

Similarily, the frequency relation for each cent is given by

(38)

As a general rule of thumb, we can satisfy a 5 cent error limit by using an integer value of u that obeys the following relation.

(39)

Our code from Fig. 11 becomes

int main(){	
	double t60=1.0, w, b, k, x=0.0, a, v=3500.0, tau, *Out, f_note=900.0, u;
	double t60_2=1.0, w2, b2, k2, x2=0.0, a2, v2=3500.0, tau2;
	double t60_3=1.0, w3, b3, k3, x3=0.0, a3, v3=3500.0, tau3;
	long length, sample_rate=44100, i, i2, i3, j;
	
	length = (long)( t60*(double)sample_rate );
	Out = new double[length];
	u = (int)ceil(24.0*(double)f_note/(double)sample_rate);
	tau = 1.0/(double)(u*sample_rate);
	w = 2.0*PI*f_note;
	b = -2.0*log(0.001)/t60;
	k = w*w + b*b/4.0;
	tau2 = tau/2.0;
	w2 = 2.0*w;
	b2 = -2.0*log(0.001)/t60_2;
	k2 = w2*w2 + b2*b2/4.0;
	tau3 = tau/3.0;
	w3 = 3.0*w;
	b3 = -2.0*log(0.001)/t60_3;
	k3 = w3*w3 + b3*b3/4.0;
	
	for( i = 0; i < length; i++){		
		for(j=0;j < u; j++){
			x += v*tau;
			a = -k*x -b*v;
			v += a*tau;		
			for(i2=0;i2<2;i2++){
				x2 += v2*tau2;
				a2 = -k2*x2 -b2*v2;
				v2 += a2*tau2;
			}			
			for(i3=0;i3<3;i3++){
				x3 += v3*tau3;
				a3 = -k3*x3 -b3*v3;
				v3 += a3*tau3;
			}	
		}
		*(Out+i) = x+x2+x3;
	}
	write_sound(Out, sample_rate, length);	
	delete Out;

return 0;
}

Download full source code: InTune.cpp InTune.java

We can make another improvement by using the leap-frog method instead of the Euler-Cromer. The leap-frog method is similar to the Euler-Cromer, but unlike the Euler-Cromer, it does not use the velocity at the current position as an input. The leap-frog method uses the velocity that occurs at the position in the middle between the current and next position. This increases accuracy because the velocity is closer to the average velocity over the entire step. To do this, all we have to do is use a different starting velocity. Specifically, we use a starting velocity that is the velocity necessary at time t = /2 for our mass to reach our desired amplitude. Because this calculation is done outside the for-loop, its implementation is also cheap. The new equation for our initial velocity is given by

(40)

We can describe this scheme with the following equations of motion.

(41a)

(41b)

(41c)

Our code for initializing the velocity will now look like the following

int main(){	
	double t60=1.0, w, b, k, x=0.0, a, v, h, A, t, tau, f_note=150.0, u, *Out;
	double t60_2=1.0, w2, b2, k2, x2=0.0, a2, v2, h2, A2, t2, tau2;
	double t60_3=1.0, w3, b3, k3, x3=0.0, a3, v3, h3, A3, t3, tau3;
	long length, sample_rate=44100, i, i2, i3, j;
	
	length = (long)( t60*(double)sample_rate );
	Out = new double[length];
	u = (int)ceil(24.0*(double)f_note/(double)sample_rate);
	
	tau = 1.0/(double)(u*sample_rate);
	h = 0.3;
	w = 2.0*PI*f_note;
	b = -2.0*log(0.001)/t60;
	k = w*w + b*b/4.0;
	t = atan(2.0*w/b) / (w);
	A = h * pow( E, b*t/2 ) / sin(w*t);	
	v = A * pow(E, -b*tau/4.0) * (cos(w*tau/2.0)*w-sin(w*tau/2.0)*b/2.0);	
	
	tau2 = tau/2.0;
	h2 = h;
	w2 = 2.0*w;
	t60_2 = t60;
	b2 = -2.0*log(0.001)/t60_2;
	k2 = w2*w2 + b2*b2/4.0;
	t2 = atan(2.0*w2/b2) / (w2);
	A2 = h2 * pow( E, b2*t2/2.0 ) / sin(w2*t2);
	v2 = A2 * pow(E, -b2*tau2/4.0) * (cos(w2*tau2/2.0)*w2-sin(w2*tau2/2.0)*b2/2.0);	
	
	tau3 = tau/3.0;
	h3 = h;
	w3 = 3.0*w;
	t60_3 = t60;
	b3 = -2.0*log(0.001)/t60_3;
	k3 = w3*w3 + b3*b3/4.0;
	t3 = atan(2.0*w3/b3) / (w3);
	A3 = h3 * pow( E, b3*t3/2.0 ) / sin(w3*t3);
	v3 = A3 * pow(E, -b3*tau3/4.0) * (cos(w3*tau3/2.0)*w3-sin(w3*tau3/2.0)*b3/2.0);
	
	for( i = 0; i < length; i++){		
		for(j=0;j < u; j++){
			x += v*tau;
			a = -k*x -b*v;
			v += a*tau;		
			for(i2=0;i2<2;i2++){
				x2 += v2*tau2;
				a2 = -k2*x2 -b2*v2;
				v2 += a2*tau2;
			}			
			for(i3=0;i3<3;i3++){
				x3 += v3*tau3;
				a3 = -k3*x3 -b3*v3;
				v3 += a3*tau3;
			}	
		}
		*(Out+i) = x+x2+x3;
	}
	write_sound(Out, sample_rate, length);	
	delete Out;

return 0;
}

Download full source code: LeapFrog.cpp LeapFrog.java

The leap-frog method has no effect on our frequency. What it does do is bring us closer to our desired value of h. When our variable u from Eq. 39 is not set properly, say u = 1 and when our frequency f_note = 8,000 Hz, the improvement is very dramatic. When u is set properly, the improvement, though very clear, is nearly negligible in practice.

7 - Conclusions

The Euler-Cromer method works surprisingly well for synthesis. It has accurate pitch and decay at low frequencies. In many circumstances it is also more efficient than using sin together with pow functions. There is a slight problem with the phase of the harmonics drifting away from one another, but this can easily be fixed by using a few extra loops.

The Euler-Cromer method becomes more problematic at higher frequencies. Pitch becomes less accurate, and the waveform has losses in terms of its periodic appearance. The decay time, however, remains accurate.

The fourth order Runge-Kutta produces musical signals with a pitch more accurate than the Euler-Cromer. It's decay time, on the other hand, is far less accurate. To get an acceptable decay time at high frequencies with the fourth order Runge-Kutta, one must perform extra iterations. Because the Runge-Kutta is computationally more expensive, the Euler-Cromer falls below our error limit sooner.

An improvement in the maximum peak value can be obtained by modifying the Euler-Cromer method so that it becomes what is known as the leapfrog method. This improvement can be very pronounced when one is willing to except a larger pitch error.

While the usefulness of this synthesis method might not be readily apparent, it does shine light on other possibilities. It is similar to filter based synthesis, in fact it is arguable that this method of synthesis can be implemented as a digital audio filter. Rather than starting with an initial velocity, one could take an input signal and feed it, sample by sample, into the velocity or the acceleration. In this sense, the spring mass simulation functions as a digital audio filter.

Spring mass systems are common place in physics. Often springs and masses will be attached end to end in a chain. These chained systems can be modeled just as the single spring mass system and these chains can in turn be implemented as filters. We can take the Laplace transform of the differential equation given in Eq. 21 and from there find a transfer function. If we then apply the bi-linear transform, we cab construct an ordinary IIR filter.

[1]	Richard G. Baldwin Java Sound, Creating, Playing, and Saving Synthetic Sounds http://www.developer.com/java/other/article.php/2226701, 2003
[1]	Behringer Racktuner BTR2000 User's Manual Behringer Spezielle Studiotechnik GmbH, Germany, 2005.
[2]	Boas, Mary L. Mathematical Methods in the Physical Sciences. Wiley, 2nd edition, 1983.
[3]	Perry Cook Real Sound Synthesis for Interactive Applications. A K Peters, Wellesley, Massachusetts, 2nd edition, 2002-2004.
[4]	Coppers, Frey, Kinsler, Sanders. Fundamentals of Acoustics. Wiley, 3rd edition, 1982.
[5]	Ehrlich, Robert Euler's Method for Coupled Differential Equations; RLC Circuits. http://35.9.69.219/home/modules/pdf_modules/m351.pdf, 3/11/2002
[6]	Gagen, Alex; Larson, Sean Coupled Oscillators online.redwoods.cc.ca.us/instruct/darnold/DEProj/Sp00/SeanAlex/finalpaper.pdf, 2000
[7]	Alanjandro Garcia Numerical Methods for Physics. Prentice Hall, 2nd edition, 1994-2000.
[8]	Halliday, Resnick, and Walker Fundametals of Physics. Addison Wesley, 5th edition, 1999.
[9]	Hyper-Physics http://hyperphysics.phy-astr.gsu.edu/hbase/oscda.html http://hyperphysics.phy-astr.gsu.edu/hbase/music/cents.html
[10]	Jack, Hugh Laplace Transforms claymore.engineer.gvsu.edu/~jackh/books/model/chapters/laplace.pdf, 1996-2003
[11]	Kent Nagle and Edward Saff Fundametals of Differential Equations and Boundry Value Problems. John Wiley & Sons Inc., 2nd edition, 1997.
[12]	Korg Inc. DTR-2000/DTR-1000 Digital Tuner Owner's Manual Korg Inc, Japan, 2002.
[13]	~pat Coupled Oscillations. http://www.scar.utoronto.ca/~pat/fun/NEWT1D/PDF/COUPLED.PDF, 1999
[14]	Sanjit Mitra. Digital Signal Processing, A Computer-Based Approach. McGraw-Hill, 2nd edition, 2001.
[15]	Robert Sedgewick and Kevin Wayne. Introduction to Programming in Java http://www.cs.princeton.edu/introcs/home/, 2007
[15]	Smith, Julius O. Physical Audio Signal Processing: for Virtual Musical Instruments and Digital Audio Effects, August 2006 Edition Center for Computer Research in Music and Acoustics (CCRMA), Stanford University http://ccrma.stanford.edu/~jos/pasp06/, Aug. 2006
[16]	Keith Symon Mechanics. Addison Wesley, 3nd edition, 1971.
[17]	Wikipedia Article: Just noticeable difference http://en.wikipedia.org/wiki/Just_noticeable_difference
[18]	Yang, Dae Ryook Laplace Transform and Transfer Function www.cheric.org/ippage/e/ipdata/2001/16/file/L5-LTxTF.pdf, 2001
[19]	Young, Hugh D. Coupled Oscillators and Normal Modes www.andrew.cmu.edu/course/33-231/coup-osc.pdf, 2003