4 · osc (part ii) — Fourier series & wavetables

The naive saw and square from part i sound great low but turn harsh and metallic above ~2 kHz. That harshness is aliasing. Fixing it teaches us the single most important idea in audio DSP — the frequency domain / Fourier series — and a classic tool, the wavetable. The math here is presented intuitively; every formula has a plain-language note.

1 · Why the naive saw aliases

In wave we learned an ADC can only capture frequencies below Nyquist ( $f_s/2$ ); anything higher folds back to a wrong, lower frequency. Our generation recipe was “imagine a waveform, sample it.” The flaw: if the waveform itself contains frequencies above Nyquist, sampling produces garbage. So the fixed recipe is:

imagine a waveform,
remove its content above Nyquist,
sample it.

Step 2 is this chapter. But to “remove high-frequency content” we first need to know what frequencies a sawtooth even contains — even though it has a single repeat rate. Enter Fourier.

2 · The frequency domain: Fourier series

The one idea: any periodic signal is a sum of sine waves whose frequencies are integer multiples of a fundamental. Those multiples are the harmonics. A sine has one harmonic; a saw has infinitely many — which is exactly why it aliases.

This is like a Taylor series, but built from sinusoids instead of polynomials — because for sound, frequency (pitch) is what matters, and harmonics are how real strings and pipes naturally vibrate. For a 1-periodic signal:

$x(\phi) = a_0 + \sum_{k=1}^{\infty}\Big(a_k\cos(2\pi k\phi) + b_k\sin(2\pi k\phi)\Big)$

The numbers $a_k, b_k$ are the Fourier coefficients — the “recipe” saying how much of harmonic $k$ to add.

The sawtooth’s recipe

For the sawtooth, the cosine terms vanish (it is an odd function) and the sine coefficients work out to a strikingly simple pattern:

$\text{saw}(\phi) = \sum_{k=1}^{\infty} (-1)^{k+1}\,\frac{2}{\pi k}\,\sin(2\pi k\phi)$

Math note — reading the saw series. Harmonic $k$ has amplitude $\frac{2}{\pi k}$ : the 2nd is half the 1st, the 3rd a third, and so on — they fade slowly. The $(-1)^{k+1}$ just flips the sign every other harmonic. “Slowly fading” is the problem: a saw at 16 kHz with $f_s = 48$ kHz already has its 2nd harmonic at 32 kHz, well above the 24 kHz Nyquist, so it folds back down and pollutes the tone. Keep only harmonics below Nyquist and the saw is reproduced cleanly.

3 · Wavetables

Summing hundreds of sin() calls per sample (a 20 Hz saw at 48 kHz needs ~1200 harmonics) is far too slow for real time. The fix: compute one cycle of the waveform once, store it in a table, and just look it up afterwards. An oscillator becomes “scale the phase to a table index, read the table.”

A sine wavetable

Sample one period of sine into N slots. Original C:

const N = 512;
var tbl: [N]f32 = undefined;

fn tblInit() void {
    const step = 2.0 * std.math.pi / @as(f32, N);
    for (&tbl, 0..) |*s, i| {
        s.* = std.math.sin(@as(f32, @floatFromInt(i)) * step);
    }
}

#define N 512
static sample_t tbl[N];

static void
tbl_init(void)
{
  float step = 2*PI_F / N;     // divide [0, 2π) into N steps
  for (size_t i = 0; i < N; ++i)
    tbl[i] = sinf(i*step);
}

Zig note — for (&tbl, 0..) |*s, i|. &tbl iterates the array by pointer so s.* can write each slot; 0.. gives the index i alongside. @as(f32, N) turns the comptime integer N into a float for the math. (N is comptime-known, so [N]f32 is a fixed-size array — no allocator needed.)

Looking up with interpolation

A phase rarely lands exactly on a stored slot. Truncating to the nearest slot sounds gritty; instead, linearly interpolate between the two neighbours. If the phase maps to index 2.3, blend 70 % of slot 2 with 30 % of slot 3:

fn tblEval(table: []const f32, phase: f32) f32 {
    const len = table.len;
    const fpos = phase * @as(f32, @floatFromInt(len));
    const i: usize = @intFromFloat(fpos); // floor (phase >= 0)
    const fr = fpos - @as(f32, @floatFromInt(i)); // fractional part
    const x0 = table[i % len];
    const x1 = table[(i + 1) % len];
    return (1.0 - fr) * x0 + fr * x1; // linear interpolation
}

static sample_t
tbl_eval(float phs)
{
  sample_t x0, x1;
  size_t i;
  float fr;

  fr = phs * N;
  i = (size_t) fr;     // integer part (floor, since phs >= 0)
  fr = fr - i;         // fractional part, 0..1

  x0 = tbl[i];
  x1 = tbl[(i+1) % N];
  return (1-fr)*x0 + fr*x1;   // linear interpolation
}

Math note — the interpolation weight. fr is “how far between the two slots are we,” from 0 to 1. The output (1-fr)*x0 + fr*x1 is a weighted average: at fr=0 you get x0, at fr=1 you get x1, and in between a straight line. It is the same fractional-blend trick we will reuse for fractional delay in chapter 7.

Zig note — @intFromFloat is the floor here. Truncating a non-negative float toward zero is the floor, and it is faster than @floor + cast. The % len keeps the “next” index in range when we wrap past the end. If len is a power of two, (i + 1) & (len - 1) is the classic faster form — but only then.

4 · A bandlimited sawtooth

Take 0 — a fixed number of harmonics

Build the table by summing harmonics instead of reading sin of the phase. Use the first 16 harmonics with the saw’s coefficients. We drop the constant $\frac{2}{\pi}$ because we normalize afterwards, so the coefficient is just $(-1)^{k+1}/k$ :

fn coeff(k: usize) f32 {
    const kf: f32 = @floatFromInt(k);
    return if (k % 2 == 0) -1.0 / kf else 1.0 / kf;
}

static float
coeff(size_t k)
{
  if (k % 2 == 0)
    return -1.f/k;
  else
    return 1.f/k;
}

fn sawTableInit(table: []f32, num_harmonics: usize) void {
    const step = 2.0 * std.math.pi / @as(f32, @floatFromInt(table.len));
    var max: f32 = 0;
    for (table, 0..) |*s, i| {
        const phi = @as(f32, @floatFromInt(i)) * step;
        var acc: f32 = 0;
        var k: usize = 1;
        while (k <= num_harmonics) : (k += 1) {
            acc += coeff(k) * std.math.sin(@as(f32, @floatFromInt(k)) * phi);
        }
        s.* = acc;
        if (@abs(acc) > max) max = @abs(acc);
    }
    if (max > 0) for (table) |*s| {
        s.* /= max;
    };
}

static void
tbl_init(void)
{
  float step = 2*PI_F / N;
  sample_t max = 0;

  for (size_t i = 0; i < N; ++i) {
    float phi = i * step;
    tbl[i] = 0;
    for (int k = 1; k <= 16; ++k)        // add up harmonics
      tbl[i] += coeff(k) * sinf(k*phi);
    if (fabsf(tbl[i]) > max) max = fabsf(tbl[i]);
  }
  if (max > 0)                            // normalize to [-1, 1]
    for (size_t i = 0; i < N; ++i) tbl[i] /= max;
}

Math note — the Gibbs phenomenon (why we normalize). A truncated Fourier sum overshoots near the waveform’s jump by about 9 % of the jump height — so the summed saw briefly exceeds ±1 (e.g. saw(0.47) ≈ 1.12). Writing that straight out would clip. So we scan for the peak and divide the whole table by it, pulling everything back into $[-1, 1]$ . The overshoot is fundamental to discontinuous waves (the series does not converge uniformly), not a bug.

Take 1 — one table per octave

Sixteen harmonics is wrong for every frequency: too few for a deep bass note, too many (aliasing again) for a high one. The classic fix is a set of tables, one per octave, each holding as many harmonics as fit below Nyquist for that octave’s top frequency. Pick the table by frequency at playback. C structure:

The Zig version, holding the whole set in one allocation:

const TBL_SIZE = 4096;
const MIN_FREQ = 20.0;

const WaveTableSet = struct {
    tables: [][TBL_SIZE]f32,
    max_freq: []f32,
    alloc: std.mem.Allocator,

    fn init(alloc: std.mem.Allocator, sr: f32) !WaveTableSet {
        const nyquist = sr / 2.0;
        // one table per octave from MIN_FREQ up to Nyquist
        const n: usize = @intFromFloat(std.math.log2(nyquist / MIN_FREQ) + 1);
        const set: WaveTableSet = .{
            .tables = try alloc.alloc([TBL_SIZE]f32, n),
            .max_freq = try alloc.alloc(f32, n),
            .alloc = alloc,
        };
        var freq: f32 = MIN_FREQ;
        for (set.tables, set.max_freq) |*tbl, *mf| {
            mf.* = freq;
            genBuf(tbl, freq, nyquist);
            freq *= 2.0; // next octave
        }
        return set;
    }

    fn deinit(self: *WaveTableSet) void {
        self.alloc.free(self.tables);
        self.alloc.free(self.max_freq);
    }

    fn genBuf(buf: *[TBL_SIZE]f32, freq: f32, nyquist: f32) void {
        const step = 2.0 * std.math.pi / @as(f32, TBL_SIZE);
        var max: f32 = 0;
        for (buf, 0..) |*s, i| {
            const phi = @as(f32, @floatFromInt(i)) * step;
            var acc: f32 = 0;
            var k: usize = 1;
            while (@as(f32, @floatFromInt(k)) * freq < nyquist) : (k += 1) {
                acc += coeff(k) * std.math.sin(@as(f32, @floatFromInt(k)) * phi);
            }
            s.* = acc;
            if (@abs(acc) > max) max = @abs(acc);
        }
        if (max > 0) for (buf) |*s| {
            s.* /= max;
        };
    }

    // first table whose max_freq exceeds the requested frequency
    fn index(self: WaveTableSet, freq: f32) usize {
        for (self.max_freq, 0..) |mf, i| {
            if (freq < mf) return i;
        }
        return self.max_freq.len - 1;
    }
};

typedef struct wave_tbl {
  float max_freq;            // valid for freq < max_freq
  sample_t buf[TBL_SIZE];
} wave_tbl;

static void
tbl_genbuf(sample_t *buf, float freq)
{
  float nyquist = sr/2.f;
  float step = 2*PI_F / TBL_SIZE;
  sample_t max = 0;
  for (size_t i = 0; i < TBL_SIZE; ++i) {
    float phi = i * step;
    buf[i] = 0;
    for (int k = 1; k*freq < nyquist; k++)   // only harmonics below Nyquist
      buf[i] += coeff(k)*sinf(k*phi);
    if (fabsf(buf[i]) > max) max = fabsf(buf[i]);
  }
  if (max > 0)
    for (size_t i = 0; i < TBL_SIZE; ++i) buf[i] /= max;
}

Math note — why one table per octave. Each octave doubles the fundamental, which halves how many harmonics fit under Nyquist. A table built for the top of its octave has no aliasing for any frequency below that — so doubling the boundary each step (20, 40, 80, …) covers the whole range with about $\log_2(\text{Nyquist}/20) \approx 11$ tables. The count is a logarithm because we are doubling.

Zig note — slices and ownership. try alloc.alloc([TBL_SIZE]f32, n) asks the allocator for n tables and returns a slice ([][TBL_SIZE]f32); the matching deinit frees it. for (set.tables, set.max_freq) |*tbl, *mf| walks two slices in lockstep — Zig checks they are the same length. try propagates an out-of-memory error to the caller.

Take 2 — crossfade between tables

Switching tables abruptly at an octave boundary clicks. Blend the chosen table with the next one, weighted by how far into the octave the frequency is:

fn evalBandlimited(set: WaveTableSet, freq: f32, phase: f32) f32 {
    const idx = set.index(freq);
    if (idx >= set.tables.len - 1)
        return tblEval(&set.tables[idx], phase); // top table, nothing to blend
    const f1 = set.max_freq[idx];
    const f0 = if (idx == 0) 0.0 else set.max_freq[idx - 1];
    const w = (freq - f0) / (f1 - f0);
    return tblEval(&set.tables[idx], phase) * (1.0 - w) +
        tblEval(&set.tables[idx + 1], phase) * w;
}

float f1 = tbl[idx].max_freq;
float f0 = idx==0 ? 0 : tbl[idx-1].max_freq;
float w = (freq[i] - f0) / (f1 - f0);
out[i] = tbl_eval(tbl[idx].buf,   phs[i]) * (1-w) +
         tbl_eval(tbl[idx+1].buf, phs[i]) * w;

Math note — same blend, different axis. This is the interpolation trick again, but now we blend across tables (frequency) instead of across samples (phase). w measures how far the frequency sits between the previous and current octave boundary; at the boundaries w is 0 or 1 so there is no jump, and in between the timbre morphs smoothly.

Driving a Phasor (chapter 3) through evalBandlimited gives a saw that stays clean from the lowest bass to the top of the spectrum.

5 · Improvements (pointers)

Table selection by binary search instead of linear scan.
Variable table size — high octaves need few harmonics, so they can be shorter, saving memory.
Faster generation — build the spectrum directly and take an inverse FFT: $O(NM) \to O(N\log N)$ .
Fixed-point phase index — a f32 index shares its 24-bit mantissa between integer and fraction; a Q12.24 fixed-point index gives the fraction far more precision for big tables.
Other shapes for free — a square is two phase-shifted saws subtracted; changing the shift changes the pulse width.

Exercises

Generate the single-table saw with 4, 16, then 64 harmonics and compare brightness. More harmonics = brighter, until they alias.
Build a bandlimited square by subtracting a phase-shifted saw from a saw. Vary the shift for pulse-width modulation.
Replace linear interpolation in tblEval with nearest-sample lookup and listen to the gritty zipper noise — that is why we interpolate.