README
¶
fft — go-fft
A pure-Go (no cgo) FFT library — the numpy.fft / scipy.fft equivalent for
Go. It computes the discrete Fourier transform of complex and real signals of
any length, with no dependency on the native FFTW3 C library.
Ruby has no cgo-free FFT (every option wraps FFTW3); gonum/dsp/fourier is pure
Go but its optimized assembly is amd64-only. This module is a fully portable
scalar core, with SIMD kernels generated across the six 64-bit Go targets
(amd64, arm64, riscv64, loong64, ppc64le, s390x) via
go-asmgen.
Status: Phase 5 — a correct pure-Go complex FFT (radix-2 Cooley–Tukey for power-of-two lengths, Bluestein's chirp-z for arbitrary lengths), the real-optimized
RFFT/IRFFT, the multi-dimensional transforms (FFT2/IFFT2,FFTN/IFFTN,RFFT2/IRFFT2), the windowing / spectral helpers (windows,FFTFreq/RFFTFreq,PSD,Spectrogram), and go-asmgen SIMD kernels (bit-identical pointwise complex multiply: SSE2 on amd64, NEON on arm64, RVV on riscv64, the vector facility on s390x — four of the six targets) behind a validated per-arch split CI, with loong64 and ppc64le on the validated scalar path (the Go assembler lacks the vector-double ops they need). The transform is also exposed to Ruby through the go-embedded-rubyFFTmodule (Phase 5). See docs/plan-fft.md for the phased roadmap.
API
import "github.com/go-fft/fft"
X := fft.FFT(x) // forward DFT of []complex128, any length
y := fft.IFFT(X) // inverse DFT, normalized by N (IFFT(FFT(x)) ≈ x)
S := fft.FFTReal(r) // forward DFT of a []float64 signal (full spectrum)
R := fft.RFFT(r) // real-input DFT, non-redundant N/2+1 bins (numpy.fft.rfft)
z := fft.IRFFT(R, len(r)) // back to []float64 (numpy.fft.irfft)
// Multi-dimensional, on flat row-major (C-order) data plus an explicit shape:
F2 := fft.FFT2(data, [2]int{rows, cols}) // 2-D DFT (numpy.fft.fft2)
d2 := fft.IFFT2(F2, [2]int{rows, cols}) // inverse 2-D DFT (numpy.fft.ifft2)
FN := fft.FFTN(data, shape) // N-D DFT (numpy.fft.fftn)
dN := fft.IFFTN(FN, shape) // inverse N-D DFT (numpy.fft.ifftn)
// Real image-style 2-D transforms (last axis keeps cols/2+1 bins per row):
RG := fft.RFFT2(img, [2]int{rows, cols}) // numpy.fft.rfft2
zG := fft.IRFFT2(RG, [2]int{rows, cols}) // numpy.fft.irfft2
// Windowing and spectral helpers:
w := fft.Hann(n) // also Hamming, Blackman, BlackmanHarris, Bartlett
f := fft.FFTFreq(n, d) // bin frequencies (numpy.fft.fftfreq)
rf := fft.RFFTFreq(n, d) // real-FFT bin frequencies (numpy.fft.rfftfreq)
p := fft.PSD(sig, d) // one-sided power spectral density (periodogram)
S := fft.Spectrogram(sig, segment, overlap, fft.Hann(segment), d) // PSD frames
// Reusable plans — precompute the twiddle tables once, amortize across calls
// (no per-call sin/cos). The convenience FFT/RFFT functions above use an
// internal per-length plan cache, so they get this for free too.
p := fft.NewPlan(n) // complex transform plan of length n
p.FFT(dst, src) // dst and src are []complex128 of length n (may alias)
p.IFFT(dst, src) // normalized inverse
rp := fft.NewRealPlan(n) // real-input transform plan
rp.RFFT(dst, src) // src []float64 (len n), dst []complex128 (len n/2+1)
rp.IRFFT(out, spec) // out []float64 (len n), spec the half spectrum
The multi-dimensional transforms are separable: the 1-D FFT is applied along
each axis in turn. The forward transforms are unnormalized; the inverses divide
by the product of the transformed axis lengths. The shape must be positive and
its product must equal len(data), else the call panics (numpy semantics).
Empty input returns an empty slice; length 1 returns a copy. RFFT keeps only
the lower N/2+1 bins because a real signal's spectrum is conjugate-symmetric
(X[N-k] = conj(X[k])); IRFFT takes the target length n explicitly.
Performance
go-fft uses a split-radix kernel for powers of two (≈⅓ fewer real
multiplies than radix-4), mixed-radix Cooley–Tukey for the other
highly-composite lengths (radix-2/3/4/5 straight-line butterflies plus a general
radix-p kernel for small primes), Rader's algorithm for primes (from N=700)
and Bluestein's chirp-z otherwise, with all twiddle factors cached per
length. It beats the pure-Go peer gonum/dsp/fourier (both CGO_ENABLED=0) at
every size — up to ~59× on primes (gonum's arbitrary-N path is O(N²)).
Against the C gold standard — native FFTW 3.3.11 plus pocketfft via
numpy.fft / scipy.fft, all single-threaded — on an Apple M4 Max (macOS
26.5), go-fft wins outright at several shapes:
| transform | go-fft | FFTW | scipy.fft | verdict |
|---|---|---|---|---|
| complex 256 | 0.79 µs | 0.42 µs | 2.3 µs | beats pocketfft ~2.9× |
| complex 1009 (prime) | 16.9 µs | 25.1 µs | 19.0 µs | beats FFTW ~1.5× |
| FFT2 1024×1024 | 2.72 ms | 9.55 ms | 5.32 ms | beats all (multicore) |
| FFT2 512×512 | 0.62 ms | 1.41 ms | 0.96 ms | beats all |
| real 1048576 | 5.54 ms | 5.14 ms | 6.37 ms | ~parity (1.08×) |
| complex 1024 | 3.44 µs | 2.06 µs | 4.47 µs | FFTW ~1.7×, beats pocketfft |
| real 1024 | 4.31 µs | 1.49 µs | 3.73 µs | FFTW ~2.9× |
| complex 10007 (prime) | 0.41 ms | 0.26 ms | 0.28 ms | FFTW ~1.6× |
go-fft wins on the large 2-D shapes (the goroutine-parallel separable path
single-threaded FFTW/pocketfft can't match), the small-N rows, and the
large-prime rows relative to FFTW's own cost. FFTW still leads the single-core
power-of-two and smooth-composite mid-range — hand-written NEON SIMD codelets a
scalar pure-Go library can't match on one core (the identified lever is
go-asmgen SIMD butterflies). Full methodology, every size, GFLOP/s, and the
per-op action items are in BENCHMARKS.md. Reproduce the
whole sweep with benchmarks/run.sh (go-fft + gonum via go test -bench, native
FFTW via a C harness, numpy/scipy via Python; correctness-gated; gonum is
isolated in the separate benchmarks/ module so the library stays
dependency-free).
Why not cgo / FFTW?
FFTW3 is a C library: binding it reintroduces a C toolchain, cross-compilation
pain, and a non-Go build. A pure-Go implementation cross-compiles to every Go
target for free and is CGO_ENABLED=0 clean — which is what makes it usable as
the FFT backend for an embedded Ruby (go-embedded-ruby) and for the wider
go-* ecosystem.
License
BSD-3-Clause. See LICENSE.
Documentation
¶
Overview ¶
Package fft is a pure-Go (cgo-free) fast Fourier transform library — the numpy.fft / scipy.fft equivalent for Go.
It computes the discrete Fourier transform (DFT) of complex and real signals of any length, with no dependency on the native FFTW3 C library. A power-of-two length uses a split-radix kernel (≈⅓ fewer real multiplies than radix-4); other lengths whose prime factors are all small use mixed-radix Cooley–Tukey (radix-2/3/4/5 straight-line butterflies plus a general radix-p butterfly for the small primes 7/11/13); a prime length uses Rader's algorithm (above a size threshold) or Bluestein's chirp-z algorithm, so any length transforms correctly and fast. Twiddle factors are precomputed and cached per length (see Plan / NewPlan), so repeated transforms of one length recompute no sin/cos.
The multi-dimensional transforms (FFT2/FFTN and their real and inverse forms) are separable: they apply a 1-D transform along each axis, and the independent lines of a large axis run in parallel across goroutines — a multicore path single-threaded references such as pocketfft cannot take. See docs/perf.md for the head-to-head benchmarks.
The forward transform follows the unnormalized convention
X[k] = sum_{n=0}^{N-1} x[n] * exp(-2πi·k·n/N)
and the inverse divides by N, so IFFT(FFT(x)) ≈ x to floating-point tolerance.
Implemented so far: FFT, IFFT and FFTReal over arbitrary lengths (Phase 0), plus the real-optimized RFFT and IRFFT (Phase 1). RFFT returns only the non-redundant N/2+1 bins of a real signal's spectrum (numpy.fft.rfft) and IRFFT inverts it back to a real signal of a caller-specified length (numpy.fft.irfft). The butterfly loops live in an internal kernels package (scalar pure-Go now) so that SIMD kernels generated by go-asmgen across the six 64-bit Go targets (amd64, arm64, riscv64, loong64, ppc64le, s390x) can drop in later without changing the public API. See docs/plan-fft.md for the full roadmap.
Index ¶
- func Bartlett(n int) []float64
- func Blackman(n int) []float64
- func BlackmanHarris(n int) []float64
- func FFT(x []complex128) []complex128
- func FFT2(data []complex128, shape [2]int) []complex128
- func FFTFreq(n int, d float64) []float64
- func FFTN(data []complex128, shape []int) []complex128
- func FFTReal(x []float64) []complex128
- func Hamming(n int) []float64
- func Hann(n int) []float64
- func IFFT(x []complex128) []complex128
- func IFFT2(data []complex128, shape [2]int) []complex128
- func IFFTN(data []complex128, shape []int) []complex128
- func IRFFT(spectrum []complex128, n int) []float64
- func IRFFT2(data []complex128, shape [2]int) []float64
- func PSD(x []float64, d float64) []float64
- func RFFT(x []float64) []complex128
- func RFFT2(data []float64, shape [2]int) []complex128
- func RFFTFreq(n int, d float64) []float64
- func Spectrogram(x []float64, segment, overlap int, window []float64, d float64) [][]float64
- type Plan
- type RealPlan
Constants ¶
This section is empty.
Variables ¶
This section is empty.
Functions ¶
func Bartlett ¶
Bartlett returns the length-n Bartlett (triangular) window:
w[k] = (2/(n-1)) · ((n-1)/2 - |k - (n-1)/2|), k = 0 .. n-1,
which rises linearly from 0 to 1 at the center and back to 0. This matches numpy.bartlett(n).
func Blackman ¶
Blackman returns the length-n Blackman window:
w[k] = 0.42 - 0.5·cos(2π·k/(n-1)) + 0.08·cos(4π·k/(n-1)).
This matches numpy.blackman(n).
func BlackmanHarris ¶
BlackmanHarris returns the length-n (4-term) Blackman–Harris window with the standard coefficients a0..a3 = 0.35875, 0.48829, 0.14128, 0.01168:
w[k] = a0 - a1·cos(2π·k/(n-1)) + a2·cos(4π·k/(n-1)) - a3·cos(6π·k/(n-1)).
This matches scipy.signal.windows.blackmanharris(n, sym=True).
func FFT ¶
func FFT(x []complex128) []complex128
FFT returns the forward discrete Fourier transform of x, using the unnormalized convention
X[k] = sum_{n=0}^{N-1} x[n] * exp(-2πi·k·n/N).
Highly-composite lengths (whose prime factors are all small) use mixed-radix Cooley–Tukey; a length with a large prime factor falls back to Bluestein's chirp-z algorithm. Twiddle factors are cached per length in an internal plan cache, so repeated transforms of the same length recompute no sin/cos. The input is not modified. Empty input returns an empty (non-nil) slice; length 1 returns a copy.
func FFT2 ¶
func FFT2(data []complex128, shape [2]int) []complex128
FFT2 is the two-dimensional FFTN: it transforms data as a row-major shape[0]×shape[1] matrix (rows then columns), matching numpy.fft.fft2.
func FFTFreq ¶
FFTFreq returns the n sample frequencies of an FFT of length n, given the sample spacing d (the inverse of the sampling rate). The result matches numpy.fft.fftfreq(n, d) exactly:
f = [0, 1, ..., (n-1)//2, -(n//2), ..., -1] / (d·n).
The DC bin is first, then the positive frequencies in increasing order, then the negative frequencies. For even n the Nyquist bin appears as -n/2; for odd n the largest positive bin is (n-1)/2. n <= 0 returns an empty (non-nil) slice.
func FFTN ¶
func FFTN(data []complex128, shape []int) []complex128
FFTN returns the forward N-dimensional discrete Fourier transform of data, interpreted as a row-major (C-order) array of the given shape. The result is a new slice of the same length; the input is not modified.
The transform is separable and unnormalized: it applies the 1-D FFT along each axis in turn, matching numpy.fft.fftn. shape must contain only positive lengths whose product equals len(data); FFTN panics otherwise (consistent with numpy raising on a mismatched shape). An empty shape transforms the single scalar unchanged.
func FFTReal ¶
func FFTReal(x []float64) []complex128
FFTReal returns the forward discrete Fourier transform of a real-valued signal as the full complex spectrum (length len(x)). It is a convenience wrapper that promotes x to complex128 and calls FFT.
func Hamming ¶
Hamming returns the length-n Hamming window:
w[k] = 0.54 - 0.46·cos(2π·k/(n-1)), k = 0 .. n-1.
This matches numpy.hamming(n).
func Hann ¶
Hann returns the length-n Hann (Hanning) window:
w[k] = 0.5 - 0.5·cos(2π·k/(n-1)), k = 0 .. n-1.
This matches numpy.hanning(n).
func IFFT ¶
func IFFT(x []complex128) []complex128
IFFT returns the inverse discrete Fourier transform of x, normalized by N so that IFFT(FFT(x)) ≈ x to floating-point tolerance:
x[n] = (1/N) * sum_{k=0}^{N-1} X[k] * exp(+2πi·k·n/N).
The input is not modified.
func IFFT2 ¶
func IFFT2(data []complex128, shape [2]int) []complex128
IFFT2 is the two-dimensional IFFTN, the inverse of FFT2, normalized by shape[0]*shape[1]. It matches numpy.fft.ifft2.
func IFFTN ¶
func IFFTN(data []complex128, shape []int) []complex128
IFFTN returns the inverse N-dimensional discrete Fourier transform of data, the inverse of FFTN. It applies the 1-D IFFT along each axis, so the result is normalized by the product of the transformed axis lengths and IFFTN(FFTN(x), shape) ≈ x. The input is not modified; the same shape rules as FFTN apply.
func IRFFT ¶
func IRFFT(spectrum []complex128, n int) []float64
IRFFT inverts RFFT, reconstructing a real-valued signal of length n from the n/2+1 non-redundant spectral bins produced by RFFT. The caller passes the desired output length n explicitly, because n and n-1 yield the same number of kept bins so the half spectrum alone is ambiguous; this mirrors numpy.fft.irfft(spectrum, n).
spectrum is read up to min(len(spectrum), n/2+1) bins; any bins beyond that (or beyond what the Hermitian mirror needs) are treated as zero. The result is normalized by n so that IRFFT(RFFT(x), len(x)) ≈ x.
The input is not modified. n <= 0 returns an empty (non-nil) slice.
func IRFFT2 ¶
func IRFFT2(data []complex128, shape [2]int) []float64
IRFFT2 inverts RFFT2, reconstructing a real row-major matrix of shape shape[0]×shape[1] from a spectrum laid out as shape[0]×(shape[1]/2+1) complex bins (the layout RFFT2 produces). The complex inverse is applied down the columns first, then the real inverse (irfft) along each row, with the target row length shape[1] supplied explicitly (since shape[1] and shape[1]-1 share a bin count). The result is normalized by shape[0]*shape[1] so that IRFFT2(RFFT2(x, shape), shape) ≈ x. This matches numpy.fft.irfft2.
shape lengths must be positive; IRFFT2 panics otherwise. data is read up to shape[0]*(shape[1]/2+1) bins; any beyond that are treated as zero. The input is not modified.
func PSD ¶
PSD returns the one-sided power spectral density (periodogram) of the real-valued signal x sampled at spacing d, using the non-redundant RFFT bins. The result has length len(x)/2+1.
Each bin is scaled to a density: |X[k]|² / (fs·N) where fs = 1/d is the sampling rate and N = len(x). To conserve total power, every bin except DC (and the Nyquist bin for even N) is doubled to fold in the discarded mirror half. This matches scipy.signal.periodogram(x, fs, window="boxcar", scaling="density", return_onesided=True). The companion frequency axis is RFFTFreq(len(x), d). The input is not modified; an empty x returns an empty (non-nil) slice.
func RFFT ¶
func RFFT(x []float64) []complex128
RFFT returns the forward discrete Fourier transform of a real-valued signal, keeping only the non-redundant spectrum. For input of length N the result has length N/2+1 (integer division): bins 0..N/2. The discarded upper half is the conjugate mirror of the lower half (X[N-k] = conj(X[k])), so it carries no new information. This matches numpy.fft.rfft.
The transform is unnormalized, consistent with FFT:
X[k] = sum_{n=0}^{N-1} x[n] * exp(-2πi·k·n/N), k = 0 .. N/2.
For even N the real signal is packed into an N/2-point complex transform and untangled — about twice as fast as a full complex FFT — using twiddle tables cached per length. The input is not modified. Empty input returns an empty (non-nil) slice.
func RFFT2 ¶
func RFFT2(data []float64, shape [2]int) []complex128
RFFT2 returns the forward 2-D DFT of a real row-major matrix of the given shape (shape[0] rows × shape[1] columns). The real transform is applied along the last axis, so each output row keeps shape[1]/2+1 non-redundant bins; the result is a row-major matrix of shape shape[0]×(shape[1]/2+1). The full complex transform is then applied down the columns. This matches numpy.fft.rfft2.
shape lengths must be positive and shape[0]*shape[1] must equal len(data); RFFT2 panics otherwise. The input is not modified.
func RFFTFreq ¶
RFFTFreq returns the n/2+1 sample frequencies of a real FFT (RFFT) of input length n, given the sample spacing d. The result matches numpy.fft.rfftfreq(n, d) exactly:
f = [0, 1, ..., n//2] / (d·n), length n//2+1.
All frequencies are non-negative because RFFT keeps only the non-redundant lower half of the spectrum. n <= 0 returns an empty (non-nil) slice.
func Spectrogram ¶
Spectrogram computes a sequence of one-sided power spectra over successive, optionally overlapping segments of the real signal x, the building block of a time–frequency display. segment is the segment length (window size) and overlap is the number of samples shared between consecutive segments (0 <= overlap < segment). Each segment is multiplied by window (which must have length segment) before its PSD is taken.
The return value is a slice of frames, each a []float64 of length segment/2+1 (the RFFTFreq(segment, d) axis); frame t covers samples [t·step, t·step+segment) with step = segment-overlap. Segments are taken only while a full window fits, so trailing samples shorter than segment are dropped (the common STFT convention). The input is not modified.
Spectrogram panics if segment <= 0, len(window) != segment, or overlap is out of range, mirroring the contract of scipy.signal.spectrogram.
Types ¶
type Plan ¶
type Plan struct {
// contains filtered or unexported fields
}
A Plan is a reusable, precomputed transform of a fixed length N. Building a plan factors N and precomputes every twiddle factor once; the resulting tables are then amortized across an unbounded number of FFT/IFFT calls, so repeated transforms of the same length pay no sin/cos cost. Plans are immutable after construction and safe for concurrent use by multiple goroutines (the transform methods write only into the caller-supplied destination and per-call scratch).
For one-off transforms the package-level FFT/IFFT functions are more convenient; they consult an internal plan cache so even they avoid recomputing twiddles for a length they have seen before.
func NewPlan ¶
NewPlan returns a transform plan for length n, precomputing all twiddle factors. n may be any non-negative integer; n == 0 and n == 1 produce a trivial plan. The same plan serves both the forward FFT and the inverse IFFT.
func (*Plan) FFT ¶
func (p *Plan) FFT(dst, src []complex128) []complex128
FFT writes the forward DFT of src into dst and returns dst. dst and src must each have length Len(); dst may alias src. src is not modified unless it aliases dst.
func (*Plan) IFFT ¶
func (p *Plan) IFFT(dst, src []complex128) []complex128
IFFT writes the inverse DFT of src into dst (normalized by N) and returns dst. dst and src must each have length Len(); dst may alias src.
type RealPlan ¶
type RealPlan struct {
// contains filtered or unexported fields
}
A RealPlan is a reusable, precomputed real-input transform of a fixed length N. Like Plan it amortizes all twiddle computation across calls; it is the real-signal counterpart used by RFFT/IRFFT. RealPlans are immutable after construction and safe for concurrent use.
For even N the plan packs the N real samples into an N/2-point complex transform and untangles the result with a precomputed table — about twice as fast as a full complex FFT. For odd N (which cannot be packed) it wraps a full complex plan of length N.
func NewRealPlan ¶
NewRealPlan returns a real-input transform plan for length n, precomputing all twiddle factors. n must be non-negative.
func (*RealPlan) IRFFT ¶
func (p *RealPlan) IRFFT(dst []float64, src []complex128) []float64
IRFFT writes the real signal of length Len() reconstructed from the half spectrum src into dst and returns dst. dst must have length Len(); src is read up to min(len(src), Len()/2+1) bins. The result is normalized by N. src is not modified.
For even N the inverse mirrors the forward packing: the half spectrum is pre-processed into an N/2-point complex spectrum, run through one N/2-point inverse complex FFT, and unpacked — half the arithmetic and memory traffic of promoting to a full conjugate-symmetric length-N inverse transform. For odd N (and the trivial N<=1) it falls back to the full conjugate-mirror inverse.
func (*RealPlan) RFFT ¶
func (p *RealPlan) RFFT(dst []complex128, src []float64) []complex128
RFFT writes the non-redundant N/2+1 spectral bins of the real signal src into dst and returns dst. src must have length Len(); dst must have length Len()/2+1. src is not modified.
Source Files
¶
Directories
¶
| Path | Synopsis |
|---|---|
|
internal
|
|
|
kernels
Package kernels holds the inner loops of the FFT — the bit-reversal permutation and the radix-2 butterfly stages — isolated from the public API.
|
Package kernels holds the inner loops of the FFT — the bit-reversal permutation and the radix-2 butterfly stages — isolated from the public API. |