import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

6.7 Implementing FM

The integrated FM formula, \(\sin(2\pi f_c t + \tfrac{D}{f_m}\sin(2\pi f_m t))\), is easy to compute directly. But there is a more flexible and more general way to implement FM that connects the pieces we have built in this book:

1. In the previous sections we built a correct time-varying oscillator that accumulates phase, osc(freq).

2. FM is just a time-varying oscillator whose frequency signal happens to be another oscillator.

Combining these, we can write a general FM oscillator whose modulating signal can be any sound at all, not just a single sinusoid. This is far more expressive than the closed-form equation: the modulator can be a chord, a noise source, or even a recorded sample.

The interactive example below builds exactly this. First, osc(freq) is our time-varying oscillator: it takes a per-sample frequency signal (in Hz), accumulates it into phase with np.cumsum, and applies np.sin. Then fm(f_c, f_m, D) implements a general FM oscillator, where f_m can be any sound with [-1, 1] amplitude, multiplied by D so the center frequency offset is in [-D, D]. We can use this primitive to implement fm_classic, the “classic” FM synthesis equation from this chapter with constant f_m and index of modulation I = D / f_m. In fm_classic, the nested osc structure is clear. Edit f_c, f_m, and I, then listen and watch the spectrum. Try to reproduce a bright harmonic tone, then an inharmonic bell.

Here osc recomputes a sine for every sample. When the modulator is itself a simple periodic waveform, this could be made more efficient by combining it with the wavetable synthesis of Chapter 3: precompute one cycle of the sine into a table and read it back at the accumulated phase, trading the per-sample np.sin for a cheap table lookup.

# hide
import numpy as np
import pyquist as pq

def osc(freq: float | np.ndarray, T: float, f_s: int = 44100) -> pq.Audio:
    """A time-varying sine oscillator. `freq` gives the frequency (in Hz) at each
    sample, so the pitch can vary over time. We accumulate frequency into phase with
    `np.cumsum` (a running total), then take the sine of the accumulated phase."""
    N = int(T * f_s)
    if np.ndim(freq) == 0:  # a single number -> a constant frequency
        freq = np.full(N, freq)
    if len(freq) != N:
        raise ValueError(f"freq has length {len(freq)}, expected {N} (= T * f_s)")
    theta = np.cumsum(2 * np.pi * freq / f_s)  # accumulated phase, in radians
    return pq.Audio(np.sin(theta).astype(np.float32), f_s)


def fm(f_c: float, f_m: np.ndarray | pq.Audio, D: float, T: float, f_s: int = 44100) -> pq.Audio:
    """General FM synthesis, built on top of `osc`. A carrier at `f_c` Hz has its
    frequency wobbled by a modulator signal `f_m`. `D` is the peak frequency deviation
    in Hertz, so the carrier's instantaneous frequency is `f_c + D * f_m`."""
    if isinstance(f_m, pq.Audio):
        f_m = f_m.as_mono().samples[:, 0]
    return osc(f_c + D * f_m, T, f_s)


def fm_classic(f_c: float, f_m: float, I: float, T: float, f_s: int = 44100) -> pq.Audio:
    """Classic FM synthesis, built on top of two calls to `osc`. Carrier at `f_c` Hz
    and sinusoidal modulation at `f_m` Hz; `I = D / f_m` is the index of modulation."""
    D = I * f_m
    return osc(f_c + D * osc(f_m, T), T, f_s)
    #return fm(f_c, osc(f_m, T), I * f_m, T)   # equivalent


T = 4.0
f_c = 440.0
f_m = 4.0
I = 20.0

basic_osc = osc(f_c, T) # basic sinusoid at f_c
fm_osc = fm(f_c, osc(f_m, T), I * f_m, T)  # classic FM with parameters f_c/f_m/I
#fm_osc = fm_classic(f_c, f_m, I, T)   # equivalent

pq.play(basic_osc)
pq.plot_spec(basic_osc)
pq.play(fm_osc)
pq.plot_spec(fm_osc)

<Axes: xlabel='Time (s)', ylabel='Frequency (Hz)'>