Documentation

Signal processing

FFT conventions, zero-padding, and how circular wrap is avoided.

Hard-won rules for DSP work in this codebase. Read before adding or changing any FFT-based op. Most of these come from incidents where the natural-looking implementation produced non-physical results.

Rule 1: Always zero-pad FFT-based convolution / filtering

The FFT is circular. A multiplication in the frequency domain is equivalent to a circular convolution in the time domain, not a linear one. If you mean linear convolution (the physical thing: a real LTI filter applied to a finite signal), you have to zero-pad before the FFT so the wrap-around can't overlap your signal.

The rule: pad the input to length

L_fft >= N_signal + N_impulse - 1

before the FFT, perform the spectral op, IFFT, then take the first N_signal (or appropriate window) of the result.

This applies whether the spectral op is: - multiplication by a filter spectrum (FFT-domain filtering) - multiplication by another signal's conjugated spectrum (cross-correlation) - multiplication by a phase ramp (fractional time-delay) - truncation of high frequencies (decimation, LPF, brick-wall)

If you omit the pad, you get circular wrap. Echo at the trailing edge of the buffer reappears at the leading edge. Energy from the last few range bins shows up in the first few. None of that is physical: it is purely an artifact of the FFT being defined on a periodic domain.

When N_impulse is "infinite" (brick-wall LPFs)

A brick-wall LPF: spectral rectangular window: has an infinite- support sinc kernel in time. There is no N_impulse that captures it exactly. The kernel amplitude at offset n decays as roughly 1/(π·n). So:

• Pad by 64 samples → wrap residual ~ 1/(π·64) ≈ −46 dB
• Pad by 256 samples → ~ −58 dB
• Pad by 1024 samples → ~ −70 dB
• Pad by 4096 samples → ~ −82 dB

Default in simulator/gpu_echo._fft_decimate: pad = max(K·256, 1024) upsample samples (~ −70 dB wrap residual). Sufficient for SAS-image DR (~40 dB). Bump if you need tighter parity against a reference.

If the brick-wall LPF residual is unacceptable, switch to a finite- support windowed-sinc filter (Hann/Hamming/Kaiser): that gives an honest N_impulse and a clean pad bound. We use windowed-sinc filters in beamformer/sas/preprocessing._fir_filtfilt_gpu (proper pad already).

Rule 2: Linear convolution length for FIR via FFT

For convolving an N_signal-length signal with an N_taps-length FIR filter, the FFT length must be ≥ N_signal + N_taps − 1. Round up to the next power of 2 (cuFFT loves smooth lengths).

Reference implementation: beamformer/sas/preprocessing.py ; _fir_filtfilt_gpu and mf_compress_all_pings_gpu both use fft_size = next_pow2(padded_len + n_taps - 1). Use them as a template.

Rule 3: Cross-correlation via FFT needs 2N padding

corr = ifft(fft(a) * conj(fft(b))) is the circular cross-correlation of a and b. The lag-L output gets contributions from lag L−N (circular wrap). To get clean linear xcorr you have to zero-pad both inputs to length ≥ 2N − 1 before the FFT.

If you Tukey-window the inputs first (as the delay-estimation code does), the wrap residual is suppressed: but only if you trust the window. Prefer explicit zero-padding when you can.

Rule 4: Fractional time-delay via spectral phase ramp is circular

y = ifft(fft(x) * exp(-1j·omega·tau)) shifts x by tau seconds, but the shift wraps around the buffer. Energy that "falls off" the trailing edge reappears at the leading edge.

To get a non-wrapping fractional delay, zero-pad x by at least ceil(|tau_max|·fs) samples (plus a guard for the sinc tail of the fractional-delay LPF) before the FFT.

Rule 5: Sampled signals can alias the source

A band-limited continuous-time signal with main support |f| < fs/2 may still have spectral wings beyond fs/2 (the support is "almost flat-top" but never exact). Sampling at fs aliases those wings back into the band. The aliasing magnitude depends on the source's spectral roll-off.

Example: the LFM chirp exp(j·π·BR·t²) over a rect-windowed support [−T/2, T/2] has spectrum approximately rect of width BW plus sinc-decaying wings. With BW = 30 kHz and fs = 40 kHz, the chirp's wings extend slightly past ±20 kHz, and direct native-fs sampling aliases ~−42 dB of energy back into the band.

Recipe: if you need a clean band-limited sampling at native fs, sample at a higher rate first (e.g., fs_up = K·fs), apply a proper- padded brick-wall LPF (Rule 1) to ±fs/2, then decimate by K. Do not trust direct analytic evaluation at native fs.

Reference implementation: simulator/gpu_echo._polyphase_chirp_at_fs builds each polyphase chirp row this way.

Rule 6: Operations that legitimately want circular FFTs

Not every FFT-based op needs zero-padding. The following are clean because circular topology is what's actually wanted:

• Terrain / heightfield synthesis by FFT shaping of a random spectrum (simulator/seafloor_builder.py, simulator/procedural_features.py, simulator/ripples.py): the terrain is meant to tile, so circular wrap = correct seamless boundary.
• Sinc interpolation by spectrum zero-padding: beamformer/sas/slc.py range/along-track interpolation. This is pure resampling: no impulse-response convolution is added, so there's nothing to wrap.
• Spectrograms / Welch periodograms: each window is its own block; wrap-within-block is by design.

When in doubt, ask: "if I made my FFT length 100× larger, would the result of the op change?" If yes, you have wrap and you need to pad. If no, you're computing something genuinely periodic.

Audit references

A full repo audit at 2026-05-24 cataloged every FFT-based op in the codebase and marked which obey Rule 1. See the audit dispatched in the 2026-05-24 polyphase echo work. Known offenders still pending fix:

• simulator/fw_engine.py:1520-1522: per-channel time-delay phase ramp without zero-pad (Rule 4). Low priority: FW backend is validation-only.
• beamformer/sas/delay_estimation.py:413, 681: FFT cross-correlation without 2N pad (Rule 3). Tukey-windowed inputs suppress wrap in practice; flag for cleanup if delay-peak quality regresses.

Default operator checklist (apply when adding new FFT ops)

1. State N_signal and the impulse-response length / phase-ramp magnitude.
2. Compute L_fft >= N_signal + N_impulse - 1 (with sinc tail margin for brick-walls).
3. Round L_fft up to a cuFFT-friendly length (next_pow2 or _next_cufft_smooth).
4. Zero-pad the input to L_fft.
5. FFT, spectral op, IFFT, divide by any scaling the op implies.
6. Discard the trailing pad: return the first N_signal samples.
7. Write a unit test that verifies parity with a reference FIR convolution at the buffer extremes (where wrap would show first).

If any of these steps don't apply to your op, write down why in the docstring next to the implementation.