3.4 Wavetable synthesis

3.4 Wavetable synthesis#

An algorithmic perspective#

We’ve now established additive synthesis as a powerful technique grounded in the Fourier series. But how efficient is it computationally? Let’s think about this from a computer science perspective.

First, recall the formula for additive synthesis, omitting \(\phi_k\) as the initial phases are perceptually negligible:

\[x(t) = \sum_{k=1}^{K} a_k \sin(2\pi [k \cdot f_0] \, t).\]

To synthesize additive synthesis in the digital domain, we can follow the recipe for synthesis from Chapter 2, which tells us that individual samples \(x[n] = x(n / f_s)\). Accordingly, to synthesize \(N\) samples of a tone with \(K\) harmonics, we compute:

\[x[n] = \sum_{k=1}^{K} a_k \sin(2\pi [k \cdot f_0] [n / f_s]) \quad \text{for } n = 0, 1, \ldots, N-1.\]

This requires \(K\) calls to sin per sample, or \(K \cdot N\) total evaluations. For one second of audio at \(f_s = 44{,}100\) with \(K = 32\) harmonics, that’s \(32 \times 44{,}100 \approx 1.4\) million sin evaluations. Modern computers may be able to keep up with this demand in real time, but what if you need to run many synthesizers in parallel (e.g., in a DAW), or what if you’re working on hardware with limited compute?

What structure can we exploit to make additive synthesis more efficient? Here’s the key insight: the output of additive synthesis is periodic. The waveform repeats every \(t_0 = 1/f_0\) seconds, or equivalently every \(f_s / f_0\) samples. So we only need to compute one cycle of the waveform, then cache the result and repeat it.

Building the wavetable#

This insight leads to wavetable synthesis. The first step is to compute a a wavetable: a single cycle of the waveform stored as an array of \(M\) samples:

\[\texttt{table}[m] = \sum_{k=1}^{K} a_k \sin\!\left(2\pi k \cdot \frac{m}{M}\right) \quad \text{for } m = 0, 1, \ldots, M - 1.\]

Here \(m / M\) maps the table index \(m \in \{0, \ldots, M - 1\}\) to the range \([0, 1)\), covering exactly one period. In code:

def build_wavetable(a: list[float], M: int = 2048) -> np.ndarray:
    K = len(a)
    a = np.array(a)
    k = 1 + np.arange(K)
    m = np.arange(M)
    # Broadcasting: (M, 1) * (K,) -> (M, K)
    table = (a * np.sin(2 * np.pi * k * m[:, np.newaxis] / M)).sum(axis=1)
    return table

Building the table costs \(O(K \cdot M)\) operations, but it runs only once for a given waveform shape.

Reading the wavetable#

To produce output at frequency \(f_0\), we need to “repeat” (read from) the table at the right rate. The table spans one period, and we want the output to complete \(f_0\) cycles per second. Working from the units:

The table has \(M\) \(\frac{\text{indices}}{\text{cycle}}\).
We want \(f_0\) \(\frac{\text{cycles}}{\text{second}}\).
The output sample rate is \(f_s\) \(\frac{\text{samples}}{\text{second}}\).

The phase increment — how far we advance through the table per output sample — is therefore:

\[\Delta m = f_0 \cdot \frac{M}{f_s} \quad \left[\frac{\text{table indices}}{\text{output sample}}\right].\]

After \(n\) output samples, we’ve accumulated \(\tilde{m} = n \cdot \Delta m\) table indices. To read the table, we wrap this phase modulo \(M\):

\[x[n] = \texttt{table}\!\left[\; \tilde{m} \bmod M \;\right].\]

A single-cycle wavetable of 256 samples (top) and its repetition over four cycles (bottom) — Fig. 6 Top: a single-cycle wavetable of \(M = 256\) indices. Bottom: the output signal produced by repeating this table at a different rate. The dashed lines mark cycle boundaries.#

However, there’s a big problem: \(\tilde{m}\) is not guaranteed to be an integer! How might we “read from” the table at a fractional index?

The animation below shows the process in motion. The gold pointer advances by \(\Delta m\) table indices per output sample and wraps modulo \(M\); because \(\Delta m\) is not an integer it usually lands between entries. The animation shows two different strategies for resolving this issue (truncating lookup and linear interpolation), which are discussed in detail below.

Truncating lookup#

The simplest implementation uses the floor operation to truncate the fractional part of \(\tilde{m} = n \cdot \Delta m\) before lookup:

\[x[n] = \texttt{table}\!\left[\lfloor \tilde{m} \rfloor \bmod M\right]\]

# hide
import numpy as np
import pyquist as pq


def build_wavetable(a: list[float], M: int = 2048) -> np.ndarray:
    K = len(a)
    a = np.array(a)
    k = 1 + np.arange(K)
    m = np.arange(M)
    table = (a * np.sin(2 * np.pi * k * m[:, np.newaxis] / M)).sum(axis=1)
    return table

def wavetable_truncate(
    table: np.ndarray, f_0: float, f_s: int, N: int,
) -> pq.Audio:
    M = len(table)
    delta_m = f_0 * M / f_s        # phase increment (table indices per sample)
    m_tilde = np.arange(N) * delta_m
    indices = m_tilde.astype(int) % M    # truncate to nearest table entry
    return pq.Audio(table[indices], f_s)

# Read a coarse M = 8 sine table back at 440 Hz. Try changing M or f_0 and listen!
table = build_wavetable([1.0], M=8)
audio = wavetable_truncate(table, f_0=440.0, f_s=44100, N=44100)
pq.play(audio)

This works, but when \(\Delta m\) is not an integer (which is common — e.g., \(f_0 = 440\), \(M = 2048\), \(f_s = 44{,}100\) gives \(\Delta m \approx 20.43\)), this is a noisy operation that introduces error into our approximation. The effect is especially audible with a small table. Compare an exact sine wave to truncating wavetable lookup with \(M = 8\):

Exact sine wave at 440 Hz

Truncating wavetable, M = 8

Exact vs. truncating wavetable sine at 440 Hz. The coarse table produces audible stepping artifacts.

Linear interpolation#

A better approach is to interpolate between adjacent table entries, assuming the signal varies linearly between them. Given a fractional index \(\tilde{m} = n \cdot \Delta m\), we split it into an integer part \(\lfloor \tilde{m} \rfloor\) and a fractional part \(\alpha = \tilde{m} - \lfloor \tilde{m} \rfloor\), then blend:

\[x[n] = (1 - \alpha) \cdot \texttt{table}\!\left[\lfloor \tilde{m} \rfloor \bmod M\right] + \alpha \cdot \texttt{table}\!\left[(\lfloor \tilde{m} \rfloor + 1) \bmod M\right].\]

In code:

# hide
import numpy as np
import pyquist as pq


def build_wavetable(a: list[float], M: int = 2048) -> np.ndarray:
    K = len(a)
    a = np.array(a)
    k = 1 + np.arange(K)
    m = np.arange(M)
    table = (a * np.sin(2 * np.pi * k * m[:, np.newaxis] / M)).sum(axis=1)
    return table

def wavetable_interp(
    table: np.ndarray, f_0: float, f_s: int, N: int,
) -> pq.Audio:
    M = len(table)
    delta_m = f_0 * M / f_s        # phase increment (table indices per sample)
    m_tilde = np.arange(N) * delta_m   # fractional table positions
    m_floor = m_tilde.astype(int)
    alpha = m_tilde - m_floor
    x = (1 - alpha) * table[m_floor % M] + alpha * table[(m_floor + 1) % M]
    return pq.Audio(x, sample_rate=f_s)

# The same coarse M = 8 sine table, now read with linear interpolation.
table = build_wavetable([1.0], M=8)
audio = wavetable_interp(table, f_0=440.0, f_s=44100, N=44100)
pq.play(audio)

Linear interpolation adds negligible computational cost (a multiply and an add per sample) but dramatically reduces the error, especially when \(M\) is small. Compare the same \(M = 8\) table with interpolation:

Exact sine wave at 440 Hz

Interpolated wavetable, M = 8

Exact vs. linear interpolation wavetable sine at 440 Hz. Even with only 8 table entries, interpolation produces a much smoother result.

In practice, most wavetable synthesizers use at least linear interpolation; some use higher-order schemes (cubic, sinc) for even better quality.

Sawtooth wave at 220 Hz synthesized via wavetable lookup with linear interpolation, using the same \(K = 32\) harmonic recipe.

Complexity#

Direct additive synthesis costs \(O(K \cdot N)\) to synthesize \(N\) samples: \(K\) sin evaluations per sample. Wavetable synthesis costs \(O(K \cdot M)\) to build the table once, then \(O(N)\) to read it — for a total of \(O(K \cdot M + N)\), where \(M << N\). Since \(M\) is a fixed constant (typically 2048 or 4096), the table-building step is a one-time cost that does not grow with the output length.

The key implication: once the wavetable is computed, the per-sample cost of synthesis is \(O(1)\) regardless of how many harmonics \(K\) went into the table. A 4-harmonic triangle wave and a 64-harmonic sawtooth are equally cheap to synthesize once their tables are built. This is in stark contrast to direct additive synthesis, where doubling \(K\) doubles the cost.

Tip

On a modern machine with NumPy, the wavetable version of a 32-harmonic sawtooth runs roughly 20–40x faster than the direct additive computation. The speedup grows with \(K\): the more harmonics in your recipe, the more work you avoid by precomputing the table.

The full implementation and timing comparison is in code/wavetable.py.