Complex Waveforms II - Direct Waveform Calculation

We can easily calculate some waveforms directly. For example, we can calculate the samples for a sawtooth wave directly from the phase.

The phase varies from [0,2π], thus the first term varies from [0,2] and subtracting 1 produces values over the range [-1,1].

phaseIncr = (twoPI / sampleRate) * frequency;
for (n = 0; n < totalSamples; n++) {
    sample[n] = (phase / PI) - 1;
    phase += phaseIncr;
    if (phase >= twoPI)
        phase -= twoPI; 
}

A faster method is to produce the sawtooth by calculating an amplitude increment directly from the number of samples per period. The amplitude increment per sample is the range of [‑1,+1] divided by the samples per period, i.e:

 

This equation produces a very efficient program that only needs to increment the amplitude for each sample and test for the end of the period.

sawIncr = (2 * frequency) / sampleRate;
sawValue = -1
for (n = 0; n < totalSamples; n++) {
    sample[n] = sawValue;
    sawValue += sawIncr;
    if (sawValue >= 1)
        sawValue -= 2; // or sawValue = -1;
}

A triangle wave is a linear increment or decrement that switches direction every π radians. Thus, the amplitude covers the range twice within a period.

This is easily implemented directly, but we can optimize the calculation by pre-calculating the value of 2/π, and then using a conditional test for the absolute value function. Likewise, we can eliminate the subtraction of π by varying the phase from [-π, π].

phaseIncr = (twoPI / sampleRate) * frequency;
twoDivPI = 2.0/PI;
for (n = 0; n < totalSamples; n++) {
    triValue = phase * twoDivPI;
    if (triValue < 0)
        triValue = 1.0 + triValue;
    else
        triValue = 1.0 - triValue;
    sample[n] = triValue;
    phase += phaseIncr;
    if (phase >= PI)
        phase -= twoPI;
}

For a square wave, we only need to check if the phase is less than or greater than the midpoint of the period and use +1 or -1 as the value. We can produce a variable duty cycle pulse wave by multiplying the midpoint by the duty cycle.

phaseIncr = (twoPI / sampleRate) * frequency;
midPoint = twoPI * (dutyCycle / 100.0);
for (n = 0; n < totalSamples; n++) {
    phase += phaseIncr;
    if (phase >= twoPI)
        phase -= twoPI; // or phase = 0;
    if (phase >= midpoint)
        sample[n] = -1;
    else
        sample[n] = 1;
}

When directly calculating the square or sawtooth waves, the waveform will transition from peak to peak at the point where the phase wraps. Because of the sudden transition at the end of the period, any accumulated round-off in the phase will shift the beginning of the waveform back and forth by one sample ("phase jitter"). In effect, we are sampling the angular velocity at a rate that is too low, and thus will produce a sub-harmonic alias frequency.  At low frequencies, the alias frequency is not heard, but at higher frequencies it is readily apparent. We can eliminate the phase jitter by wrapping the phase to zero each time it overflows, but this produces a slight frequency error and the oscillator goes "out of tune" as the frequency increases. The amount of error for any given frequency is the fractional part of the sample count divided by the number of samples in a period. The period is longer at low frequencies and the frequency error is small, but increases as frequency increases. The accurate audio range for these oscillators is only a portion of the audible frequency range. Thus, direct calculation of a square or sawtooth wave should be limited to use as an LFO for vibrato, tremolo, panning effects, etc.

One other caution should be noted. A straight up or down change in a waveform indicates a large number of partials, including frequencies beyond the Nyquist limit. This can create problems when the resulting signal is passed through various signal processing functions. This is another reason why waveforms directly calculated in this manner are best limited to low frequencies or as control signals.

Links:

Dan Mitchell's Personal Website