Audio DSP & the FFT

01Two domains

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. Frequency analysis is what lets you build an equalizer, spectral effects, and any analyzer or visualizer.

02The DFT & the FFT

The DFT of N samples produces N complex bins; each bin's magnitude = √(re²+im²) is that frequency's amplitude and its phase = atan2(im, re) its offset. 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].

Sweep a sine across frequency (and detune it between bins) or switch to a square wave; the live magnitude spectrum is below. The square wave shows the odd-harmonic series:

// In-place radix-2 decimation-in-time. N MUST be a power of two.
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 ang = -2.0f * 3.14159265f / (float)len;     // -2pi/len (forward)
        std::complex<float> step(cosf(ang), sinf(ang));
        for (size_t s = 0; s < n; s += len) {
            std::complex<float> tw(1, 0);
            for (size_t k = 0; k < len/2; ++k) {           // the butterfly
                auto even = samples[s+k], odd = samples[s+k+len/2] * tw;
                samples[s+k] = even + odd; samples[s+k+len/2] = even - odd;
                tw *= step;
            }
        }
    }
}

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 ang = -2.0 * std::f32::consts::PI / len as f32;
        let step = Complex32::new(ang.cos(), ang.sin());
        let mut s = 0;
        while s < n {
            let mut tw = Complex32::new(1.0, 0.0);
            for k in 0..len/2 {                       // the butterfly
                let (even, odd) = (samples[s+k], samples[s+k+len/2] * tw);
                samples[s+k] = even + odd; samples[s+k+len/2] = even - odd;
                tw *= step;
            }
            s += len;
        }
        len <<= 1;
    }
}

03Windowing & leakage

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].

A single sine whose frequency you can push between bins; switching windows drops the leakage as the peak widens:

04The STFT & spectrogram

One FFT is a snapshot. The Short-Time Fourier Transform slides a windowed FFT along the signal, taking a transform of each chunk, producing 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 windows overlap and you don't miss transients).

05Digital filters

A linear filter is a difference equation, each output sample a weighted sum of recent inputs and recent outputs:

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, can be exactly linear-phase (no phase distortion), but needs many taps, 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 workhorse. 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 cw = cosf(w0), alpha = sinf(w0)/(2.0f*q);
        float a0 = 1.0f + alpha;                // normalize everything by a0
        b0 = ((1-cw)*0.5f)/a0; b1 = (1-cw)/a0; b2 = ((1-cw)*0.5f)/a0;
        a1 = (-2*cw)/a0; a2 = (1-alpha)/a0;
    }
    void processBlock(float* s, int count) {
        for (int i = 0; i < count; ++i) {
            float x0 = s[i];                  // the difference equation
            float y0 = b0*x0 + b1*x1 + b2*x2 - a1*y1 - a2*y2;
            x2=x1; x1=x0; y2=y1; y1=y0;       // shift history (y feedback is serial)
            s[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 (cw, alpha) = (w0.cos(), w0.sin()/(2.0*q));
        let a0 = 1.0 + alpha;                 // normalize by a0
        self.b0 = ((1.0-cw)*0.5)/a0; self.b1 = (1.0-cw)/a0; self.b2 = ((1.0-cw)*0.5)/a0;
        self.a1 = (-2.0*cw)/a0; self.a2 = (1.0-alpha)/a0;
    }
    fn process_block(&mut self, s: &mut [f32]) {
        for v in s.iter_mut() {
            let x0 = *v;                       // difference equation
            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
            *v = y0;
        }
    }
}

The biquad's frequency response, recomputed live from the RBJ coefficients; sweep cutoff, Q, and type, and push Q high to see the resonant peak:

06Convolution reverb

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: FFT both, multiply the spectra, inverse-FFT, O(N log N)^[1].

07Real-time game DSP

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

That's the signal processing under the audio engine: see sound as frequencies, read the spectrum honestly, shape it with filters, and place it in a room. Back to the series hub for the rest of the engine, or the Audio Engine module this builds on.

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.
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. Schroeder/Moorer comb+allpass reverb; Freeverb's 8 combs + 4 all-passes.
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 DSP Architecture. fmod.com. Game audio as a directed graph of DSP effect nodes, per-voice and per-bus.