Creating complex waveforms through Fourier / Additive Synthesis

This post belongs to a sub-series on complex waveforms, in my bigger series on developing a digital Synthesizer in C++. This post will be about creating complex waveforms through Fourier Synthesis, also known as Additive Synthesis. As introduced in the parent post, this method uses the summation of a finite series of sine waves. Please read about the basic idea behind this type of synthesis in the parent post if you haven't yet, as I will really focus on concrete examples and applications now and only explain what partials must be added with what properties.

EDIT: I have published a revised version of the code, however I do all the explaining here, so I recommend reading this post first.  

The code used and discussed here presupposes knowledge of how a sine wave is created in code. I will introduce a function for additively creating a square wave in the first part of this post, to which we will then add functionality step by step for each waveform so that at the end of this article you will have a fully generic function that can create any waveform imaginable with a variable number of partials.

Waveform I - The square wave

In terms of Fourier Synthesis, square waves are created by summing all odd partials (3rd, 5th, 7th etc.) of the fundamental pitch with the respective amplitude being the reciprocal (1 / x) of the partial number (1/3, 1/5, 1/7 etc.).


Here are some pictures of square waves with 2, 4, 8, 16, 32 and finally 64 partials:

As you can see, the more partials, the closer to perfection is the waveform. Here is a function to create a sample buffer for such waveforms:

double * genWave(double fundFreq, int maxPartial, int bufferLen = 44100)
{
    // the sample buffer
    double * buffer    = new double [bufferLen];
    
    // increment by one because for our calculations, the
    // fundamental pitch is a partial too
    maxPartial++;
    
    // holds the current phase for all partials
    double * phase     = new double [maxPartial];
    
    // their phase increments
    double * phaseIncr = new double [maxPartial];
    
    // their amplitudes
    double * amp       = new double [maxPartial];
    
    // the current partial, starts at one (fundamental)
    int currPart = 1;
    
    // the increment of the fundamental pitch
    double fundIncr = (twoPI/bufferLen) * fundFreq;
    
    for (int p = 0; p < maxPartial; p++)
    {
        // partials usually start at a phase of 0,
        // however you could offset some partials if
        // you want to
        phase[p]     = 0.0;
        phaseIncr[p] = fundIncr * currPart;
        amp[p]       = 1.0 / currPart;
        
        // the step by which to increment the partial number
        // because we want all odd partials, we add two to
        // get to the next odd partial e.g. 3 -> 5 -> 7
        currPart += 2;
    }
    
    // calculate each sample
    for (int sample = 0; sample < bufferLen; sample++)
    {
        double value = 0.0;
        
        // for each sample, get all partial values
        // and sum them together
        
        for (int p = 0; p < maxPartial; p++)
        {
            // add partial value and multiply by
            // respective amplitude
            value += sin(phase[p]) * amp[p];
            
            // increment respective phaes
            phase[p] += phaseIncr[p];
            
            // check if we have to reset phasor
            if (phase[p] >= twoPI)
                phase[p] -= twoPI;
        }
        
        buffer[sample] = value;
        
    }
    
    delete [] phase;
    delete [] phaseIncr;
    delete [] amp;
    
    return buffer;
}

The algorithm above fills dynamic arrays of size maxPartial + 1 with initial phase, phase increment and amplitude according to the description I gave above on how square waves are formed. We must increment maxPartial because, if you recall, the fundamental pitch is sort of seen as the first partial for this algorithm, so that it has a frequency of 1 * fundamental frequency and 1 as it's amplitude. You could also just fill the first entries for the arrays with the fundamental frequency’s items intially, so you set phaseIncr[0] to fundIncr and amp[0] to 1 and set the part variable to 2, however I think this method is clearer. We increment the partial number for each array entry by 2, because we want all odd partials: 1 + 2 = 3, 3 + 2 = 5 etc. 

When calculating the sample buffer, we calculate the value for each partial like we did for a single sine wave and add them all together. 

If you call this function and write the sample buffer to a WAV file, it might interest you to see how the number of partials affects the shape of the waveform. The more partials, the more it will look like a square wave. Calling this function with maxPartial set to 1 will yield a waveform where you can really distinguish how this partial affects the fundamental one. For the waveform the nearest to perfection I recommend 64 partials, as a higher partial number is really unnecessary and yields little improvement.

In my last post I mentioned something about a “hack” to undo the so called “Gibbs phenomenon” (these horns at the ends of the additively created waveforms), namely the Lanczos sigma factor. It is quite an interesting little factor, so I devoted another blog post to it. I recommend that you read about the theory in that very short post so that you have a better understanding of what I will show you now.

The sigma factor is defined as:

σ = sin (x) / x 

x being:

x = nπ / M

Where n is the current partial number and M the number of total partials. This little factor σ we will then multiply each partial's amplitude by. 

We can refactor the second equation so that we can calculate π / M in advance. We will then multiply this constant value by the partial number in each round of the for loop. Firstly, the constant part, declared before the first for loop in the genWave function:

double sigmaConst = PI / maxPartial;

Inside the for loop, we then always need to calculate the partial's 'x'

double sigma = sigmaConst * currPart;

And then finally calculate the actual σ factor and multiply the current amplitude by it:

amp[p] *= sin(sigma) / sigma

This gives us the following changed for loop in the genWave function:

// constant part of the sigma factor
double sigmaConst = PI / maxPartial;

for (int p = 0; p < maxPartial; p++)
{
    // partials usually start at a phase of 0,
    // however you could offset some partials if
    // you want to
    phase[p]     = 0.0;
    phaseIncr[p] = fundIncr * currPart;
    amp[p]       = 1.0 / currPart;
    
    // variable part of the sigma factor
    double sigma = sigmaConst * currPart;
    
    // multiply the amplitude by it
    amp[p] *= sin(sigma) / sigma;
    
    // the step by which to increment the partial number
    // because we want all odd partials, we add two to
    // get to the next odd partial e.g. 3 -> 5 -> 7
    currPart += 2;
}

If you compile the new code and store the buffer in a wave file, you will notice how amazingly straight the waveform is. However, if you take a look at what I talk about towards the end of my post about the Lanczos sigma factor, you will see why the Lanczos sigma factor should only be used for waveforms you want really near to perfection, as they really straighten out every waveform, even those of which you want the imperfect, curvy shape. Therefore, we should add an if clause in the for loop that computes and adds the sigma factor only if the total partial number is higher than a certain number. As I see a waveform with 64 partials as the waveform I want to be nearest to perfection, this number at which we start using the sigma factor is 64. 

The part where we add the sigma factor inside the for loop should therefore look like this now:

if (maxPartial >= 64)
{
    // variable part of the sigma factor
    double sigma = sigmaConst * currPart;
    
    // multiply the amplitude by it
    amp[p] *= sin(sigma) / sigma;
}     

  Waveform II - The Sawtooth wave

The Sawtooth waveform is arguably the most straight-forward waveform to create additively. It differs in the algorithm used for square waves only in the fact that instead of using only all odd partials, we use every partial. Therefore, we will not increment the partial count by 2, but simply by 1, so the 1st, the 2nd, the 3rd partial and so on. The amplitude is still the reciprocal of the partial number. So 1/2, then 1/3, then 1/4 etc.

Here are some pictures of sawtooth waves with 2, 4, 8, 16, 32 and finally 64 partials (last one with σ factor):

We will now add the first element that makes the previous function more generic. Instead of statically incrementing the partial number by 2, we will use a step parameter.

Add a step parameter to the function header:

From:

double * genWave(double fundFreq, int maxPartial, int bufferLen = 44100)

To:

double * genWave(double fundFreq, int maxPartial, int step, int bufferLen = 44100)

And change where we increment the step in the first for loop:

From:

currPart += 2;

To:

currPart += step;

Now you can call the function with step set to 1 and you will get a nice sawtooth waveform. Setting it to 2 gives a square wave again.

Waveform III - Ramp wave

Ramp waves differ from sawtooth waves in the way that they ascend from low to high instead of descending from high to low as do saw waves. The only change in code that we have to make, therefore, is to make all amplitudes negative. Technically. Our aim is it to end up with ONE FUNCTION TO RULE THEM ALL! Wow, what was that. I meant one function to generate them all, of course. Therefore, we will also introduce an ID parameter, with which we can add some extra features to our genWave function. 

Ramp waves with 2, 4, 8, 16, 32 and 64 partials (last one with σ factor):

As mentioned, we must set all amplitude values to negative, the way we calculate it stays the same as for the other two waveforms. Because we don't really have any means of distinguishing a ramp wave from a sawtooth wave purely from the parameters, we must add another parameter to the function.

Change the function header from:

double * genWave(double fundFreq, int maxPartial, int step, int bufferLen = 44100)

To:

double * genWave(double fundFreq, int maxPartial, int step,
                   int ID, int bufferLen = 44100)

Now that we have a way of distinguishing between waveforms, we can say that the square wave will have the ID 0, the sawtooth wave the ID 1 and the ramp wave then the ID 2. With this, we can make
the this minor add-on to our function. Change:

amp[p] = 1.0 / currPart;

To:

// Ramp
if (ID == 2) amp[p] = -(1.0 / currPart);

else amp[p] = 1.0 / currPart;

in the first for loop, in which we set the partial parameters. If you now call the function for example like this:

double * buffer = genWave(220, 64,1,2);

we will get the sample buffer for a ramp wave of 220 Hz frequency and 64 partials.

Waveform IV - The Triangle wave

Triangle waves are a bit different to the ones I showed you how to make until now, mainly in respect to it's amplitude. We used the formula 1 / n to compute the amplitude for all our waveforms until now and only the ramp waved differed in the fact that the amplitude was all negative.

For one, the amplitude for triangle waves is no longer 1 / n, but 1 / n^2. Moreover, the amplitude sign changes for every partial, so we must check the sign of the previous partial's amplitude on each iteration and use the opposite sign for the current partial. The partials we use are the same as for the square wave, namely all odd partials.

Here a picture of triangle waves with 2, 4, 8, 16, 32 and finally 64 partials (last one with the sigma factor). Note that they do not change a lot and that even two partials already give a neat triangle wave:

Since we introduced the ID parameter to our function for the ramp wave, we can have the triangle wave have the ID 3. As mentioned, we must for one change the way the amplitude is calculated and also somehow change the amplitude's sign on every iteration. This leads us to the following if clause:

// Triangle
else if (ID == 3)
{
    // amplitude = 1 / n ^ 2
    double val = 1.0/(currPart*currPart);
    
    // Change sign
    // First make sure we have a previous
    // element to check
    if (p > 0 && amp[p-1] > 0)
    {
        amp[p] = -val;
        
    } else {
        
        amp[p] = val;
    }
    
}

The first part of this if clause declares and initialises a variable val and computes the formula 1/ n ^ 2 in accordance to the current partial number n. Then, we check if the previous partial, so the previous entry in the amp array, was positive or negative. If it was positive, we set the amplitude for the current partial to the negative val, else to the positive val. We need to do check first if p is greater 0, since amp[-1] is invalid. We combine this with the other conditions so that the fundamental frequency has a positive amplitude.

Outro

The adding of the functionality for the triangle wave was our last step in creating our super awesome function. We now have a function that can create a variable waveform with a variable number of partials. I will now post the entire function, however I changed a few things that should be noted. I removed the step parameter, since, as we have the ID parameter, we know what waveforms need what steps. Sine waves now have the ID 0, square waves the ID 1, sawtooth is 2, ramp 3 and the triangle wave has the ID 4. We now just have to check for the ID and set the step accordingly.

Voilà:

double * genWave(int ID, double fundFreq, int maxPartial, int bufferLen = 44100)
{
    // the sample buffer
    double * buffer    = new double [bufferLen];
    
    int step;
    
    // SINE
    if (ID == 0) step = 0;
    
    // SAW and RAMP
    else if (ID == 2 || ID == 3) step = 1;
    
    // TRIANGLE and SQUARE
    else step = 2;
    
    // increment by one because for our calculations, the
    // fundamental pitch is a partial too
    maxPartial++;
    
    // holds the current phase for all partials
    double * phase     = new double [maxPartial];
    
    // their phase increments
    double * phaseIncr = new double [maxPartial];
    
    // their amplitudes
    double * amp       = new double [maxPartial];
    
    // the current partial, starts at one (fundamental)
    int currPart = 1;
    
    // the increment of the fundamental pitch
    double fundIncr = (twoPI/bufferLen) * fundFreq;
    
    // constant part of the sigma factor
    double sigmaConst = PI / maxPartial;

    for (int p = 0; p < maxPartial; p++)
    {
        // partials usually start at a phase of 0,
        // however you could offset some partials if
        // you want to
        phase[p]     = 0.0;
        phaseIncr[p] = fundIncr * currPart;
        
        // Ramp
        if (ID == 3) amp[p] = -(1.0 / currPart);
        
        // Triangle
        else if (ID == 4)
        {
            // amplitude = 1 / n ^ 2
            double val = 1.0/(currPart*currPart);
            
            // Change sign
            // First make sure we have a previous
            // element to check
            if (p > 0 && amp[p-1] > 0)
            {
                amp[p] = -val;
                
            } else {
                
                amp[p] = val;
            }
            
        }
        
        else amp[p] = 1.0 / currPart;
        
        if (maxPartial >= 64)
        {
            // variable part of the sigma factor
            double sigma = sigmaConst * currPart;
            
            // multiply the amplitude by it
            amp[p] *= sin(sigma) / sigma;
        }

        // the step by which to increment the partial number
        // because we want all odd partials, we add two to
        // get to the next odd partial e.g. 3 -> 5 -> 7
        currPart += step;
    }
    
    // calculate each sample
    for (int sample = 0; sample < bufferLen; sample++)
    {
        double value = 0.0;
        
        // for each sample, get all partial values
        // and sum them together
        
        for (int p = 0; p < maxPartial; p++)
        {
            // add partial value and multiply by
            // respective amplitude
            value += (sin(phase[p])) * amp[p];
            
            // increment respective phaes
            phase[p] += phaseIncr[p];
            
            // check if we have to reset phasor
            if (phase[p] >= twoPI)
                phase[p] -= twoPI;
        }
        
        buffer[sample] = value;
        
    }
    
    delete [] phase;
    delete [] phaseIncr;
    delete [] amp;
    
    return buffer;
}

Feel free to comment if you have any questions or suggestions. You should also check out my post on creating complex waveforms through direct calculations!