# 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;
}
```