BFFT

G. Bruun's real-coefficient Fourier transform of 1978, exact and perfectly conditioned.

Read On Bruun, Revisited (PDF, 12 pages)Source on GitHub

The transform

In 1978 G. Bruun factored the real Fourier transform using real coefficients alone, and the method was set aside as numerically unreliable. The trouble was the coordinate each stage worked in. Rewritten in a normalized frame, every stage is a rotation followed by a sum-and-difference merge, the whole transform has condition number one, and it equals the DFT up to a stated reordering.

BFFT is a MIT-licensed, SIMD-optimized C and C++ library built on that proof, with Python bindings. It computes the ordinary discrete Fourier transform of a real signal: the standard output holds the familiar bins from DC through Nyquist, and the inverse returns the original samples. Beneath that interface lie the normalized form of Bruun's factorization, three orders for walking it and direct access to its native residue coordinates, so a caller can stop inside the transform, filter there, or resume it.

import numpy as np
import bfft

x = np.random.randn(1024)
X = bfft.rfft(x)       # ordinary bins: 0 … N/2
y = bfft.irfft(X)      # returns the real samples

assert np.allclose(y, x)

Bruun's factor tree

x(z) = Σ x[n]zⁿ reduce through the real factors of zᴺ − 1 z²ᴹ + a zᴹ + 1, |a| ≤ 2 finish at conjugate roots of unity X[k] = x(e−2πik/N)

A DFT evaluates the polynomial x(z) at the N-th roots of unity, so it can be computed by reducing x(z) modulo a tree of factors of zN − 1. For a real input the conjugate bins carry the same information, and Bruun's tree keeps every coefficient real until the leaves, aligning the arithmetic with that symmetry instead of running a complex transform and discarding half of it.

The usual recurrence keeps each quadratic leaf in the raw monomial coordinates (1, z), with coefficients such as 2 cos θ. Near particular angles that basis is badly conditioned and the inverse formulas suffer subtractive cancellation, which is the unreliability for which the method was set aside. BFFT keeps the tree and writes every leaf in the same normalized complex plane, so each twiddle is a point on the unit circle and every merge obeys one energy law.

The normalized frame

The normalized butterfly rotates the odd child by (cos θ, sin θ), then adds it to and subtracts it from the even child. The rotation preserves the odd pair's energy, and the sum-and-difference merge doubles the total energy uniformly, so no direction receives a special scale. Move the angle and the odd child below; the stage ratio stays at exactly two.

The raw monomial chart (1, z) has κ₂ = cot(θ/2) for 0° < θ ≤ 90°. The normalized chart stays at κ₂ = 1.

Input energy
—
Output energy
—
Stage ratio
2.000000×
Raw chart κ₂
—
Normalized chart κ₂
1.000000

(eₐ, eᵦ, oₐ, oᵦ) → (eₐ + r, eᵦ + i, eₐ − r, i − eᵦ)

Here r = cos θ·oₐ − sin θ·oᵦ and i = sin θ·oₐ + cos θ·oᵦ. The inverse applies the adjoint cell with a factor of ½ at each level, and across log₂N levels those factors become the ordinary 1/N inverse normalization.

The proof

On Bruun, Revisited proves four linked results: every reduction edge has an explicit inverse, every conjugate-pair quotient has a unit complex chart, every normalized cell has the same metric, and evaluation at the leaves is exactly the real DFT after one stated permutation.

Two-sided reconstruction

GNFN = I

A closed-form Chinese-remainder merge reverses every sibling split, so reduction and reconstruction are mutual inverses, locally and globally.

Unit quadratic chart

eθ = (z − cos θ) / sin θ, eθ2 = −1

Each real conjugate-pair quotient becomes an ordinary Euclidean complex line, with no complex polynomial arithmetic.

Uniform cell metric

TθTTθ = 2I

Every nontrivial cell is a rotation followed by a sum-and-difference merge. Only the angle changes from cell to cell, so every cell obeys the same energy law.

Exact Fourier boundary

PNBN = DNR

Reduction commutes with evaluation at the roots, so BFFT is a factorization of the DFT itself.

BNTBN = NI

After the standard packing weights, every singular value is √N and the condition number of the whole transform is exactly one. In binary64 the largest error against a direct DFT is 6.66 × 10⁻¹⁶ at N = 16 and 7.11 × 10⁻¹⁴ at N = 4096, and the measured condition number is at most 1.000000000000001.

History

  1. Cooley and Tukey popularize the recursive FFT

    Their machine algorithm made the factor-and-reuse strategy standard. It was one of several ways to factor a DFT, and it became the vocabulary against which later algorithms were described.

    Original paper

  2. G. Bruun starts from z-transform DFT filters

    Bruun showed that a DFT could be implemented as a filter bank with fewer coefficients. The real-signal form kept real coefficients through the factor tree and, in Bruun's accounting, used half as many real multiplications as the classical FFT.

    Bruun, “z-transform DFT filters and FFT's”

  3. Real-data FFTs receive a systematic comparison

    Sorensen, Jones, Heideman and Burrus described construction methods across the major FFT families and presented a real split-radix implementation with a lower operation count. Bruun's branch remained mathematically interesting and did not become the practical default.

    Real-Valued Fast Fourier Transform Algorithms

  4. The algorithm families are surveyed as a mature field

    Duhamel and Vetterli's tutorial review placed Cooley–Tukey, split-radix, prime-factor, Winograd and polynomial-factor approaches in one account. The open question had become which arithmetic and implementation structure deserved to survive.

    Tutorial review

  5. Planning and memory hierarchy become part of the algorithm

    FFTW demonstrated that operation count alone does not determine speed on modern processors. Codelets, runtime planning, SIMD and memory movement became first-class design objects.

    The Design and Implementation of FFTW3

  6. Bruun remains an instructive alternate factorization

    Steven G. Johnson's MIT FFT notes presented the Bruun construction directly, as a recursive factorization of zN − 1 into real quadratic factors whose constants remain real until the final step.

    MIT FFT notes

  7. BFFT writes every leaf in a normalized complex frame

    The raw quadratic residue basis gives way to unit-circle frames, with a matched inverse, decimation-in-time and diagonal walks, and C, C++ and Python interfaces on the same core. On Bruun, Revisited proves the result.

    BFFT repository

A physically realizable transform

Every stage of BFFT can be carried out with quantities that are never negative, as mass moved between two rails or light divided by lossless couplers, and the signal is simply the difference between the rails. That makes the transform something a physical system could perform, not only a computer.

Every signed BFFT stage can be lifted onto two nonnegative rails, x ↦ (x⁺, x⁻), whose difference is the signal. A backward diagonal gauge makes each lifted stage column-stochastic, a conservative transport that keeps the total mass on the rails, and a forward gauge makes it row-stochastic, so that each output is a convex combination with a unique equilibrium at the transformed state. Subtracting the rails at the boundary recovers the stage exactly, ΠC(M) = MΠ, and in the paper's N = 512 audit the stochasticity defect is 2.22 × 10⁻¹⁶. Before real folding, every normalized DIP butterfly is a unit phase delay followed by a lossless balanced coupler. Each stage of the Fourier transform is therefore exactly realizable with ideal passive elements: nonnegative transports, convex junctions and lossless phase couplers.

Three walks

One normalized geometry can be traversed in three orders. They share the same unit-frame leaves and differ in the order of transport.

DIF, factor-tree descent

Natural-order samples enter. One node angle applies across a whole block, so the interior is an isoclinic rotation that widens cleanly across SIMD lanes. Working sets shrink as the walk descends; the spectrum emerges in Bruun's residue order and can be packed into standard FFT order.

DIT, spectral ascent

Small seed blocks enter in bit-reversed block order and merge upward, with angles that vary by position. The standard spectrum order falls out of the walk without a final permutation, in exchange for some full-array stage sweeps.

DIP, diagonal phase packets

The walk resolves one time bit and one frequency bit together. At every level its N scalars form a two-dimensional lattice with a fine-frequency row axis and a fine-time column axis, and it can be paused, inspected, aligned or resumed to the same final Fourier product.

The DIP lattice

The diagonal intermediate phase walk is an ordering of the transform itself, and each of its stages is an exact transform state. After stage t, let e = 2t and q = N/e. The state is an e × q finite Zak lattice,

Bt[δ, j] = Σr x[j + rq] e−2πiδr/e

Rows δ carry the low frequency bits and columns j carry the low time bits. Coarse frequency survives as phase along a row, and coarse time as phase down a column. Stage 0 is the time signal and stage 8 is the final spectrum of these 256 samples.

Rows, frequency refinement
16
Columns, time refinement
16
Stored complex cells
256

Time-frequency motion

A two-axis lattice translation can align diagonal motion, such as a chirp or a shifted observation, which no frequency-only shift of the final spectrum can express.

Pause and finish

A recovered or modified intermediate state continues through the remaining feed-forward cells. Re-entering the transform needs no inverse to time and no second FFT.

Fixed geometry

Every stage is a scaled unitary map and stores exactly N complex cells, so adjoints, momentum and comparison across stages all operate at known conditioning.

Fractional channels

The half-bin BODFT is the α = ½ case of the same twisted cell family.

The library

Core
Power-of-two real transforms, N ≥ 4, in double and single precision, forward and inverse, with SIMD paths built for the host processor
Layouts
Standard FFT order, native Bruun order and real residue coordinates
Interfaces
A stable C ABI, a C++17 RAII wrapper, a Python and NumPy package, and Numba-callable plans
Streaming
Reusable workspaces, STFT plans, and magnitude-only and magnitude-and-phase paths
Related transforms
The exact odd-frequency DFT at half bins (BODFT) and the forward-only Fast Correlated Transform
Filtering
Standard frequency responses converted to residue filters, applied and inverted directly
Build
Make, CMake and pkg-config, on Linux, macOS and Windows
Tests
Round-trip, direct-DFT, NumPy, layout, SIMD, standards-mode and sanitizer checks

Experiments

The repository is also where other work built on the transform lives. BFFT Cartoon splits an image into a smooth cartoon layer and signed texture with a Meyer–Bregman decomposition solved on BFFT's spectral path, fast enough to run on live video as a filter for OBS. FlowCells takes its image geometry from one such decomposition. Inside the library, the same normalized cells give the exact half-bin transform, BODFT, and a forward-only Fast Correlated Transform that reports each bin at its most correlated leading-edge slice.

Sources