Recursive Synthesis of Triangle, Sawtooth and Square waves

First, let's look at the common aspects of digital synthesis. It will help us deduce some formulae. Surely, you know all about the concepts of digital sound, but in spite of this I recommend you not to skip it. So, we have a ceaseless analog signal - perfect and smooth (for example a sine wave):

Fig. 1

This signal comes from "real life". We should transfer it to the
computer world - the world of numbers. How is this done? There is a parameter called
*discretization frequency*. It means that we divide each second into a big
amount of parts or *samples*. The discretization frequency is the number of
such parts per second. The value of each sample is simply the voltage of the
original signal taken from the corresponding time period and converted by the soundcard's
ADC (Analog-to-Digital-Converter) to signed computer words - this is 16-bit recording.
That is, we replace the original ceaseless signal with a discrete signal and therefore
the higher discretization frequency the better. So, we get a discrete signal:

Fig. 2

Terrible, isn't it? You can say: "We have never seen such distortions in sound editing programs like Sound Forge or Cool Edit!" Right, because these programs use interpolation to draw the wave on the screen. But this interpolation is only graphical, the actual signal is always discrete as in fig. 2.

Well, we're going on. Now we can come close to the topic of this
article. There's a question: can we produce for example a sine wave with a frequency of
6,000 Hz if the discretization frequency will be 11,025 Hz? Answer: no, we can't.
A well known Nyquist theorem claims: "If we have a discrete signal with
discretization frequency ** frq **then the highest frequency we can
store in this signal is

Now look at the following: (sorry for small pictures)

Fig. 3

There are some symbols: *Xn* - these are the values of time when the
corresponding sample values are taken. Let's define a function *f(x)* - it
determines the waveform of our sound. In our case *f(x)*=sin *x.*
Obviously, the samples in the digital signal will be *f(X0)*,
*f(X1)*, *f(X2)*, ..., *f(Xn)*

And now we define a value called *step* - this is the interval
of time between two neighboring samples. If you think a bit, you'll notice that
*step*=1/*discretization_frequency*.

I think that you see that X1=X0+*step*,
X2=X1+*step*, ..., Xn=Xn-1+*step*

This is a recurrence relation, isn't it? But this is not enough
to say about recursive synthesis. However,
IF WE FIND THAT WE CAN **DERIVE** *f(Xn)* FROM SOME AMOUNT
OF **PREVIOUS** SAMPLES (PREFERABLY FROM *f(Xn -1)*)
THEN WE

*f(Xn)*=*f(Xn-1)*+*A*

(*) or

*f(Xn)*=*f(Xn-1)***A*

(**) or

a combination of this (***), where *A* is 1. constant (the
best case) or 2. some equation with changing variables. The second case is less preferable
than first 'cos we need to calculate this *A* at every step of calculating
*f(x)*.

You might say "This is crap, I should perform multiplication or addition to produce every sample! Using the Look-Up Tables (LUT) is much easier! I should only read a number from an array, not having to calculate it!" But I should tell you: ATTENTION! LOOK-UP TABLES ARE INAPPLICABLE IN SOUND SYNTHESIS!!! You simply receive wrong waves with them! If you are interested why, you can send an e-mail to me and I'll explain this fact to you. Also with LUTs you can't get waves with changing frequencies. Well, you can, but you MUST recalculate LUT for every new frequency! It eats a lot of time. With recursive synthesis changing a wave's frequency is as easy as the fact that you can't live for ever at this Earth.

So, to check the capability of your function to be recursive you should use the following algorithm: (I think that the copyright is mine 'cos I haven't ever seen this)

1. Write down the full formula of your function *f(x)*, for
example, *f(x)*=2x (quite a simple case, isn't it?)

2. Then take the value of your formula not from *x* but
from *x+step*: *f(x+step)*=2(*x+step*)=2*x*+2*step*

3. Try to find the way of deriving *f(x+step)* from *f(x)*
- try to find an equivalent formula which has a form (*) or (**) or (***). In
our case we have already found such a formula: 2*x*+2*step*, because
2*x* is a *f(x)* and 2*step* is a constant - this is the (*)
case. In other formulae, i.e. for *cosine* or *sine*, this step is
a bit harder, even a big bit harder, so you should know some school
and university formulae perfectly. :)

When you have got a recursive formula you can jump and yell, "I've done it!" :) If you have read this article from the beginning you should already have some ideas about how to code this.

If you can't derive such a formula you have two ways:

1. Say "This formula doesn't have a recursive form!"

2. Try to find this form once more. ;-(

*[Ed.: Iliks now knows a more intelligent way of doing this. It will be explained in Hugi 23.]*

That is the whole algorithm. The following text contains just examples of application of this method in order to synthesize common waveforms. Be sure to read it.

Synthesis of some common waveforms

Uff, yeah! The theoretical part has been done! Now let's look at the short index of the following:

1. Synthesis of sine wave

2. Synthesis of sawtooth wave

3. Synthesis of triangle wave

4. *Bonus*. Synthesis of square wave

5. Pulse Width Parameter and how to implement it.

Sine wave

Well, I must give a credit. I saw the final formula
of the following deduction in Franky's article in Hugi 16 but he didn't
prove it at all. He simply said that he deduced it using derivations.
Well, I agree with him that *(cos wx)'=-w*sin wx*, but he simply wrote
this statement without proving that *sin=sin+wcos*. Unlike him, I *proved*
this formula, and there's no need of using derivations at all, more, I have
no idea how to prove it with derivations. Franky, if you are reading this article
now then please contact me and give me your proof of this, if you have it, of
course! So, let's start. We will use my algorithm mentioned above.

We have

where *A* is the amplitude, or volume of our wave (0-32767), and the symbol after
*Pi* is the frequency of our wave, not the discretization frequency.

If you don't know physics, you can wonder from where

appears. Well, there is a concept of *cyclic frequency* in physics, which
shows us how many periods of wave are restricted in the interval of time equal
to *2*pi* seconds - 6, 28 seconds approximately. I won't explain all this
stuff, it will be better to say that if you want to produce a sine wave with the frequency
*frq* then your function should have the following form: *f(x)=sin(2Pi*frq*x)*

Then we should derive *f(x+step)* from *f(x).* Ok,
let's do it. I think it's obvious that

Now we should use the following school formula:

We'll get

Next, we should use more complicated mathematical stuff, the theory of limits. If you don't know what's the hell this means I can't help you - it's quite a big part of math analysis. Simply enter university and you'll know all these things.

For further reasoning we need two limits:

(1)

and

(2).

I can
deduce this limits but I won't do it because it will take too much space
and this is not our main topic. Have you got any idea about how to use these
equations in the previous formula? !Tada-dam! This is simple. Somewhere in the
beginning of this article I wrote that *step*=1/*discretization_frequency.*
The discretization frequency is quite a big number and therefore its reciprocal
value will be a very small number, very close to zero. And therefore if we multiply
any number with this number we will get also a number which is very close to zero.
Thus, we can say that

And only after this fact we can use limits (1) and (2). So, our final formula will be the following:

or

We have deduced the recursive formula. However, we see that we Also have to calculate. Which is why we should also deduce one more formula.

But it is very easy after aforecited reasoning. You know that

We should do the same deductions as for sine. And this is already your job if you wanna see the following proved. I just give the final recursive formula for cosine:

We have got all the necessary things to implement this technique
In the computer. Next I'll give you a procedure in *C* that synthesizes a sine wave,
but now I want to give you some advice. This algorithm has its own limitations,
like all things do. For good fidelity of waveform we should use the calculating
features of the computer very carefully, because in the recursive formula any
following value of the function depends on the previous ones, so if we do inaccurate
calculations we will get poor waveform quality in the end. Thus, just avoid
possible loss of data when converting, dividing, multiplying numbers! Meet:
The Source.

voidsinus(floatfrq) {floattemp=2.0*M_PI*frq/discr_frq;floatcos=16384.0;floatsin=0.0;shortoutput[length_of_sound]; output[0]=sin; for (inti=1; i<=length_of_sound-1; ++i) { sin+=temp*cos; cos-=temp*sin; output[i]=int(0.5+sin); } }

I'm sure you know that *int(0.5+sin)* means *Round(sin)*.
Doubtless, rounding gives us a higher quality waveform. The amplitude of the signal
is the first value of *cos*. The main fact is that you can change the amplitude
and the frequency in the inner loop! However, Frankie wrote about it in Hugi
16. All you have to do for this is to recalculate *temp* according to the
changes.

Sawtooth wave

If you don't know what a sawtooth wave looks like then look at the following:

NB! The waveform comes at the first onset to minimal level after measuring up to the maximal level, not through the zero! Mathematically, there should be no vertical lines in the graph, but sound editors simply draw lines from the previous point to the following. Let's deduce the recursive formula for the sawtooth function.

I think it is obvious that a sawtooth wave is simply a line. Yep,
ingenious things are simple! And all we have to do is to remember the equation of the line.
From school you should know that this is the first-order equation and its
form is *f(x)=kx+b*.

Look at the two points: (0;-A) and (x1;0). We see that they belong to our graph. Thus, their coordinates match our equation. Now we ought to solve system of two equations:

We have found that our sawtooth consists of lines that are set by

equation. Let's find *f(x+step).*

Now about the implementation of this. It is obvious that we should have
a counter which will count the calculated samples and check whether we have got
through the current period of function or not. When we've achieved the minimal value
of the saw, we set the value of the current sample to -A (our amplitude) and
continue calculating. It's quite clear that *x1=T*/2 (half-period of our
function) and *T=1/frq* where *frq* is the frequency of the function. Thus,
*x1=1/(2*frq)*. Next, how to calculate whether we get through a period or
not? If you think, you'll find that discretization_frequency / frequency_of_wave
value is the period of our function but in samples. Catch the point? So, the
source is coming.

voidsawtooth(floatfrq)

{temp1=discr_frq/frq;

float

floattemp=2.0*20000/temp1;

floata;

intcounter=1; a=-20000; //amplitude output[0]=a; for (i=1; i<=length_of_sound-1;++i) { output[i]=output[i-1]+temp+0.5; ++counter; if (counter>int(0.5+temp1)) // !!! { output[i]=a; counter=0; } } }

As in the previous example, as well as in the following ones, you can easily change the frequency and the amplitude, all you have to do is to recalculate the dependent variables.

I repeat once more: you MUST be accurate in your calculations!
You shouldn't say "Who cares!" There is a bottleneck in the line marked
with "!!!". If you don't round the value of *temp1* it will be simply truncated.
That is why you will get a wrong waveform! However, the distortions are not significant
- you won't see them in a sound editor, only in your program which prints the data
of the wav file. But this is so only in this concrete case. Who knows what you'll
get in other applications if you aren't careful! By the way, experiment with
this, check what will be if you don't round *temp1*!

We're going on!

Triangle wave

I won't draw its waveform - I think these things should be obvious - it is simply a composition of two sawtooth waves. I have explained all this stuff in the previous chapter so you should get outside of this all alone. The maximum I can give you is my source.

voidtriangle(floatfrq) {intcounter=int(0.5+discr_frq/(4*frq)); output[0]=0; for (i=1;i<=length_of_sound-1;++i) { if (counter<int(0.5+discr_frq/(2*frq))) { ++counter; output[i]=output[i1]+int(0.5+(4*20000*frq/float(discr_frq))); } else { ++counter; if (counter>float(discr_frq)/frq) { counter=0; output[i]=-20000; } else output[i]=output[i-1]-int(0.5+(4*20000*frq/float(discr_frq))); } } }

Square wave

Well, it's just a bonus because it doesn't fit
the recursive synthesis principles. The square wave is a composition of two horizontal
lines - the first half of the period is a line *y=A* (A-amplitude) and
the second and the last half is a line *y=-A.* Like with the sawtooth wave there are
no null samples at all! Of course if A doesn't equal to zero.

Any further explanation is not necessary. See the program!

voidsquare(floatfrq) {floattemp=discr_frq/frq;floattemp1=temp/2;intcounter=1; output[0]=20000; for(inti=1;i<=length_of_sound-1;++i) { if (counter<=temp1) { output[i]=20000; ++counter; } else { output[i]=-20000; ++counter; if (counter>temp) { counter=0; } } } }

Pulse Width Parameter and how to implement it

The Pulse Width parameter (in the following simply PW) is a fractional value between 0 and 1 which defines "proportions" of our wave. As I know, PW is defined for square and triangle waves. The sawtooth wave is just a triangle wave with PW equal to one. For an exact definition of PW look at the pictures.

For square:

For triangle:

where T is a period of our function.

We have to do only minimal changes in our routines so I just give the complete source to you.

1. For square:

voidsquare(floatfrq) {floatpulsewidth=0.2;floattemp=discr_frq/frq;floattemp1=temp*pulsewidth;intcounter=1; output[0]=20000; for(inti=1;i<=length_of_sound-1;++i) { if (counter<=temp1) { output[i]=20000; ++counter; } else { output[i]=-20000; ++counter; if (counter>temp) { counter=0; } } } }

2. For triangle:

voidtriangle(floatfrq) {floatpulsewidth=0.7;intcounter=int(0.5+(pulsewidth*discr_frq)/(2.0*frq));inttemp=int(0.5+(2*20000*frq/(pulsewidth*discr_frq)));inttemp1=int(0.5+(2*20000*frq/((1-pulsewidth)*discr_frq)));inttemp2=int(0.5+float(discr_frq)/frq);inttemp3=int(0.5+(pulsewidth*discr_frq/frq); output[0]=0; for (inti=1;i<=length_of_sound-1;++i) { if (counter<temp3) { ++counter; output[i]=output[i-1]+temp; } else { ++counter; if (counter>temp2) { counter=0; output[i]=-20000; } else output[i]=output[i-1]-temp1 } } }

If you have no ideas how I've done it then simply reread the explanation for the sawtooth wave, then deduce the formulae for the triangle wave and then prove the formulae for the case with pulsewidth. It's really simple!

Now I have found that sine with pulse width isn't impossible. I think I'll do it, however with more complicated ideas! Sure, not tomorrow but maybe in the next issue of Hugi, who knows? :)

In the end of this article I wanna say the following things:

1. If you have read it from the beginning to the end then I congratulate you! Thanks. I'm sure this information won't be useless for you.

2. If you like this article, I can write several more, especially
for *Hugi*. I can light up the following themes: 1) The complete description
of the *wav* format. 2) *Wow!* Distortion and morphing (like *Smorphi*
does!) 3) Working with *DirectSound* (beginner's course).

How can you force me to write one of this? Simply send me an email in which you write which one you need the most. If I receive more than 10 emails I'll write the next article with the topic which will have the biggest number of votes. I think that Adok and the other members of the Royal Family won't be against this. Come on! You are welcome!

**iliks (Ilja Palopezhencev)
iliks@chat.ru **(preferable)

or