BLEP

Overview

The BLEP algorithm implements a set of oscillators using Polynomial Bandlimited Step functions, also known as PolyBLEP.

Much of this is based on the implementation found on blog post by Martin Finke, with a few adjustments.

Due to their bandlimited properties, these oscillators are great things to reach for when working for a sound source to use with subtractive synthesis techniques. Bandlimited things sound better because they reduce aliasing, an audible artifact in the sound that happens when a signal plays frequencies that are a little too high (there's a lot of resources on aliasing, so this is pretty much all I'm going to say on this).

Algorithm Overview

BLEPs aim to create better versions of what we would call wavetable oscillators, or table-lookup oscillators like osc or oscf. You can think of these methods as taking a single waveform and repeating it a bunch of times to produce a sound. Changing the length of that waveform changes the frequency. Changing the shape of the waveform changes the timbre. Famous oscillator names like sawtooth, square, pulse, and triangle all their names from the appearance of their waveform.

For reasons beyond the scope of this document, these table-lookup wavetable oscillators often produce a great deal of unwanted noise known as aliasing. A lot of sources of aliasing occur when large discontinuities, or jumps, happen in the waveform. Square and pulse waves, for example, make a giant jump from a high state to a low state. BLEPs work by finding these large discontinunuities in classic waveform shapes, and then smoothing them out a little bit using a polynomial curve. It's not a perfect process, but it does a pretty decent, especially from a perceptual standpoint.

Tangled Files

blep.c and blep.h.

<<blep.c>>=

#include <math.h>

#define SK_BLEP_PRIV
#include "blep.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

<<static_funcdefs>>
<<funcs>>

<<blep.h>>=

#ifndef SK_BLEP_H
#define SK_BLEP_H

#ifndef SKFLT
#define SKFLT float
#endif

<<typedefs>>
<<funcdefs>>

#ifdef SK_BLEP_PRIV
<<structs>>
#endif

#endif

The PolyBlep

This algorithm centers around a tiny function called polyblep. It takes in a time value t representing the position in the waveform in a normalized range, and the value dt, the delta time between samples.

This will apply two different polynomial curves if the position is at the beginning or ends of the position.

<<static_funcdefs>>=

static SKFLT polyblep(SKFLT dt, SKFLT t);

<<funcs>>=

static SKFLT polyblep(SKFLT dt, SKFLT t)
{
    if (t < dt) {
        t /= dt;
        return t + t - t * t - 1.0;
    } else if(t > 1.0 - dt) {
        t = (t - 1.0) / dt;
        return t * t + t + t + 1.0;
    }

    return 0.0;
}

Initialization and Struct

sk_blep is the struct.

<<typedefs>>=

typedef struct sk_blep sk_blep;

<<structs>>=

struct sk_blep {
    <<sk_blep>>
};

Initialize with sk_blep_init.

<<funcdefs>>=

void sk_blep_init(sk_blep *blep, int sr);

<<funcs>>=

void sk_blep_init(sk_blep *blep, int sr)
{
    <<init>>
}

Components

Frequency Value

The frequency uses parameter caching.

<<sk_blep>>=

SKFLT freq;
SKFLT pfreq;

<<init>>=

sk_blep_freq(blep, 1000);
blep->pfreq = -1;

Onedsr

The onedsr constant is 1/sr.

<<sk_blep>>=

SKFLT onedsr;

<<init>>=

blep->onedsr = 1.0 / sr;

Phasor Values

Like any good oscillator, under the hood there is a phasor. The phs keeps track of the phase, and the inc incrementor keeps track of the increment.

<<sk_blep>>=

SKFLT inc;
SKFLT phs;

<<init>>=

blep->inc = 0;
blep->phs = 0;

This is another small change from Finke's original implementation. Using a normalized phasor range instead of one that goes between 0 and 2 pi simplifies the computation.

Leaky Integrator

For the triangle wave, a leaky integrator will be used. We will use a very small pole value of 100ms as the filter coeffiecient A. This value was empirically chosen as a reasonably close value to 1.

<<sk_blep>>=

SKFLT A;
SKFLT prev;

<<init>>=

blep->A = exp(-1.0/(0.1 * sr));
blep->prev = 0;

Note: Finke's original implementation uses the increment value as the filter's coefficient, and it's unclear to me why. So I've gone with something I can better understand and reason with.

DC Blocker

That pesky triangle! The leaky integrator it uses introduces some serious DC. A DC blocking filter is used to remove this.

<<sk_blep>>=

SKFLT R, x, y;

The DC blocking coefficient R has been chosen to be close to 0.99 (a common DC blocker coefficient value) when the sampling rate is 44.1kHz.

<<init>>=

blep->R = exp(-1.0/(0.0025 * sr));
blep->x = 0;
blep->y = 0;

Setting The Frequency

The frequency of the oscillator is set with sk_blep_freq.

<<funcdefs>>=

void sk_blep_freq(sk_blep *blep, SKFLT freq);

<<funcs>>=

void sk_blep_freq(sk_blep *blep, SKFLT freq)
{
    blep->freq = freq;
}

Core Tick Function

The core computation is done with a static function called tick. It's a generalized function that takes in a callback for each waveform.

<<static_funcdefs>>=

static SKFLT tick(sk_blep *blep,
                  SKFLT (*wave)(sk_blep *, SKFLT));

<<funcs>>=

static SKFLT tick(sk_blep *blep,
                  SKFLT (*wave)(sk_blep *, SKFLT))
{
    SKFLT out;

    out = 0.0;

    <<update_increment>>
    <<compute_wave>>
    <<update_phasor>>

    return out;
}

To begin, the increment value is updated if the frequency is changed.

<<update_increment>>=

if (blep->freq != blep->pfreq) {
    blep->pfreq = blep->freq;
    blep->inc = blep->freq * blep->onedsr;
}

The wave callback gets used to compute the actual BLEP'd sample.

<<compute_wave>>=

out = wave(blep, blep->phs);

To wrap up, the internal phasor is updated.

<<update_phasor>>=

blep->phs += blep->inc;

if (blep->phs > 1.0) {
    blep->phs -= 1.0;
}

Sawtooth

A sawtooth BLEP is produced with sk_blep_saw.

<<funcdefs>>=

SKFLT sk_blep_saw(sk_blep *blep);

<<funcs>>=

SKFLT sk_blep_saw(sk_blep *blep)
{
    return tick(blep, blep_saw);
}

The sawtooth is the most straightforward BLEP. The phasor value already produces the correct shape. It just needs to be scaled to be in range -1 to 1. This signal is then fed into the blep function to smooth out the edges.

<<static_funcdefs>>=

static SKFLT blep_saw(sk_blep *blep, SKFLT t);

<<funcs>>=

static SKFLT blep_saw(sk_blep *blep, SKFLT t)
{
    SKFLT value;

    value = (2.0 * t) - 1.0;
    value -= polyblep(blep->inc, t);

    return value;
}

Square

A square wave BLEP is computed sk_blep_square.

<<funcdefs>>=

SKFLT sk_blep_square(sk_blep *blep);

<<funcs>>=

SKFLT sk_blep_square(sk_blep *blep)
{
    return tick(blep, blep_square);
}

The square shape is first derived by splitting the phasor signal in half: lower half is -1, the upper half is 1.

The blep is rounded on both edges of the square, so the BLEP gets called twice. The fmod(t + 0.5) is a trick to offset the value in the right position.

<<static_funcdefs>>=

static SKFLT blep_square(sk_blep *blep, SKFLT t);

<<funcs>>=

static SKFLT blep_square(sk_blep *blep, SKFLT t)
{
    SKFLT value;

    value = t < 0.5 ? 1.0 : -1.0;
    value += polyblep(blep->inc, t);
    value -= polyblep(blep->inc, fmod(t + 0.5, 1.0));

    return value;
}

Triangle

A triangle BLEP is generated with sk_blep_triangle.

<<funcdefs>>=

SKFLT sk_blep_triangle(sk_blep *blep);

<<funcs>>=

SKFLT sk_blep_triangle(sk_blep *blep)
{
    return tick(blep, blep_triangle);
}

A triangle wave BLEP is computed by calculating the integral of a square wave. When the square wave is at the lower state, it goes down. When it is at the higher state, it goes up.

Integration is a fancy way of saying "sum it all up". Left to their own devices, a integrated square wave would produce triangle waves with huge amplitudes proportional to their wavelength (in samples). Dividing by this wavelength will normalize the samples.

Integration in floating point numbers can eventually result in weird numerical errors. As a preventative measure, The summation is designed to "forget" about previous values over time, creating what is known as a leaky integrator.

BUT, this numerical forgetfulness comes at a cost of some initial DC offset at the beginning. This can be mostly mitigated with a DC blocking filter.

<<static_funcdefs>>=

static SKFLT blep_triangle(sk_blep *blep, SKFLT t);

<<funcs>>=

static SKFLT blep_triangle(sk_blep *blep, SKFLT t)
{
    SKFLT value;

    /* compute square */
    value = t < 0.5 ? 1.0 : -1.0;
    value += polyblep(blep->inc, t);
    value -= polyblep(blep->inc, fmod(t + 0.5, 1.0));

    /* scale and integrate */
    value *= (4.0 / blep->freq);
    value += blep->prev;
    blep->prev = value;

    /* dc blocker */
    blep->y = value - blep->x + blep->R*blep->y;
    blep->x = value;

    return blep->y * 0.8;
}