BFFTThe Bruun transform,
repaired.

In 1978, G. Bruun described a real-coefficient route through the discrete Fourier transform. It was remembered as clever, awkward, and numerically suspect. BFFT asks a different question: what if the factor tree was sound and the local coordinates were standing wrong?

Read On Bruun, RevisitedPDF · 12 pages ↗ Inspect the implementation Source, tests, and history ↗
Historical origin
Bruun, 1978
Public implementation
C / C++17 / Python
Forward routes
DIF · DIT · DIP
Mathematical result
Exact · condition one
License
MIT

01 / THE OBJECT

A real FFT whose intermediate states remain meaningful.

BFFT computes the ordinary discrete Fourier transform of a real signal. Its standard output has the familiar bins from DC through Nyquist, and its inverse returns the original samples. The unusual part lies beneath that interface: a normalized form of Bruun’s real factorization, several lawful orders for walking it, and direct access to the native residue coordinates.

That distinction matters. Most users need only a reliable transform. A researcher may also want to pause inside it, retain its two-dimensional intermediate lattice, filter in residue space, follow the difference between decimation-in-frequency and decimation-in-time, or build a transform whose product includes support and timing information. BFFT makes those possibilities branches of one object instead of unrelated post-processing tricks.

In one sentence

BFFT takes Bruun’s recursive real-polynomial factorization, expresses every quadratic leaf in a normalized complex frame, and turns that repaired geometry into a usable forward, inverse, and intermediate-state library.

The paper closes the algebraic loop.

The implementation showed that the repaired transform worked. On Bruun, Revisited now proves why: every reduction edge has an explicit inverse, every normalized cell has a uniform metric, and terminal evaluation is exactly the ordinary real DFT after one declared permutation.

Two-sided reconstructionGNFN = I

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

Unit quadratic charteθ = (z − cosθ) / sinθ
eθ2 = −1

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

Uniform cell metricTθTTθ = 2I

Every nontrivial cell is a rotation followed by a Hadamard merge. The angle changes; the energy law does not.

Exact Fourier boundaryPNBN = DNR

Reduction commutes with evaluation at the roots. Therefore BFFT is a factorization of the DFT itself, not a numerically similar transform.

Global consequenceBNTBN = NI

After the standard packing weights, every singular value is √N and the transform condition number is exactly one.

On Bruun, RevisitedRead the proofs, examples, audits, and appendices · PDF ↗

02 / HISTORICAL LINE

From filter bank to forgotten branch to working library.

  1. Cooley and Tukey popularize the recursive FFT

    Their machine algorithm made the factor-and-reuse strategy standard. It was not the only way to factor a DFT, but it became the dominant 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 but 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 state-of-the-art account. By then the question was no longer whether many FFTs existed, but 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 changes the coordinate system

    The public BFFT repository began on June 10. The implementation replaced the raw quadratic residue basis with normalized unit-circle frames, completed a matched inverse, added DIT and diagonal walks, and built C, C++, Python, Numba, STFT, BODFT, magnitude, phase, native-layout, and residue-filter interfaces around the same core.

    BFFT repository and history

03 / BRUUN’S FACTOR TREE

The alternative is in the factorization, not the answer.

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

Why it was attractive

For a real input, conjugate frequency bins contain redundant information. Bruun’s factor tree works with real coefficients until the leaves, aligning the arithmetic with that symmetry rather than running a full complex transform and throwing half of it away.

Why it looked troublesome

The familiar raw recurrence uses coefficients such as 2 cos θ inside quadratic residue coordinates. Near particular angles, the basis can become badly conditioned and inverse formulas inherit subtractive cancellation. A numerical defect in those coordinates was easily read as a defect in the factorization itself.

BFFT’s wager

Do not tune the unstable coordinates harder. Replace each local quadratic pair with the same normalized complex plane, so every twiddle is a point on the unit circle and every merge has the same energy law.

04 / THE COORDINATE REPAIR

Every leaf becomes the same little complex plane.

The normalized butterfly rotates the odd child by (cos θ, sin θ), then combines it with the even child. Rotation preserves the odd pair’s energy; the plus/minus merge doubles total energy uniformly. No direction receives a special scale.

The raw monomial chart (1, z) has κ₂ = cot(θ/2) for 0° < θ ≤ 90°. The normalized chart remains 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ᵦ)

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

The same signed map admits conservative physical representatives.

The paper does not claim that silicon has become unnecessary. It proves a narrower and more useful statement: any signed BFFT stage can be lifted to nonnegative rails, gauged into a mass-preserving transport or a convex equilibrium map, and projected back to the exact Fourier observable.

Encodex ↦ (x⁺, x⁻)

A signed coordinate becomes two nonnegative rails whose difference is x.

TransportC(M) ≥ 0

Positive and negative matrix parts route mass without changing the signed projection.

ObserveΠC(M) = MΠ

Subtract the rails at the boundary and recover the original stage exactly.

Conservative transport

A backward diagonal gauge makes each lifted stage column-stochastic. Total two-rail mass is conserved; equal positive and negative common mode can be removed without altering any Fourier coordinate.

Convex equilibrium

A forward gauge makes the stage row-stochastic. Each output becomes a convex combination, with a unique clamped quadratic-energy equilibrium at the desired transformed state.

Passive phase packets

Before real folding, every normalized DIP butterfly is a unit phase delay followed by a lossless balanced coupler. The folded form adds only the declared real reflection.

05 / THREE WALKS THROUGH ONE GEOMETRY

DIF, DIT, and DIP differ by transport order.

DIF

Factor-tree descent

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

Best understood as: finish one region, then descend.
DIT

Spectral ascent

Small seed blocks enter in bit-reversed block order and merge upward. Angles vary by position rather than node. Standard spectrum order falls out of the walk without a final permutation, exchanging some full-array stage sweeps for a simpler exit.

Best understood as: build local spectra, then merge.
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. It can be paused, inspected, aligned, or resumed to the same final Fourier product.

Best understood as: keep time and frequency alive together.

06 / SPECIAL BREAKOUT

The Diagonal Intermediate Phase transform.

DIP is not another display made after an FFT. It is an ordering of the transform itself. 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. Columns j carry the low time bits. Coarse frequency survives as phase along a row; coarse time survives as phase down a column.

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

A continuum of useful interiors

Stage 0 is the time signal. Stage 8 is the final spectrum. The interior stages are exact transform states, not approximations between the two.

Time-frequency motion has a native place

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

Pause and finish are ordinary operations

A recovered or modified intermediate state can continue through the remaining feed-forward cells. No inverse-to-time and second FFT are required merely to re-enter the transform.

The stage stack has fixed geometry

Every stage is a scaled unitary map and stores exactly N complex cells. Adjoints, momentum, and cross-stage comparison therefore operate on known conditioning.

Support can travel inside the transform

The Zak boundary law turns a leading edge into complete comb rows plus one partial row. That creates an intrinsic route for timing, support, and certified phase-disk selection rather than an endpoint heuristic added afterward.

Fractional channels are members of one family

Half-bin BODFT is the α = ½ case of the same twisted cell family. Finer dyadic offsets can share transport, opening a path toward support-aware and correlated transforms.

07 / THE PUBLIC LIBRARY

A repaired algorithm still has to be ordinary software.

CorePower-of-two real transforms, N ≥ 4; double and float32; forward and inverse
LayoutsStandard FFT order, native Bruun order, and real residue coordinates
InterfacesStable C ABI, C++17 RAII wrapper, Python/NumPy package, Numba-callable plans
StreamingReusable workspaces, STFT plans, magnitude-only and magnitude/phase paths
Related transformsExact odd-frequency DFT at half bins, BODFT interface, and forward-only Fast Correlated Transform
FilteringConvert standard responses to residue filters, apply them, and invert directly
BuildMake, CMake, pkg-config, package config, Linux, macOS, and Windows-oriented builds
VerificationRound-trip, direct-DFT, NumPy, layout, SIMD, standards-mode, sanitizer, and benchmark probes
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)

08 / RECEIPTS AND BOUNDARIES

What is established, and where to inspect it.

Established in the implementation

  • The normalized forward and matched inverse return ordinary Fourier bins and round-trip at floating-point precision.
  • DIF and DIT share the same unit-frame leaves and differ in ordering and transport pressure.
  • DIP stage states satisfy the finite Zak identity numerically through the checked sizes.
  • Public APIs, builds, tests, examples, and source are available under the MIT license.

Active research beyond the base library

  • Intrinsic phase-disk support selection is exact in its reference and joined-kernel experiments, but its strongest shared-transport complexity target remains research.
  • Mid-level fusion and super-resolution results are demonstrations of the DIP geometry, not a universal replacement for ordinary spectral methods.
  • Physical delay-line and optical interpretations are designs over the same walk, not shipped hardware.
1978G. Bruun, z-transform DFT filters and FFT’s 1987Sorensen et al., real-valued FFT algorithms 1990Duhamel & Vetterli, tutorial review 2005Frigo & Johnson, FFTW3 2008MIT polynomial-factor account 2026On Bruun, Revisited — complete proofs and audits (PDF) CodeBFFT source, history, tests, and documentation