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.

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.