fft

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-ruby FFT module (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
  • func NewPlan(n int) *Plan
  • func (p *Plan) FFT(dst, src []complex128) []complex128
  • func (p *Plan) IFFT(dst, src []complex128) []complex128
  • func (p *Plan) Len() int
• type RealPlan
  • func NewRealPlan(n int) *RealPlan
  • func (p *RealPlan) IRFFT(dst []float64, src []complex128) []float64
  • func (p *RealPlan) Len() int
  • func (p *RealPlan) RFFT(dst []complex128, src []float64) []complex128

Functions

func Bartlett

func Bartlett(n int) []float64

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

func Blackman(n int) []float64

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

func BlackmanHarris(n int) []float64

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

func FFTFreq(n int, d float64) []float64

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

func Hamming(n int) []float64

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

func Hann(n int) []float64

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

func PSD(x []float64, d float64) []float64

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

func RFFTFreq(n int, d float64) []float64

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

func Spectrogram(x []float64, segment, overlap int, window []float64, d float64) [][]float64

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

func NewPlan(n int) *Plan

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.

func (*Plan) Len

func (p *Plan) Len() int

Len reports the transform length the plan was built for.

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

func NewRealPlan(n int) *RealPlan

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) Len

func (p *RealPlan) Len() int

Len reports the real-input length the plan was built for.

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.