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.