Creating complex waveforms through Fourier / Additive Synthesis - Version 2

In the last post in my "Developing a digital synthesizer in C++" series, I described how Fourier Synthesis, also known as Additive synthesis, could be used to create complex waveforms of any shape or form. In the recent months I have continued to work on my synthesizer and progressed a lot with my knowledge, therefore I have revised and improved lots of my code, including my functions for creating complex waveforms. Because the source code I provided in the last post has lots of problems, I decided to publish my revised code. I will only explain why and how I changed the previous code, so I recommend you have a look at the previous post if you haven't yet, especially if you want to learn about Fourier synthesis in general.

The main problem with the previous code was that it was not generic at all and actually didn't follow one of the basic rules about functions: they should do one thing and one thing only. The function might be good enough if you want to create only a handful of different waveforms, say square, sine and saw, but as soon as we're talking about creating any waveform we want to, using the same function, we'd be getting into problems pretty soon. The solution I found to this problem is quite nice, but it requires a few more initial steps.

In Fourier Synthesis, a partial is generally described by its frequency, its amplitude and optionally its phase offset. Because a partial is an integer multiple of a fundamental frequency and given the fact that wavetables usually store waveforms with a period of exactly 1 Hertz, the frequency is directly related to the number of the partial, so that's all we need.

Now that we know what information we need about each partial, we can create such a struct:

struct Partial
{
    Partial(unsigned short number, double ampl, double phsOffs = 0)
    : num(number), amp(ampl), phaseOffs(phsOffs)
    { }
    
    const unsigned short num;
    double amp;
        
    double phaseOffs;
};

Each partial has a number, which stays constant, an amplitude that we might want to change and also a variable phase offset, in radians, which we'll usually just set to 0 (no offset).

If you recall from my introductory post on complex waveforms and also from my last post on Fourier synthesis, every waveform has a distinct "algorithm" which determines how we get from a bunch of sine waves to, say, a square wave.

Speaking of square waves, they are formed by adding all odd partials (the 1st, 3rd, 5th etc.). To avoid amplitude overflow, we must also set each partial's amplitude to the inverse of their respective partial number, so the 3rd partial has 1/3rd the maximum amplitude of the fundamental sine wave for example. Usually the fundamental sine wave has an amplitude of 1, so the amplitude of the partials is literally just the inverse of their partial number.

We can therefore have a vector of Partials for any waveform we want, by just changing the settings of the partials. For a square wave with 64 partials:

std::vector<Partial> vec;
    
for (int i = 1; i < 128; i += 2)
{
    vec.push_back(Partial(i, 1.0/i));
}

Simple! We construct the first 64 odd-numbered partials with their number and their amplitude.

A saw wave would be:

for (int i = 1; i <= 64; ++i)
{
    vec.push_back(Partial(i,1.0/i));
}

For a ramp wave, change the amplitude's sign in the above code (-1.0/i).

And a triangle wave (odd partials, alternating sign for the amplitude):

double amp = -1; // So that the first amplitude is positive

for (int i = 1; i < 128; i += 2)
{
    amp = (amp > 0) ? (-1.0/(i*i)) : (1.0 / (i*i));
    
    vec.push_back(Partial(i,amp));
}

Now that we have these vectors for our waveforms, here is the function that creates a wavetable for them:

double * Wavetable::genWave(const partVec& partials,
                            double masterAmp,
                            unsigned int wtLen,
                            bool sigmaAprox)
{
    double * wt = new double [wtLen];
    
    double * amp = new double [partials.size()];
    double * phase = new double [partials.size()];
    double * phaseIncr = new double [partials.size()];
    
    // constant sigma constant part
    double sigmaK = PI / partials.size();
    
    // variable part
    double sigmaV;
   
    // 2*pi / wavetable is the fundamental increment (gives one period
    // of a pure sine wave in the wavetable)
    const double fundIncr = 6.2831853071795865 / wtLen;

    // fill the arrays with the respective partial values
    for (unsigned short p = 0; p < partials.size(); p++)
    {
        phase[p] = partials[p].phaseOffs;
        
        phaseIncr[p] = fundIncr * partials[p].num;
        
        amp[p] = partials[p].amp * masterAmp;
        
        if (sigmaAprox)
        {
            sigmaV = sigmaK * partial.num;
            
            amp[p] *= sin(sigmaV) / sigmaV;
        }
    }
    
    for (unsigned int n = 0; n < wtLen; n++)
    {
        double value = 0.0;
        
        // do additive magic
        for (unsigned short p = 0; p < partials.size(); p++)
        {
            value += sin(phase[p]) * amp[p];
            
            phase[p] += phaseIncr[p];
            
            if (phase[p] >= twoPI)
                phase[p] -= twoPI;
        }
        
        wt[n] = value;
        
    }

    delete [] phase;
    delete [] phaseIncr;
    delete [] amp;
    
    return wt;
    
}

I recommend reading the last post on Fourier synthesis, if you haven't yet, to understand what the function does. Once your done, I'm sure you'll agree with me that this method is a lot nicer. We have gotten rid of our ridiculous ID parameter and can now pass on a vector of partials containing information for any waveform we want. Have fun with it and comment below if you have any questions or suggestions.

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!

Creating complex waveforms through direct calculation

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 direct calculation, so following mathematical formulas as opposed to Fourier / Additive synthesis. 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. A lot of what I use here is based on what Daniel R. Mitchell writes about in his book "BasicSynth", which is highly recommended. 

The code used and discussed here presupposes knowledge of how a sine wave is created in code. 

Waveform I - The square wave

The square wave is undoubtedly the easiest to create through calculation, as it is literally nothing else than 2 values, usually 1 and -1 (normalised amplitude scale), interchanging at certain intervals. 

We have two possibilities now, we can either think in terms of samples or in terms of time. The good thing about time is that you can break time into smaller and smaller pieces, however a sample is a sample and there are no half samples. Thus, thinking in terms of time will create less problems with rounding. 

We must calculate the time for one sample and add it to a count until we reach half the period of the waveform, after which we change the volume from maximum to minimum. When we reach the time for one period, we reset the count:

double * square(double freq, int bufferLen = 44100)
{
    // the sample buffer
    double * buffer = new double [bufferLen];
    
    // time for one sample
    double sampleTime = 1.0 / 44100;
    
    double period = 1.0 / freq;
    
    // the midpoint of the period
    double mid = period / 2;
    
    double curr = 0;
    
    // fill the sample buffer
    for (int n = 0; n < bufferLen; n++)
    {
        // if the current time is less than
        // half the period, set the amplitude
        // to 1 (full)
        if (curr < mid)
            buffer[n] = 1;
        
        // if the current time is more than
        // the midpoint, set to -1
        else
            buffer[n] = -1;
        
        // add the time for one sample
        curr += sampleTime;
        
        // and check whether we have yet reached the
        // time of one period, after which we must
        // reset our count
        if (curr >= period)
            curr -= period;
    }
    
    return buffer;
}

Waveform II - The sawtooth wave

Sawtooth waves descend from the full amplitude to the lowest amplitude linearly and then jump right back to the maximum. There is a mathematical formula for this, however if we think logically, this involves nothing we need a complex formula for. We simply calculate the period in samples, so the samplerate divided by the frequency, and then divide the range over which we wish to iterate by this value. In our case we want to go from 1 to -1, so we divide 1 - (-1) = 2 by the period in samples to give us the increment, or decrement in this case, per sample. We then go from 1 to -1 and jump right back when we reach -1. 

double * saw(double freq, int bufferLen = 44100)
{
    // the sample buffer
    double * buffer = new double [bufferLen];
    
    // the period, in samples = samplerate / frequency
    double period = 44100 / freq;
    
    // how much we must decrement the count
    // by at each iteration
    // 2.0 because the range is from 1 to -1
    double decr = 2.0 / period;
    
    double curr = 1;
    
    for (int n = 0; n < bufferLen; n++)
    {
        buffer[n] = curr;
        
        // decrement
        curr -= decr;
        
        // reset if we reached the lowest
        // amplitude
        if (curr <= -1)
            curr = 1;
    }
    
    return buffer;
}

Waveform III - The ramp wave

The ramp wave differs in respect to the sawtooth only in the way that it increments from minimum to maximum instead of decrementing from max. to min.

Note: Many people interchange the terms sawtooth waves and ramp waves, so you may see these terms describing the other of the two. I find it more logical to call a ramp wave the one that increments, since it should be/look like a ramp. If ramps would be like a wall and descend behind it, the number of skateboard accidents would increase quite drastically, no?

In terms of the algorithm, the only thing we have to change is our starting point, namely the low-point and then add instead of subtract the increment.

double * ramp(double freq, int bufferLen = 44100)
{
    // the sample buffer
    double * buffer = new double [bufferLen];
    
    // the period, in samples = samplerate / frequency
    double period = 44100 / freq;
    
    // how much we must increment the count
    // by at each iteration
    // 2.0 because the range is from 1 to -1
    double incr = 2.0 / period;
    
    double curr = -1;
    
    for (int n = 0; n < bufferLen; n++)
    {
        buffer[n] = curr;
        
        // decrement
        curr += incr;
        
        // reset if we reached the lowest
        // amplitude
        if (curr >= 1)
            curr = -1;
    }
    
    return buffer;
}

Waveform IV - The Triangle wave

Triangle waves ascend linearly for half the period and then descend for the other half. We have two options for creating this waveform. The first uses the linear incrementing / decrementing methods we just used for the saw and ramp waves and simply switches the increment's sign at the midpoint. The problem is, that we again have the problem that we cannot really represent the the midpoint in samples as there will always be a slight (!) round-off error. A program for this could look like this:

double * tri(double freq, int bufferLen = 44100)
{
    // the sample buffer
    double * buffer = new double [bufferLen];
    
    // the period, in samples = samplerate / frequency
    double period = 44100 / freq;

    // round-off error here
    double incr = 2.0 / (period/2);
    
    double curr = -1;
    
    for (int n = 0; n < bufferLen; n++)
    {
        buffer[n] = curr;
        
        // decrement
        curr += incr;
        
        if (curr >= 1)
        {
            curr = 1;
            incr = -incr;
        }
        
        else if (curr <= -1)
        {
            curr = -1;
            incr = -incr;
        }
    }
    
    return buffer;
}

However, as mentioned, there will be a round-off error. The other method involves a certain mathematical formula which does not introduce any rounding errors and is thus more precise, albeit not as "straightforward". Full courtesy for this algorithm goes to Daniel R. Mitchell, the author of "BasicSynth".

The formula for calculating a triangle wave is:

s = 1 - [ (2 * | ϕ - π |) / π]

Where ϕ is the current phase. 

“We can optimize the code by pre-calculating the value of 2/π and eliminate the subtraction of π by varying the phase from [-π,π].” (Mitchell, 2008, p. 69) 

The code:

double* tri(double freq, int bufferLen = 44100)
{
    double* buffer = new double[bufferLen];
    
    double phase = 0;
    
    // usual phase increment
    double phaseIncr = (twoPI / 44100) * freq;
    
    double triValue;
    
    // precompute this constant part
    double twoDivPi = 2.0/PI;
    
    for (int n = 0; n < bufferLen; n++)
    {
        // non-constant part
        triValue = (phase * twoDivPi);
        
        if (phase < 0) triValue = 1.0 + triValue;
        
        else triValue = 1.0 - triValue;
        
        buffer[n] = triValue;

        phase += phaseIncr;

        if (phase >= PI)
            phase -= twoPI;
        
    }
    
    return buffer;
}

This last function produces a very nicely looking waveform, without the round-off error from the first method.

To conclude, I want to reiterate my point from the main post on this sub-series on complex waveforms, namely that while computing waveforms directly may seem simpler and neater than through Fourier / Additive synthesis, it is nowhere near to what natural tones sound like. Nature is not perfect, nature is not directly computed, which is why you should also check out my post on additive synthesis! However, direct computation is a must for low frequency oscillators or any other modulation, as in these cases even minor imperfections can lead to unwanted modifications of a carrier sound. Therefore, it is really important to know of both methods I describe and use them accordingly.

Feel free to comment questions or suggestions.

The Lanczos sigma factor

When creating complex waveforms through Fourier synthesis, also known as Additive synthesis, it may have occurred to you that the waveforms have a slight ripple and "horn" or overshoot on the left and right of a peak or trough. This is known as the Gibbs phenomenon and is a result of summing only a finite series of sine and cosine waves, opposed to an infinite series as suggested by Joseph Fourier's theorem. Luckily, Hungarian scientist Cornelius Lanczos came up with a solution, namely the so-called Lanczos sigma ( σ ) factor, also known as sigma-approximation. This little factor, added to the amplitude of each partial of a waveform, reduces the Gibbs phenomenon almost entirely (enough for most cases). 


Here is an image of a square wave created by summing 64 partials. You can clearly see the little horns as well as the ripples at the ends of the peaks and throughs.

The sigma factor is defined as:

σ = sin (x) / x

x being:

x = nπ / M

Where n is the current partial number and M is the total partial number.

So, say you want to calculate this factor for the first partial and you want to add 32 partials in total. Take into account here that the fundamental pitch is seen as partial number one, so for this algorithm the total partials are then 33 and the first "real" partial is actually number 2. 

x = 2 * π / 33

x = 0.19039955476

σ = sin (x) / x

σ = 0.99396894386

This means that the amplitude of partial number 2, which we would have had set at something like 0.5 (normalised between  1 and -1), will now be 0.5 * 0.99396894386 = 0.49698447193. If we do this for all the other partials as well, we will almost entirely make the Gibbs phenomenon disappear.

Here is a picture showing the same waveform (64 partials) as in the picture above, but now with the Lanczos factor added:

Convinced? Nevertheless, however great this might seem, it is important to point out that this perfection is not always wanted. I am not saying that the sigma factor takes away from the realism of the sound, as you will still get a very natural, nicely sounding sound for waveforms with high partials, when you are trying to get as close as possible to a perfect square wave while still retaining the natural sound as much as possible. The “problem” with the Lanczos factor is, however, that it straightens out every sound wave, no matter how many partials. Even waveforms with only two partials will be straightened out heavily. 

Here a square wave consisting of a fundamental pitch and two partials:

And here the same waveform with the Lanczos sigma factor:

Wanted? Rarely. Therefore, I suggest that the sigma factor be used only when you really want to get a near perfect waveform. Whenever you’re actually out for the sound of a square wave with so and so many partials, I recommend not using it, as it really then destroys the shape of the additively synthesized waves. 

If you have any questions, suggestions or free cake, feel free to comment. 

Developing a digital synthesizer in C++: Part 3 - Complex waveforms

In the first part of this series, "Developing a digital synthesizer in C++", I showed you how you could create a simple sine wave in code and store it in a sample buffer, after which I also outlined how you can then store this data buffer into a WAV file and actually listen to the data. However, after a while, this sine wave will get a little boring to listen to. The fact is, that there are an infinite number of other sound waves we could have created, a few, as well as our friend, the sine wave, are shown in this picture:

Source:http://cdn.ttgtmedia.com/WhatIs/images/waveform.gif

However, these are not as straightforward to make as calling the sine function on a phase value. Add to that the fact that complex waveforms can be created using various different methods, each yielding more or less the same result. I will explain two frequently used methods that are quite opposite to each other in the approach taken and the outcome received, namely the creation of complex waveforms through addition of sinusoids (sine waves), also called Fourier synthesis, followed by the creation of complex waveforms via direct calculation. As this subject is quite big and equally important, I will split up the implementations of each method into other blog posts, providing the general introduction for each of them here.

At the end of each part of this sub-series on complex waveforms, I will include functions both for summing a variable number of sine waves to give any waveform imaginable through Fourier synthesis, as well as a number of functions for direct calculation of waveforms, without even demanding cake or ice cream in return. 

A quick note on a few naming conventions used in this blog post:

• A partial, so an integer multiple of a fundamental pitch, is the same thing as a harmonic or an overtone. These are just different terms used. I use the word partial. 

• Sawtooth waves and ramp waves are often confused with one another. For me, a sawtooth wave descends from the full amplitude down to the lowest amplitude linearly, before going back to the full amplitude from one sample to the next. Ramp waves, as I refer to them, do the opposite in ascending linearly from the lowest amplitude to the maximum before dropping back down. 

Fourier Synthesis

In the early 1800s French mathematician Joseph Fourier came up with a number of equations that are nowadays known as the Fourier Series. What Monsieur Fourier found was, basically, that any complex waveform can be created simply by summing sine and cosine waves. This is quite astonishing, if you think that even sharp-cornered square waves can be created just by summing together lots of sine waves! The mathematics behind this are naturally quite complex, however the perk of us being programmers is that we do not necessarily need to think in terms of mathematical formulas, but only in terms of applying the theory handed down to us by our fellow mathematicians for us to use pragmatically. 

Here is a picture showing the product of summing a certain number of sine waves to give a rough looking square wave, with a perfect one behind:

To go into more detail, the properties of the above sine waves are that each wave (partial, more precisely) following the fundamental pitch/frequency (first picture) is an odd partial, so 3 times the fundamental frequency, then 5 times the fund. freq., then 7 times etc. and the amplitude of each partial is 1 divided by the partial number (the reciprocal). So partial number one (in this case the fundamental frequency) has an amplitude of 1/1 times the maximum amplitude. The next partial, in this case the next odd partial, which is 3, has an amplitude of 1/3 times the maximum. The next partial of 1/5 etc. 

Brief digression about why the amplitude must decrease with an increasing partial number

The amplitude cannot be the same for each partial because adding sine waves at the same amplitude would quickly lead to overflow of the data type as they would eventually become values that could no longer be represented by the primitive types used to hold them. For example, WAV files use 16 bit signed integers to represent sample values, so ± 32767. If you have an amplitude at that maximum so, e.g. 32767 and add another one at the same amplitude, you would end up with a value that is larger than 32767. As you might know, when signed primitive types overflow, they change their sign, therefore, your sine wave would be normal until it overflows, after which you suddenly have a negative peak!

Theoretical example:

Take a fundamental pitch, so just a plain old sine wave, at a certain frequency, say 20Hz. Add to this the next odd partial so in this case the third partial, which is at 20 * 3 = 60 Hz ( then 5 * 20, 7 * 20, 9 * 20 etc.). As explained, we cannot have all partials at the full amplitude, or else the integer used to hold the sample values is likely to overflow. Therefore, we must have this partial at the reciprocal of the partial number. The amplitude is therefore 1 / 3 the maximum amplitude.

Thus, for any nth partial, given a fundamental frequency f (above I used 20 Hz) and a maximum amplitude a:

• Its frequency is: f * n
• Its amplitude is: a / n 

With each partial added, the waveform created looks more and more like a real, perfect square wave. Here is an image of what the square waves look like depending on how many sinusoids/partials you add (partial number in the legend):

Source: http://www.fourieroptics.org.uk/wp-content/uploads/2012/02/squarewave2.png

I usually use 64 partials for the square waves I want nearest to perfection. It gives extremely crisp and nicely sounding square waves without having to take the partial adding to an extreme. I would not recommend using 200 partials as shown above, as 64 is really more than enough and already more than many other implementations use, however it is finally up to you, of course. Here is what my square wave looks like when summing 64 partials:

Pretty close, right ?! Notice, however, the "horns" at the end of some of the peaks? This is the so called Gibbs phenomenon and is a result of summing a finite series of sine waves and not an infinite series, as Fourier's formula suggests. However, there is a really neat work-around called the Lanczos sigma factor with which we can counter-act it.

Follow this link to get to the more detailed explanation and implementation of creating sine waves through Fourier synthesis.

Direct calculation

Another, apparently easier, method of creating complex waveforms is through direct calculation, so through certain mathematical formulas that compute the exact (I will explain shortly why this is not always good) waveforms wanted. 

The simplest waveform to calculate directly is, of course, the square wave. I'm sure that if you give it a thought, you can already guess how this waveform can be created in code, purely through mathematics. It is nothing else than setting the samples of the first half of the period to the full amplitude and the other half to the lowest amplitude. Here is an example of creating a square wave in the Python programming language (Python looks so close to pseudo code that I'm sure you understand it with any background):

buffer = []

period = 100

midPoint = int(period/2)

for n in range(period):
    
    if n < midPoint:
        buffer.append(1)
        
    else:
        buffer.append(-1)

This code fills the buffer with the maximum amplitude until it reaches the mid point, after which we add one sample at half the amplitude and then continue with the lowest amplitude for the rest of the period.

As you can see, direct calculation is often a lot easier, though this comes at a cost. 

Follow this link to see how to create a few other waveforms directly.

VIP - Very important point (I am so funny)

Finally, I want to point out something very important concerning the creation of complex waves. The fact is, nature is not perfect. Nature is chaos. Therefore, it is essential to understand that the most naturally sounding waveforms, closest to what we would hear in reality, are never the most perfect ones. Thus, Fourier synthesis will always create more naturally sounding waveforms and should usually be consulted if you want to use the waveforms for audio, as the methods involving direct calculation, while from a mathematical standpoint being perfect, will not sound as familiar to the human ear. There are, however, uses for directly calculated waveforms, namely whenever you want to use them for modulation, like with an LFO (low frequency oscillator), as any minor imperfection will make your tone sound different which is therefore usually unwanted. 

Developing a digital synthesizer in C++: Part 2 - Wavefile output

After reading the last part in my series, in which I showed you how to program a basic sine wave and store it in a buffer, you may be wondering how you can get from this array of integers inside your code to ... actually hearing something!

We have two options here:

1. Sending the data buffer to your sound card for it to output it directly. Or:
2. Storing and encoding the data in a certain file format.

While the first option may appeal more as you get faster results, it is also a lot more complex. In this post I will show you how to do the 2nd, namely storing the data in a file. In our case, as the title indicates, we will be dealing with the WAV format.

The WAV format

 The Waveform Audio File Format is a lossless, uncompressed audio file format. It is part of the Resource Interchange File Format (RIFF)  file format family and thus stores it's information in a certain number of chunks. It was originally considered a Windows-only file format, though nowadays this is not true anymore and it is one of the most popular file formats out there.

From a listener's standpoint, this is due to the fact that the audio data stored in the format is uncompressed, meaning every single sample we calculated will be stored in exactly that form, unchanged for virtually all eternity. This results in very good audio quality no matter how often you re-process the samples (given you do not change them). The downside to this increased audio-quality is naturally that also the file size is usually quite a bit higher than that of compressed formats such as the MPEG-3 (mp3) files.

Now from the programmer's viewpoint, the popularity of the WAV file format stems from the fact that it is insanely simple. A .wav file is basically made of three "chunks" of data:

1. RIFF chunk: 
1. ID chunk: a size-4 char array holding the word "RIFF" to identify the file as part of the RIFF family.
2. Chunk size: the size of the entire file in bytes (not just the data, also the "RIFF" char array for example).
3. File format chunk: another size-4 char array that identifies the specific file format of the RIFF Family, in our case this will hold the string "WAVE". 
2. FMT (format) chunk:
1. ID chunk again: size-4 char array "FMT " (notice the trailing space)
2. FMT chunk size (the size of this chunk)
3. - 6, other meta-information about the data buffer like samplerate, number of channels, byte rate and a few other things.
3. DATA chunk:
1. "DATA" char array
2. Chunk size
3. The raw data buffer

I'm not kidding you, this is not some XML format I made up for my latest side-project, but one of the most widely-used and most popular file formats on this planet. If you want to read more on the exact specifications, have a look at this website.

The code

All of these chunks can be implemented in code as the following struct:

    struct waveheader
    {
        uint8_t riff_id[4]; // 'R' 'I' 'F' 'F'
        
        uint32_t riff_size; // chunks size in bytes
        
        uint8_t wavetype[4]; // 'W' 'A' 'V' 'E'
        
        uint8_t fmt_id[4]; // 'f' 'm' 't' ''
        
        uint32_t fmt_size;
        
        uint16_t fmt_code; // 1 = pulse code modulation
        
        uint16_t channels;
        
        uint32_t samplerate;
        
        uint32_t byterate; // bytes per second
        
        uint16_t align; // bytes per sample * channel
        
        uint16_t bits; // one byte per channel, so 16 bits per sample
        
        uint8_t wave_id[4]; // 'd' 'a' 't' 'a'
        
        uint32_t wave_size; // byte total
    };

Our task is it now to insert our data into the members of this struct and write it to a file along with the data buffer we created here. Below is a function that does exactly that (the struct I shared above is declared outside the function). Note that we have to write each sample to the file twice in a row, once for each channel (have a look at the exact specification from the link I shared above to see why).

void writeWav(std::string fname,int16_t * buff, uint32_t dur, uint16_t sr, uint8_t ch)
{
    waveheader wh;
    
    // copy the string "RIFF" into the riff_id char array
    memcpy(wh.riff_id, "RIFF", 4*sizeof(uint8_t));
    
    // same for the wavetype char array
    memcpy(wh.wavetype, "WAVE", 4*sizeof(uint8_t));
    
    // take notice of the trailing whitespace in the string
 memcpy(wh.fmt_id, "fmt ", 4*sizeof(uint8_t));
    
    // fmt_size is 16, the size (in bytes) of the rest of the subchunk
 wh.fmt_size = 16;
    
    /*
    // fmt_code is the audio format used for storing the data. This
    // should usually be 1, which is used for Pulse Code Modulation (PCM).
    // PCM is basically how the file connects the samples we created to 
    // make an actual waveform (curve) out of it.
    */
    
 wh.fmt_code = 1;
    
 wh.channels = (uint16_t) ch;    // 1 = mono, 2 = stereo
    
 wh.samplerate = sr;
    
    // The size of the individual samples. Since we are using
    // 16bit signed integers to represent the sample values,
    // we will use 16 here.
    
    wh.bits = 16;
    
    // Align is the total number of bytes for all channels
 wh.align = (wh.channels * wh.bits) / 8;
    
    // How many bytes per second, there are samplerate
    // (44100) samples per second, each being of size
    // align ((ch * bits )/ 8)
    wh.byterate = (wh.samplerate * wh.align);
    
    // copy the "data" string to the chunk id array
 memcpy(wh.wave_id, "data", 4*sizeof(uint8_t));
    
    uint32_t total_samples = sr * dur;
    
    // number of data bytes in the file, the two
    // stands for the 2 bytes per sample sent per
    // channel
    uint32_t byte_total = total_samples * 2 * ch;
    
    // now that we have ALL the info, we can calculate
    // the total size
    wh.riff_size = byte_total + sizeof(wh) - 8;
    
    // assigning byte_total
    wh.wave_size = byte_total;
    
    // Samples are stored channel for channel, so
    // a single sample has to be written twice in a
    // row, once for each channel.
    
    total_samples *= 2;
    
    int16_t * end_buffer = new int16_t[total_samples];
    
    for (uint32_t n = 0; n < total_samples; n += 2)
    {
        int16_t val = buff[n/2];
        
        end_buffer[n] = val;
        end_buffer[n + 1] = val;
    }
    
    std::ofstream f;
    
    f.open(fname, std::ios::binary | std::ios::trunc);

    if (! f.write(reinterpret_cast(&wh), sizeof(wh)) ||
        ! f.write(reinterpret_cast(end_buffer), byte_total))
        
        throw std::runtime_error("Thy digital oscillation datum did not write properly!");
    
    delete [] end_buffer;
}

The code together with the comments should be pretty self-explanatory and as far as I'm concerned, this code works perfectly. I hope it does so for you too, but please feel free to comment if you are having problems or have any questions about this topic. I am currently working on creating ADSR envelopes as well as directing the sound to the sound card for direct playback and will post either of these steps in the following part of this series.

Developing a digital synthesizer in C++: Part 1 - Generating a sine wave

In this first part of my series "Developing a digital synthesizer in C++", I'll show you how to generate a sine wave and store it into a buffer. With "generating", I only mean calculating the actual values, storing them into an audio file such as .wav is a totally different subject and will be discussed in my following post.

Below you can find the code. I commented on almost everything and also the variable names should more or less tell you what's going on. 

#include <cmath>
#include <stdint.h>

#define twoPI 6.28318530717958 // 2*pi = 360˚ = one full cycle

// Using standard typedefs for portability, you can change them to normal data types if you like

// @param dur: I was indecisive about whether the duration
// should be 32 or 16 bits aka 130 years or 18 hours, but then
// again you never know what sound installations people come up with

int16_t* gen(float freq, uint32_t dur, float vol=1.0);

int main(int argc, const char * argv[])
{
    int16_t* buff = gen(440,3);
    
    /* Do something with it */
    
    delete [] buff;
}

int16_t* gen(float freq, uint32_t dur, float vol)
{
    uint32_t samplerate = 44100;    // samples per second
    
    // initial phase, you could offset it, but then again you could not
    double phase = 0;
    
    // The phase increment is the phase value the phasor increases by
    // per sample.
    double phaseincr = twoPI/samplerate * freq;
    
    // The amount of samples the buffer must hold
    uint32_t total_samples = samplerate * dur;
    
    int16_t *buffer = new int16_t [total_samples]; // grab a new array with size for the entire buffer
    
    for (int i = 0; i < total_samples; i++) // fill the buffer
    {
        // the factor 32767 comes from the fact that .wav files store samples
        // as 16 bit signed integers, before this the values are normalized
        // between -1 and + 1 (the sine values). If you want to do something
        // else with these values before storing them somewhere I recommend leaving        // the factor away.
        buffer[i] = 32767*(sin(phase) * vol);
        phase += phaseincr;
        
        // when the phasor hits 2Pi/360˚/full circle we have to reset the phase
        if (phase >= twoPI) phase -= twoPI;
    }
    
    return buffer;
}

You could change the values for the sample rate, initial phase or volume, but especially changing the samplerate is not recommended as 44.1 kHz is really quite standard. I actually wanted to write a bit about the code, but there's nothing that's not already explained in the comments. If you do have any issues/questions or the code doesn't build for you, feel free to comment.

Developing a digital synthesizer in C++: Part 0 - Introduction

I will be graduating from my high school in about one and a half years and for my diploma, I am required to write some form of research paper called a "Vorwissenschaftliche Arbeit" or "Pre-scientific paper", to give you a bad English translation. This is due to a major education reform that's taking place in Austria and my year is the first to be affected by it. Anyway, my topic is, as the title states, "Developing a digital synthesizer in C++". I'm planning on finishing the actual software as of 2014, and then I have another half year to write about it. In the meantime, I will write lots of articles on the different experiences I encounter and post code-snippets under "an" open source license , should anyone seek to do something similar.