Making DPW-ish integrated wavetable synthesis work well with PM

By Joel K. Pettersson. Added 2022-08-31. Updated 2022-09-01.

An older anti-aliasing technique for audio oscillators uses integration of the underlying wave type followed by differentiation for each sample to be produced. Most well-known is "differentiated parabolic wave" (DPW) sawtooth wave synthesis, where the integral of a sawtooth – a parabola – is differentiated. With a 6 dB per octave aliasing reduction result, this generally does less at a greater computational cost than the more modern PolyBLEP anti-aliasing approach – but it can be used in some additional ways, including for arbitrary stored waveforms handled by the same code, applied to angle modulation in an oscillator (FM and PM), and also for waveshaping.

It can also be layered several times (i.e. using higher-order integration and dfferentiation), each time increasing the aliasing reduction by 6 dB per octave, but I'll focus on the simple version here. The simple version becomes more complicated when tweaking a wavetable oscillator using it to correctly handle arbitrary through-zero angle modulation – and trying to combine that with several layers is even trickier, and doesn't seem worth it. Here I'll describe my experience getting a then-new oscillator in saugns v0.3.9 to work well (I missed a few bugs and the v0.3.12 version will be more definite).

For more general information, here's some articles about the technique in general. The first-listed author, T.Rochebois, relates having used it for wavetable synthesis in the 1990s and it being a known "trick" back then; he also refers to it here and there on the KVR Audio forum, without describing how to iron out the remaining issues for more fully general modulation use. The other academic sources don't deal with use for modulation at all, and make more general audio use sound a bit simpler than it is. Outside audio, the same basic thing is also done in image processing and is then referred to as using a "summed area table".


Simple waveform reproduction

For each output sample, with this technique two look-up table (LUT) values are compared and the difference used. The difference is scaled down according to the distance, or difference in phase, between the values. Each LUT value represents the sum of the original non-integrated values up to that point, and differences and a division are used to effectively average out what changes between two parts of a waveform instead of picking two isolated values and missing the impact of what's in-between. But this averaging, sometimes called box filtering, only reduces aliasing to a certain extent, it does not eliminate it.

It is easiest to get an oscillator of this type to produce a waveform when the frequency is not very low. The DAFx papers don't mention one reason to always use at least linear interpolation for reading values from the look-up table (LUT), what happens as frequency drops far; with no interpolation, nearby phase values will round to the exact same LUT position, and the difference between successive values – and so the output – will be zero most of the time, with rarer bursts of great amplitude when the full difference between values suddenly hits. This matters a lot when arbitrary frequencies are to be allowed, but it often doesn't matter for restricted ranges.

Other than that, reproducing a waveform with no modulation in this way is pretty forgiving. It's natural that for the highest frequencies, the amplitude becomes a little lower, the extent of this varying per waveform. The half-sample delay added to the output can be compensated for (technically, overcompensated) by incrementing phase prior to using it each sample instead of incrementing it after using it; the result "looks" much more like no delay. On numerical precision, unlike when doing something more fancy it seems enough to use single-precision floating point, for both LUTs and the calcuations afterward, though there may be edge cases rarely hit which can only be fixed by making the oscillator modulation-proof with some double-precision usage.

For multi-layered versions giving greater aliasing reduction, eventually both fuller use of double precision and higher-quality interpolation become a must, and even with that done it remains trickier to also get angle modulation working stably. Trying to do more can also become computationally expensive enough to seem pointless for most practical purposes, and I gave up the project of trying to do everything mentioned in my saugns program. I stuck to one layer and modulation-friendly behavior.

Restoring removed DC offset

The technique requires the absence of DC offset in the stored waveforms, prior to integration and by extension also after. An oscillator like this requires some help when an output waveform with DC offset is actually wanted – such as storing the removed DC offset as one look-up value for each possible LUT selection, and adding it per sample after differentiation. This may be worthwhile for oscillators that are meant to produce an output within a specific amplitude bounds, e.g. +/- 1.0, so as to fully use the range.

A related question is, what happens for a frequency of exactly 0 Hz? An oscillator like this can only differentiate when the phase distance is non-zero, special handling being needed for the zero case. You could e.g. have a 0 amplitude output for that special case, but it would not be consistent with how the oscillator functions at nearby frequencies (and would also give strange results if PM is thrown into the mix). For through-zero behavior, the result should depend on the current phase; for example, for a square wave, there's two main 0 Hz output values depending on the phase, 0 Hz meaning to keep repeating the same output.

To repeat the result when there's no phase difference, the simple option is to just store the result of differentiating when differentiating, and copy the stored value when not. To handle edge cases including first results from the oscillator, an initial fallback value can be made for it before running after configuring it. (A fallback differentiation result is easily calculated by comparing two adjacent LUT values, once the LUT and starting phase are decided.)

However, details such as these could be tweaked differently, depending on what's wanted from an oscillator and the feature combination expected to work consistently. In my case, I wanted to keep most of the behavior of an older naive oscillator, and reproduce the same set of waveforms with centered peak amplitudes.

Reducing aliasing for angle modulation

Two values at some distance (the phase difference) are used each output sample, to retrieve the difference between the values, then divided by the phase difference; for a constant frequency, that phase difference will remain constant, while for a varying frequency, it will vary. With phase modulation, there's two different ways to do it here; to also reduce aliasing for PM, the PM input should be included in the difference handled in the differentiation; alternatively, the PM input can be placed outside that difference (by moving both of the values compared), to keep PM in particular naive.

If the PM phase offset is applied for just one of the two values diffed each sample, and the corresponding scaling factor is also updated to match, then there will be reduced aliasing for the PM. This actually makes for the shortest code, but any tiny mistake can make samples sometimes pop with mismatched amplitude. Such mistakes include, I've concluded from my testing, continuing to use single-precision floating point for the calculation; it's still fine to use for the LUTs themselves, but beginning with the interpolation when reading values from a LUT, double precision is needed until after the differentiation is done. This problem can also become noticeable with FM, which by varying the frequency likewise varies how diffed values are related in phase.

It may be easier (but I don't recommend it) to avoid such click/pop glitches for PM in particular by not anti-aliasing it, each sample applying the PM input to the phase for both values to compare and diff at the same time (the scaling factor then being unchanged as the phase difference between the two values is unchanged). But this requires another interpolated table look-up each sample (since both values move instead of one being shared across sample pairs) apart from it accomplishing less in terms of not reducing PM aliasing. The possibility of doing this is mostly a curiosity, making both fully aliased and reduced-aliasing angle modulation possible, even potentially at the same time.

Dealing with large modulator amplitudes

After getting the first solution working a year ago in saugns v0.3.9, I still missed one remaining severe problem with the oscillator. When PM modulator amplitudes were larger, significant distortion resulted in the oscillator handling that input, sounding very similar to bitcrushing distortion, and looking on visual examination as if a bitcrushed shape had been mixed with a smooth shape. The greater the PM modulation amplitude, the stronger the effect, giving a very "retro" sound of the wrong kind for my own purposes.

This problem completely disappears when using higher-quality interpolation for reading from the LUTs; specifically, I changed it from linear to Hermite, 4 values. Various other better-than-linear interpolation choices may also work well. Unfortunately, with that much-needed change, performance dropped significantly, and the oscillator become slow enough that I've been considering experimenting with something else in its place later.

In conclusion, for use with modulation

This type of oscillator is interesting in part because the aliasing reduction also applies to angle modulation. But the whole reduction is comparatively weak by modern standards, and with all adjustments made which are needed for problem-free output with wide-ranging parameters, it becomes attractive to experiment with other solutions, mainly to remove the need for at least one instance of Hermite interpolation or similar drawing on 4 LUT value reads per sample. (Audio rendering in my very oscillator-centric saugns v0.3.x program is generally a little less than half as fast with this oscillator, compared to naive oscillator code which also uses simpler linear interpolation. Other aliasing reduction schemes could do more for that cost, except likely not work for modulation.)

With only 6 dB/octave aliasing reduction, decent aliasing reduction depends on using a high-enough sample rate. For general audio purposes, a sample rate of at least 96 kHz seems sensible if there is to be much of a point – and then the difference for very intense angle modulation is also easier to hear, beyond the most obvious "as if high-pass filtered" effect on noise from aliasing. (And for some types of sound design, where the aliasing from FM & PM is taken into account and expected, this may be the wrong type of oscillator to use.) A higher sample rate also pushes the fairly mild high frequency roll-off that goes with the aliasing reduction upwards in the frequency spectrum (which would be even more important for some possible alternatives like PolyBLEP).