Short-time Fourier transform and spectrogram

Deep dive

How it works

In: signal x, window length N, hop, window function · Out: matrix X[k, t] of complex bins (N/2+1 per frame for real input) · State: none for offline use; one frame buffer when streaming · Cost: O(N log N) per frame, so about (L/hop)Nlog2(N) for L samples

x:  |---------------------------------------------------|
    [ frame 0: N ]
         [ frame 1 ]      shifted by hop
              [ frame 2 ]
each frame: x[t*hop + i] * w[i]  ->  FFT  ->  column t of X[k, t]
spectrogram = 10*log10(|X|^2 + eps)        (dB, floored)

1. Pick N (power of two for radix-2), hop and window. Periodic Hann (0.5 - 0.5 cos(2 pi i / N)) is the default.
2. For frame t: multiply x[t*hop .. t*hop+N] by w, run the FFT, keep bins 0..=N/2 (real input is conjugate-symmetric).
3. Power or magnitude: |X|^2 or |X|; take a log for display and most ML features. Add a floor (eps) so silence is not -inf.
4. Optional: pool bins onto a mel or other filterbank (a sparse matrix product per frame).
5. Inverse (weighted overlap-add): inverse-FFT each frame, multiply by the synthesis window, add at t*hop, divide by sum_t w^2[n - t*hop]. Reconstruction is exact where that sum is non-zero.
6. COLA: if the windows alone sum to a constant (sum_t w[n - t*hop] = const), plain overlap-add of the frames reconstructs. Periodic Hann is COLA at hop N/2 and N/4. Dividing by sum w^2, as above, makes reconstruction exact for any hop with non-zero window energy, and is the right normalization after you modify the spectrum.

Code

use std::f32::consts::PI;

/// In-place radix-2 complex FFT (len power of two). `inverse` omits the 1/n scale.
fn fft(re: &mut [f32], im: &mut [f32], inverse: bool) {
    let n = re.len();
    let mut j = 0;
    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 { re.swap(i, j); im.swap(i, j); }
    }
    let sign = if inverse { 1.0 } else { -1.0 };
    let mut len = 2;
    while len <= n {
        let ang = sign * 2.0 * PI / len as f32;
        for start in (0..n).step_by(len) {
            for k in 0..len / 2 {
                let (s, c) = (ang * k as f32).sin_cos();   // precompute a table in real code
                let (a, b) = (start + k, start + k + len / 2);
                let tr = re[b] * c - im[b] * s;
                let ti = re[b] * s + im[b] * c;
                re[b] = re[a] - tr; im[b] = im[a] - ti;
                re[a] += tr;        im[a] += ti;
            }
        }
        len <<= 1;
    }
}

/// Periodic Hann window: satisfies COLA at hop = n/2, n/4, ...
pub fn hann(n: usize) -> Vec<f32> {
    (0..n).map(|i| 0.5 - 0.5 * (2.0 * PI * i as f32 / n as f32).cos()).collect()
}

/// Frames of n/2+1 complex bins (re, im). No padding: frame t covers x[t*hop .. t*hop+n].
pub fn stft(x: &[f32], n: usize, hop: usize) -> Vec<Vec<(f32, f32)>> {
    let w = hann(n);
    let mut out = Vec::new();
    let mut t = 0;
    while t + n <= x.len() {
        let mut re: Vec<f32> = (0..n).map(|i| x[t + i] * w[i]).collect();
        let mut im = vec![0.0; n];
        fft(&mut re, &mut im, false);
        out.push((0..=n / 2).map(|k| (re[k], im[k])).collect());
        t += hop;
    }
    out
}

/// Inverse by weighted overlap-add: sum(w*frame) / sum(w^2). Exact where the window sum is non-zero.
pub fn istft(frames: &[Vec<(f32, f32)>], n: usize, hop: usize, len: usize) -> Vec<f32> {
    let w = hann(n);
    let (mut y, mut norm) = (vec![0.0f32; len], vec![0.0f32; len]);
    for (f, spec) in frames.iter().enumerate() {
        let mut re = vec![0.0; n];
        let mut im = vec![0.0; n];
        for k in 0..=n / 2 {                       // rebuild the conjugate-symmetric half
            re[k] = spec[k].0; im[k] = spec[k].1;
            if k > 0 && k < n / 2 { re[n - k] = spec[k].0; im[n - k] = -spec[k].1; }
        }
        fft(&mut re, &mut im, true);
        for i in 0..n {
            y[f * hop + i] += re[i] / n as f32 * w[i];
            norm[f * hop + i] += w[i] * w[i];
        }
    }
    for i in 0..len { if norm[i] > 1e-6 { y[i] /= norm[i]; } }
    y
}

/// Power spectrogram in dB, floored: 10*log10(|X|^2 + eps).
pub fn power_db(frames: &[Vec<(f32, f32)>]) -> Vec<Vec<f32>> {
    frames.iter().map(|f| f.iter().map(|&(r, i)| 10.0 * (r * r + i * i + 1e-10).log10()).collect()).collect()
}

Tests (run): a 1 kHz tone at fs 8 kHz, N=512, hop 128 peaks in bin 64 in every frame with magnitude 0.25*N (A * sum(w) / 2); random noise reconstructs to within 1e-3 at hop N/4 away from the first/last N samples; a bin-64 sine gives 10*log10((N/4)^2) dB. The FFT recomputes sin_cos in the inner loop for brevity; production code uses a twiddle table.

Complexity and constants

• Per frame: a complex FFT costs order of 5 N log2 N flops; a real FFT about half. N=512: on the order of 11,000-23,000 flops, microseconds per frame.
• Frames per second of audio: fs/hop (16 kHz, hop 128: 125 frames/s). Output per second at N=512: 125 x 257 bins = 32 k complex values, 257 kB as complex f32 or 128 kB as f32 magnitude.
• Resolution: bin spacing fs/N, effective resolution of Hann about 1.5 bins (noise bandwidth); time resolution about the window length N/fs.
• Redundancy: N+2 real numbers (N/2+1 complex) per hop input samples, about N/hop: 2x at hop N/2, 4x at hop N/4.

Efficiency

The work is the FFT; everything else is bandwidth.

• Real FFT: pack N real samples into an N/2-point complex FFT plus a post-processing pass: about half the flops and memory of a complex FFT of the same size (realfft documents this).
• Hop: cost is inversely proportional to hop. Do not use hop N/8 unless you need smooth time sampling; hop N/2 halves the cost of hop N/4 for analysis-only features.
• Batch frames: one plan (twiddles, bit-reversal) shared for all frames; all frames independent, so parallelise over frames or run one batched FFT on GPU. Reuse one scratch buffer instead of allocating per frame.
• No extra work: skip the log and the phase when only energy is needed; pool bins (mel) immediately so you store 80 values per frame instead of 257. Whisper uses 80 mel bins from N=400, hop 160 at 16 kHz.
• Cheaper alternative: for a handful of known frequencies, a Goertzel filter costs about N multiply-adds per tone per frame, which beats the FFT when the tone count is well below log2 N.

Implementation notes

• Layout: frame-major ([t][k]) is what the FFT writes; some consumers want [k][t]: transpose once at the end.
• Frame count: 1 + (L - N) / hop with no padding; to center frame 0 at sample 0, pad N/2 samples on both ends (librosa and torch default to center=True), else the first N/2 samples are under-weighted.
• Window choice: Hann for general use. Rectangular only for transients shorter than N. Blackman-Harris (4-term, about -92 dB sidelobes) if a strong tone masks a weak one; its noise bandwidth is about 2 bins (Harris 1978).
• Scale for amplitude readout: divide by sum(w)/2 to get tone amplitude, by sqrt(sum w^2) for noise power density comparisons.
• N need not be a power of two: libraries (rustfft, pocketfft) handle 400, 1000, and so on with mixed radix, at a modest cost; prefer sizes with small prime factors if free to choose.
• f32 is enough for audio features; pick the log floor above f32 rounding noise.

Pitfalls

• Spectral leakage: a rectangular window smears strong tones over every bin at -13 dB falling slowly; use a tapered window.
• Symmetric vs periodic windows: scipy.signal.get_window default is periodic (fftbins=True), np.hanning is symmetric; mixing them breaks COLA and gives ripple in resynthesis.
• Wrong time axis: frame t is centered at t*hop + N/2, not t*hop; centered padding shifts it by N/2.
• Phase after modification: editing magnitudes and reusing the old phase leaves inconsistent frames (the STFT is redundant); the output has "phasiness". Griffin-Lim iterations reduce it.
• Zero-padding is not resolution: padding the frame to 2N interpolates bins, it does not separate close tones; only a longer window does.
• Log of zero: digital silence gives -inf without a floor.

Variants and alternatives

In production

Speech models (Whisper log-mel front end), audio codecs and fingerprinting (peaks of the spectrogram: Wang), vibration condition monitoring, radar micro-Doppler (Chen et al.), spectrum analyzers and SDR waterfalls, music-source separation.

Implementations

Sources

• J. B. Allen, L. R. Rabiner, A unified approach to short-time Fourier analysis and synthesis, Proc. IEEE 65(11), 1977, 1558-1564. doi
• D. W. Griffin, J. S. Lim, Signal estimation from modified short-time Fourier transform, IEEE Trans. ASSP 32(2), 1984, 236-243. doi
• F. J. Harris, On the use of windows for harmonic analysis with the discrete Fourier transform, Proc. IEEE 66(1), 1978, 51-83. doi
• Avery Wang, An industrial-strength audio search algorithm, ISMIR 2003. paper
• A. Radford et al., Robust Speech Recognition via Large-Scale Weak Supervision (Whisper), 2022. arXiv · code