Audio DSP & the FFT

01Two domains

Looking at a waveform tells you how loud a sound is moment to moment, but not whether it's a bass rumble or a hiss. For that you need a second view of the same samples. A sampled signal has two equivalent representations: a sequence of amplitude samples over time (the time domain), and a set of sinusoidal components, amplitude and phase at each frequency, that sum to it (the frequency domain). The DFT is the exact, invertible bridge between them. Equalizers are designed, and spectral effects and analyzers are built, in terms of the frequency domain.

02The DFT & the FFT

An analyzer has to do this for every block, many times a second, so the cost of the transform matters. The DFT of N samples produces N complex bins; each bin's magnitude = √(re²+im²) is proportional to that frequency's amplitude and its phase = atan2(im, re) is its offset^[3]. Computed directly it's O(N²). The FFT (Cooley-Tukey) computes the same result in O(N log N) by recursively splitting the transform into even/odd halves combined with twiddle-factor butterflies^[1][2].

Each bin k is the signal compared against a spinning arrow at that bin's frequency:

In words: for bin k, multiply every sample by an arrow that spins k full turns over the block and add up the results. If the signal contains a frequency that keeps pace with the arrow, the products line up and the sum is large; anything else cancels out. Doing that for all N bins directly is N × N multiplications, which is what the FFT avoids. Hover any term to see what it stands for.

The widget computes the magnitude spectrum of a test signal every frame. The square wave shows the odd-harmonic series:

// In-place radix-2 decimation-in-time. N MUST be a power of two.
// Needs <complex>, <cmath>, <utility>, <vector>.
void fftRadix2(std::vector<std::complex<float>>& samples) {
    const size_t n = samples.size();
    for (size_t i = 1, j = 0; i < n; ++i) {        // bit-reversal permutation
        size_t bit = n >> 1;
        for (; j & bit; bit >>= 1) j ^= bit;
        j ^= bit; if (i < j) std::swap(samples[i], samples[j]);
    }
    for (size_t len = 2; len <= n; len <<= 1) {          // stages: 2, 4, 8, ..., N
        float angle = -2.0f * 3.14159265f / (float)len;   // -2pi/len (forward transform)
        std::complex<float> twiddleStep(std::cos(angle), std::sin(angle));
        for (size_t blockStart = 0; blockStart < n; blockStart += len) {
            std::complex<float> twiddle(1, 0);         // e^(-i2πk/len), advanced by one step per k
            for (size_t k = 0; k < len/2; ++k) {           // the butterfly
                auto even = samples[blockStart+k], odd = samples[blockStart+k+len/2] * twiddle;
                samples[blockStart+k] = even + odd; samples[blockStart+k+len/2] = even - odd;
                twiddle *= twiddleStep;
            }
        }
    }
}

use num_complex::Complex32;
// In-place radix-2 DIT. n MUST be a power of two. Production: the `rustfft` crate.
fn fft_radix2(samples: &mut [Complex32]) {
    let n = samples.len();
    let mut j = 0usize;
    for i in 1..n {                                  // bit-reversal permutation
        let mut bit = n >> 1;
        while j & bit != 0 { j ^= bit; bit >>= 1; }
        j ^= bit; if i < j { samples.swap(i, j); }
    }
    let mut len = 2;
    while len <= n {                                // stages
        let angle = -2.0 * std::f32::consts::PI / len as f32;   // -2pi/len (forward transform)
        let twiddle_step = Complex32::new(angle.cos(), angle.sin());
        let mut block_start = 0;
        while block_start < n {
            let mut twiddle = Complex32::new(1.0, 0.0);   // e^(-i2πk/len), advanced by one step per k
            for k in 0..len/2 {                       // the butterfly
                let (even, odd) = (samples[block_start+k], samples[block_start+k+len/2] * twiddle);
                samples[block_start+k] = even + odd; samples[block_start+k+len/2] = even - odd;
                twiddle *= twiddle_step;
            }
            block_start += len;
        }
        len <<= 1;
    }
}

03Windowing & leakage

Your analyzer shows a single clean tone as a wide smear covering a dozen bins, hiding a quieter tone nearby. The cause is the block boundary. Pulling a finite N-sample block from a longer signal implicitly multiplies by a rectangular window. Unless the signal is exactly periodic in the block, the edges are discontinuous, and that smears each frequency's energy across many bins: spectral leakage. The rectangular window's first side lobe is only 13 dB down, so leakage is severe. A window function (Hann, Hamming, Blackman) tapers the block toward zero at the edges, applied before the FFT, cutting the side lobes^[4].

The widget analyzes a single sine with each window:

04The STFT & spectrogram

A music visualizer has to show the bass hitting on each beat as it happens. A single FFT describes one block of time. The Short-Time Fourier Transform slides a windowed FFT along the signal, transforming each chunk into a column of a 2D time × frequency map; plot magnitude (in dB) as brightness and you have a spectrogram. The hop size is how far the window advances between frames; typical analysis hops are 25 to 50% of the window, so consecutive windows overlap. Without overlap, a tapered window all but ignores transients that fall near a frame's edges.

05Digital filters

A sound behind a wall loses its high frequencies, so muffling it means turning the highs down smoothly while the lows pass through. That's a filter. A linear time-invariant digital filter of finite order is a difference equation, each output sample a weighted sum of recent inputs and recent outputs:

In words: read it as a mixing desk on a sliding window. The new output y[n] is a weighted blend of the current and past inputs (the b terms, the feed-forward path) minus a weighted blend of past outputs (the a terms, the feedback path)^[7]. A subscript like b₁ just labels a coefficient; [n−1] means "one sample ago". With no feedback (no a terms) it is an FIR filter; with feedback it is IIR. Hover any symbol to see what it stands for.

• FIR (no feedback): always stable and can be exactly linear-phase (no phase distortion), but a sharp response needs many taps, so more compute.
• IIR (has feedback): achieves a given response with far fewer coefficients (cheap), but can be unstable and has nonlinear phase.

The biquad is a second-order IIR section (two poles, two zeros), the standard building block of audio EQ. The RBJ "Audio EQ Cookbook" gives ready coefficients for low/high/band-pass, notch, peak, and shelf filters^[6]. Conceptually, zeros pull the response down, poles push it up (resonance), and a pole on or outside the unit circle is unstable.

struct Biquad {
    float b0,b1,b2,a1,a2;                 // coefficients (a0 normalized to 1)
    float x1=0,x2=0,y1=0,y2=0;            // state: last two inputs/outputs
    void setLowpass(float cutoffHz, float q, float sampleRate) {
        float w0 = 2.0f*3.14159265f*cutoffHz/sampleRate;
        float cosW0 = std::cos(w0), alpha = std::sin(w0)/(2.0f*q);   // RBJ's cos(w0) and alpha
        float a0 = 1.0f + alpha;                // normalize everything by a0
        b0 = ((1-cosW0)*0.5f)/a0; b1 = (1-cosW0)/a0; b2 = ((1-cosW0)*0.5f)/a0;
        a1 = (-2*cosW0)/a0; a2 = (1-alpha)/a0;
    }
    void processBlock(float* samples, int count) {
        for (int i = 0; i < count; ++i) {
            float x0 = samples[i];            // x0 = x[n], x1 = x[n-1], y1 = y[n-1] ...
            float y0 = b0*x0 + b1*x1 + b2*x2 - a1*y1 - a2*y2;   // the difference equation
            x2=x1; x1=x0; y2=y1; y1=y0;       // shift history (y feedback is serial)
            samples[i] = y0;
        }
    }
};

struct Biquad { b0:f32,b1:f32,b2:f32,a1:f32,a2:f32, x1:f32,x2:f32,y1:f32,y2:f32 }
impl Biquad {
    fn set_lowpass(&mut self, cutoff_hz: f32, q: f32, sample_rate: f32) {
        let w0 = 2.0*std::f32::consts::PI*cutoff_hz/sample_rate;
        let (cos_w0, alpha) = (w0.cos(), w0.sin()/(2.0*q));  // RBJ's cos(w0) and alpha
        let a0 = 1.0 + alpha;                 // normalize by a0
        self.b0 = ((1.0-cos_w0)*0.5)/a0; self.b1 = (1.0-cos_w0)/a0; self.b2 = ((1.0-cos_w0)*0.5)/a0;
        self.a1 = (-2.0*cos_w0)/a0; self.a2 = (1.0-alpha)/a0;
    }
    fn process_block(&mut self, samples: &mut [f32]) {
        for sample in samples.iter_mut() {
            let x0 = *sample;                  // x0 = x[n], x1 = x[n-1], y1 = y[n-1] ...
            let y0 = self.b0*x0 + self.b1*self.x1 + self.b2*self.x2 - self.a1*self.y1 - self.a2*self.y2;
            self.x2=self.x1; self.x1=x0; self.y2=self.y1; self.y1=y0;   // shift history (y feedback is serial)
            *sample = y0;
        }
    }
}

The widget plots the biquad's frequency response, recomputed from the RBJ coefficients:

06Convolution reverb

A footstep recorded in a dry studio sounds wrong in a cathedral level, and re-recording every sound in every space isn't an option. Convolution reverb convolves the dry signal with a measured room impulse response (the recorded response of a real space to an impulse); the output sounds as if played there. Direct convolution is O(N·M), far too slow for a seconds-long IR. The convolution theorem (time-domain convolution equals frequency-domain multiplication) is why FFT convolution is fast: zero-pad both to at least N+M−1 so the FFT's circular convolution equals the linear one, FFT both, multiply the spectra, inverse-FFT, O(N log N)^[1].

07Real-time game DSP

Sixty-four voices each need a distance filter, and the whole mix needs reverb, inside a callback with a few milliseconds to spare. All of this runs inside the hard-real-time audio callback (or worker threads feeding it): no heap alloc, no locks, no I/O, no unbounded work. Filters and reverb are pre-allocated with fixed state; parameter changes arrive via the lock-free ring, not a mutex. Process whole blocks per effect (cache- and SIMD-friendly).

Wrong answers, and why: the FFT is exact (not approximate) and real input gives N/2+1 bins (not N); and windowing trades leakage for resolution (no best window) and is applied before the FFT.

08Pitfalls

09What's next

Next in the series is the first milestone, the 2D capstone, which wires the loop, the sprite renderer, input, asset loading and this audio engine into a small playable game. The Audio Engine module this builds on covers the callback and mixer; the full path is on the series hub.

1. Steven W. Smith. The Scientist and Engineer's Guide to Digital Signal Processing, ch. 12 (the FFT), 9 (the real DFT, N/2+1), 18 (FFT convolution). dspguide.com. The FFT is exact and more precise; bin counts; the convolution theorem.
2. James Cooley and John Tukey. "An Algorithm for the Machine Calculation of Complex Fourier Series." Math. of Computation 19(90), 1965. overview. The O(N log N) divide-and-conquer FFT.
3. Julius O. Smith III. Mathematics of the DFT (with Audio Applications). CCRMA, Stanford. ccrma.stanford.edu. The DFT definition, magnitude and phase, for audio.
4. Julius O. Smith III. Spectral Audio Signal Processing, "Spectrum Analysis Windows." dsprelated.com. Window side-lobe/main-lobe figures (rect −13, Hann −31, Hamming −43, Blackman −58 dB) and the tradeoff.
5. Fredric J. Harris. "On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform." Proc. IEEE 66(1), 1978. doi.org. The canonical window catalog and performance parameters, including the DFT-even (periodic) form for spectral analysis.
6. Robert Bristow-Johnson. "Cookbook formulae for audio EQ biquad filter coefficients" (Audio EQ Cookbook). w3.org/TR/audio-eq-cookbook. The RBJ biquad coefficients for every filter type; Direct Form I recommendation.
7. Julius O. Smith III. Introduction to Digital Filters (with Audio Applications). CCRMA. ccrma.stanford.edu. FIR vs IIR, poles and zeros, stability, the difference equation.
8. Julius O. Smith III. Physical Audio Signal Processing, "Freeverb" / "Artificial Reverberation." ccrma.stanford.edu. Freeverb: eight parallel Schroeder-Moorer filtered-feedback combs into four series Schroeder all-passes per channel.
9. Matteo Frigo and Steven Johnson. FFTW ("Fastest Fourier Transform in the West"). fftw.org. The production FFT: any size, O(N log N), Cooley-Tukey + Rader/Bluestein. (KissFFT, rustfft are lighter alternatives.)
10. Firelight Technologies. FMOD Engine User Manual 2.03, "Core API: Using DSP Effects." fmod.com. The directed DSP graph from generators to the head node; effects on channels (per voice) or channel groups (per bus).
11. William G. Gardner. "Efficient Convolution without Input-Output Delay." JAES 43(3), 1995. aes.org. Non-uniformly partitioned FFT convolution plus a direct-form head for zero-latency convolution reverb.