**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.

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) / sigmaThis 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

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

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:

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:

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!

## No comments :

## Post a Comment