PolyBLEP: two samples that fix oscillator aliasing

Abstract

A sawtooth generated as a raw ramp contains a discontinuity, and a discontinuity in discrete time contains harmonics above the Nyquist frequency that fold back into the audible band as inharmonic partials. Measured on a 65536-point FFT at 48 kHz, a naive ramp at 1199.71 Hz carries alias energy 15.15 dB below its wanted harmonics; adding PolyBLEP, a two-sample polynomial correction applied at each wrap, moves that to 31.71 dB for roughly six arithmetic operations per sample — better than ideal four-times oversampling, which reaches 27.91 dB at four times the oscillator cost.

Write an oscillator the obvious way and it sounds wrong. Not dull, not thin — wrong, in a way that gets worse as you play higher, and that no filter afterwards can remove. The cause is four characters long: the wrap at the end of the ramp.

Every number in this article is recomputed on each build of this site by rendering the oscillators and taking a 65536-point FFT at a 48 kHz sample rate. The code is the synthesiser in the demoscene engine.

The naive oscillator, and why it is wrongPermalink to “The naive oscillator, and why it is wrong”

/* the whole of a naive sawtooth */
v = 2.0f * ph - 1.0f;
ph += dt;
if (ph >= 1.0f) ph -= 1.0f;

The ramp is fine. The wrap is not: once per cycle the output jumps from +1+1 to 1-1 instantaneously, and an instantaneous jump has a spectrum that falls off only as 1/k1/k and never reaches zero. In continuous time that is what a sawtooth is, and it is harmless. In discrete time every harmonic above the Nyquist frequency has to go somewhere, and where it goes is back into the band, mirrored around fs/2f_s/2:

falias=kf0nfsf_{\text{alias}} = \left| k f_0 - n f_s \right|

Nothing in the code is wrong. The signal being sampled simply has more bandwidth than the sample rate can carry, and sampling it anyway folds the excess back on top of the part you wanted.

Where the images landPermalink to “Where the images land”

Take a note at 1199.71 Hz against a 48 kHz sample rate. Twenty harmonics fit below Nyquist: the twentieth sits at 23994.14 Hz. The twenty-first would sit at 25193.85 Hz, which does not exist at this sample rate, so it appears instead at

4800025193.85=22806.15 Hz48000 - 25193.85 = 22806.15\ \text{Hz}

That number is the whole problem in one line. It is not a multiple of the fundamental — the nineteenth harmonic is at 22794.43 Hz — so the image lands 11.72 Hz away from a real partial, and two tones 11.72 Hz apart beat at 11.72 Hz. The ear hears a slow warble on a note that should be steady, and a shimmer that changes character with every note played, because each note folds its images to a different place.

This is why aliasing is not the same complaint as “too bright”. Brightness is harmonic and sits where the ear expects it. Aliasing is inharmonic, moves in the opposite direction to the note played, and is the single most reliable way to make a synthesiser sound like software.

Measuring itPermalink to “Measuring it”

Claims about aliasing are easy to make and easy to fake, so the measurement here is pinned down:

  • The oscillator runs for 65536 samples at 48 kHz, and the frequency is chosen to fall exactly on an FFT bin, so no window is needed and there is no leakage to hide behind.
  • Every partial, wanted or folded, lands on an exact bin, so the spectrum separates cleanly into two sets: harmonics genuinely below Nyquist, and images of harmonics above it.
  • The reported figure is total alias energy relative to total harmonic energy, in decibels.

As a check that the analysis is measuring what it claims, the fundamental of the naive ramp comes out at 0.318310, which is 1/π1/\pi to six decimals — the value the Fourier series of a unit sawtooth requires, and a number no bug in the binning would produce by accident.

NoteHarmonics in bandNaive rampWith PolyBLEPGain
199.95 Hz12022.95 dB39.12 dB16.17 dB
599.85 Hz4018.17 dB34.49 dB16.32 dB
1199.71 Hz2015.15 dB31.71 dB16.56 dB
2400.15 Hz911.66 dB26.27 dB14.61 dB
3999.76 Hz69.87 dB27.69 dB17.81 dB
6999.76 Hz36.81 dB22.87 dB16.06 dB

Read the columns as “how far the rubbish sits below the music”. At 199.95 Hz the naive ramp is tolerable; at 6999.76 Hz the alias energy is only 6.81 dB below the harmonics, which is to say the oscillator is producing almost as much noise as signal. That is the mechanism behind the folk knowledge that software synthesisers “fall apart in the top octave”.

The fix, and where the polynomial comes fromPermalink to “The fix, and where the polynomial comes from”

The band-limited step, or BLEP, is the properly filtered version of a jump: the step you would get if the discontinuity had been low-passed before sampling. Subtract the ideal step from it and what remains is a short residual concentrated around the discontinuity. Adding that residual to the naive signal turns the raw jump into a band-limited one.

The exact residual is infinitely long. PolyBLEP truncates it to two samples and fits a quadratic to each half, which is the cheapest approximation that gets the first-order behaviour right:

static inline f32 polyblep(f32 t, f32 dt)
{
    if (t < dt)        { t /= dt;         return t + t - t * t - 1.0f; }
    if (t > 1.0f - dt) { t = (t - 1.0f) / dt; return t * t + t + t + 1.0f; }
    return 0.0f;
}

t is the phase in [0,1)[0,1) and dt is the phase increment, so dt is one sample expressed in cycles — at 1199.71 Hz and 48 kHz it is 0.024994 of a cycle. The two branches are the same parabola seen from either side of the wrap. Writing uu for the position in samples relative to the discontinuity:

before (u[1,0)):(u+1)2,after (u[0,1)):(1u)2\text{before } (u \in [-1,0)): \quad (u+1)^2, \qquad \text{after } (u \in [0,1)): \quad -(1-u)^2

The residual therefore rises smoothly from 0 to 1 over the sample before the jump and continues from 1-1 back to 0 over the sample after it, which is exactly what replaces a vertical edge with one that takes two samples to fall. That is the whole idea: the discontinuity is not removed, it is given a slope the sample rate can represent.

Applied to the sawtooth, it is one subtraction:

case DM_OSC_SAW:
    return 2.0f * ph - 1.0f - polyblep(ph, dt);

What it costs, and what it buysPermalink to “What it costs, and what it buys”

A branch, a divide and four arithmetic operations per sample, on the two samples per cycle where the correction is non-zero and a single failed comparison everywhere else. For that, alias energy at 1199.71 Hz drops from 15.15 dB below the harmonics to 31.71 dB below them — 16.56 dB of improvement — and the loudest single image, the one at 22806.15 Hz, falls from 36.39 dB below full scale to 45.13 dB below it.

The obvious alternative is oversampling: run the oscillator faster, filter, decimate. Held to the same measurement, with an ideal brick-wall filter that no real implementation achieves:

MethodAlias energy below harmonicsCost per output sample
Naive ramp15.15 dB1 oscillator
Ideal 2× oversampling21.77 dB2 oscillators, plus a filter
Ideal 4× oversampling27.91 dB4 oscillators, plus a filter
PolyBLEP31.71 dB1 oscillator, plus about six flops

PolyBLEP beats ideal four-times oversampling at this note while costing a twentieth as much, because it attacks the shape of the spectrum rather than pushing the fold point further out. That is not an argument against oversampling — the two compose, and a production synthesiser under heavy modulation usually wants both — but it settles which one to write first.

The correction is also nearly free in the sense that matters musically: the fundamental measures 0.318310 without it and 0.317656 with it, a change of 0.018 dB, so the note does not get quieter or duller. What it removes is only what should never have been there.

The square wave needs two of themPermalink to “The square wave needs two of them”

A square has two discontinuities per cycle, one at the start and one at the pulse width, so it takes two corrections with opposite signs:

case DM_OSC_SQUARE: {
    f32 v  = ph < pw ? 1.0f : -1.0f;
    f32 p2 = ph + 1.0f - pw;
    if (p2 >= 1.0f) p2 -= 1.0f;
    return v + polyblep(ph, dt) - polyblep(p2, dt);
}

p2 is the phase rotated so that the second edge sits at the wrap, which is the only place polyblep looks. Getting this wrong is the classic pulse-width-modulation bug: the correction tracks the leading edge, the trailing edge stays raw, and the aliasing that survives sweeps with the modulation.

Why the triangle gets nothingPermalink to “Why the triangle gets nothing”

A triangle is continuous. Only its slope jumps, so its harmonics fall as 1/k21/k^2 rather than 1/k1/k, and the images that fold back are already 20 dB or more further down. The engine generates it with no correction at all:

case DM_OSC_TRI:
    return 4.0f * fabsf(ph - 0.5f) - 1.0f;

The corresponding fix, if the top octave ever demands one, is a BLAMP — a band-limited ramp, the same construction integrated once more, correcting a slope discontinuity instead of a value discontinuity. It has not been needed here, which is worth saying plainly rather than implying the omission is an oversight.

What PolyBLEP does not fixPermalink to “What PolyBLEP does not fix”

  • Very high notes. At 6999.76 Hz the corrected oscillator is still only 22.87 dB clean, because a two-sample window is a smaller fraction of a shorter cycle. Above a few kilohertz a wavetable with per-octave band-limited tables is the better structure.
  • Hard sync and phase distortion. These create discontinuities at phases the oscillator does not know about in advance. Corrections there need the exact sub-sample position of the event, which is what minBLEP is for.
  • Frequency modulation at audio rates. The residual assumes the phase increment is constant across the two samples it spans.
  • Anything after the oscillator. A distortion stage or a hard clipper generates its own harmonics past Nyquist, and no amount of care upstream helps. That is what the oversampling argument is actually for.

The general lessonPermalink to “The general lesson”

The reason this is worth 1500 words in a demo engine is that it is the same lesson as choosing a step size for a raymarcher and the same one as fixing the summation order in a solver: the defect is not in the code, it is in the gap between the continuous object being described and the discrete machine describing it. Nothing in the naive oscillator is a bug. It is a correct implementation of a signal that cannot be sampled, and the fix is to change the signal rather than to hunt through the implementation.

Measure it before and after, in decibels, on a spectrum you can point at. Six flops is a cheap answer, but only once you know what question it was answering.

Related