Derive a radix-two FFT and count its work

Try

A Fourier transform produces every frequency bin from every input coordinate. Does that require repeating all the arithmetic separately for each bin? At length eight, start with the real input x=(1,2,0,0,0,0,0,0). Its defining transform is particularly transparent: Xₖ=1+2ω₈^k. We will recover these same eight values by sharing smaller transforms, then prove why that sharing works for arbitrary complex inputs of power-of-two length.

Our convention is fixed throughout: indices begin at zero, ωₙ=exp(−2πi/n), and Xₖ=Σ from j=0 to n−1 of xⱼωₙ^(jk). There is no forward normalization. The inverse has a positive exponent and a factor 1/n. A twiddle factor is simply one of the complex numbers ωₙ^k used in a butterfly; the name does not introduce a different mathematical operation.

Separate the input into its even positions e=(x₀,x₂,x₄,x₆)=(1,0,0,0) and odd positions o=(x₁,x₃,x₅,x₇)=(2,0,0,0). Their four-point transforms are E=(1,1,1,1) and O=(2,2,2,2). Predict how their frequency bins might combine without computing either smaller transform again.

Predict

For a length n=2m input, split the defining sum into even and odd j. The even contribution uses powers ωₙ^(2j k), equal to ωₘ^(jk) because ωₙ²=ωₘ. The odd contribution has one additional factor ωₙ^k. Thus for 0≤k<m,

Xₖ=Eₖ+ωₙ^k Oₖ.

For bin k+m, the child transforms repeat: E and O have frequency period m. Meanwhile ωₙ^(k+m)=−ωₙ^k because ωₙ^m=−1. Therefore

Xₖ₊ₘ=Eₖ−ωₙ^k Oₖ.

Compute t=ωₙ^k Oₖ once, then return Eₖ+t and Eₖ−t. This two-output calculation is a butterfly. Predict the pair at k=2 in the eight-point example. Since ω₈²=−i, the product is −2i and the bins are X₂=1−2i and X₆=1+2i. A positive exponent would reverse these imaginary signs.

Experiment

Select Eight samples in Input samples to match the worked input (1,2,0,0,0,0,0,0). Before changing Root butterfly k, predict its two output positions k and k+n/2 and the shared product. Inspect k=2, then compare the other pairs with the exact table below. Changing the input preset returns k to zero. Then select Eight alternating samples and predict its nonzero bin before reading the result. Inspect all bins and the reconstruction and Parseval checks. The explorer uses bounded floating-point arithmetic; the symbolic table gives exact reference values for the Eight samples preset.

Follow one butterfly into two output bins

The forward transform uses exp(−2πijk/n), with no scale factor. The inverse uses the opposite sign and divides by n once. Predict the two bins before selecting a butterfly.

Input samplesFour samples: (1, 2, 0, −1)Eight samples: (1, 2, 0, 0, 0, 0, 0, 0)Eight alternating samples: (1, −1, …) Root butterfly k01 Reset explorer

Even-index inputs: (1, 0). Odd-index inputs: (2, -1). Indices start at zero. Each list is transformed recursively before this final combine step.

Even output E[0] = 1

Odd output O[0] = 1

Twiddle ω⁽ᵏ⁾ = 1; multiply once to get T = 1.

Upper branch: X[0] = E + T = 2

Lower branch: X[2] = E − T = 0

Each butterfly performs one complex multiplication and two complex additions/subtractions. The full recursive FFT used 4 multiplications and 8 additions/subtractions. Straight direct summation uses 16 multiplications and 12 additions under the stated dense convention. These are arithmetic counts, not measured timings; trivial twiddles are counted and root generation, storage and display work are excluded.

All bins in natural index order

Index	Input	FFT output	Direct DFT	Inverse reconstruction
0	1	2	2	1
1	2	1 − 3i	1 − 3i	2
2	0	0	0	0
3	-1	1 + 3i	1 + 3i	-1

Maximum FFT-versus-direct difference: 6.66e-16. Maximum inverse reconstruction error: 4.31e-16. Input energy Σ|x|² = 6; spectral energy Σ|X|²/n = 6.000000.

Complex values are rounded to six decimal places; tiny residual components can disappear in the table while remaining in the error measurements. The calculations use floating-point arithmetic. Agreement in these bounded cases checks the implementation; the lesson's induction proves the general algorithm.

Try the alternating eight-sample input. Predict which bin carries all the energy, then compare the k = 0 butterfly with the other root butterflies. Changing examples resets k to zero. These temporary controls do not record completion.

Compute the last column independently, then check it against the table. Each row uses Eₖ=1 and Oₖ=2 for the fixed input introduced above. Read √2 as the positive square root of two.

The right-hand output column contains bins 4,5,6,7, not bins 1,3,5,7. This distinction is an indexing contract, not a cosmetic display choice. Rearranging the two columns by visual alternation would mix frequencies.

Understand

The recursive algorithm accepts n=2^r for an integer r≥0. At n=1, return the single value: the defining sum has one term and no arithmetic operation is needed. Otherwise create the even and odd subsequences in increasing index order, recursively transform both, and combine matching child frequency k with the butterfly above. Assign the sum to output[k] and the difference to output[k+m]. Each output index occurs exactly once.

This proves correctness by induction. The base is the one-term definition. Assuming both child arrays are their correct transforms in increasing frequency order, the two identities prove that every parent output equals its defining n-point Fourier sum. Recursion terminates because the positive length halves to one. Nothing in this proof assumes real inputs, sparsity or a special waveform.

Track actual input indices to see why a recursive implementation does not need an unexplained output permutation. At length eight, the even child receives original indices 0,2,4,6; its even grandchild receives 0,4 and its odd grandchild receives 2,6. The odd child receives 1,3,5,7 and similarly splits into 1,5 and 3,7. Each child uses local indices starting at zero. Its result returns in natural local frequency order, and the parent places the combined bins in natural parent frequency order. Some iterative FFT implementations use an explicit bit-reversal permutation. That is an implementation choice requiring its own index explanation; it is not an additional unexplained step in this recursive algorithm.

Now count the stated algorithm, including products by trivial twiddles such as 1 and products whose data operand happens to be zero. A length-n call performs n/2 butterflies, each using one complex multiplication and two complex additions or subtractions. With A(1)=M(1)=0,

A(n)=2A(n/2)+n, and M(n)=2M(n/2)+n/2.

Every nontrivial recursion level totals n additions and n/2 multiplications across its calls. There are log₂n such levels, so A(n)=n log₂n and M(n)=(n/2)log₂n. At n=8 this is 24 additions and 12 multiplications. A direct DFT that computes every term product and uses n−1 additions per bin instead uses 64 multiplications and 56 additions. At n=1 our FFT returns its input directly, whereas that stipulated direct loop would still multiply its one term by 1; the counting conventions explain the boundary difference.

These are complex-arithmetic counts, not timings. They exclude computing or looking up twiddles, copying subsequences, indexing, memory traffic and recursion overhead. Counting a complex multiplication as one operation also hides its real arithmetic implementation. Special handling of trivial factors can reduce counts, but must be stated rather than silently mixed into one side of a comparison.

The inverse is obtained by conjugating the input spectrum, running this same forward algorithm, conjugating the result, and dividing every coordinate by n. Conjugation turns the negative exponential into the positive one, so this is exactly the previously derived inverse formula. It adds conjugation and scaling work outside the forward FFT counts. Omitting the division reconstructs n times the input, not the input itself.

Apply

Use a changed input with its only nonzero coordinate x₂=1 at length eight. The even child is (0,1,0,0), the odd child is zero. Derive the result X=(1,−i,−1,i,1,−i,−1,i) both recursively and from the single surviving term ω₈^(2k). This case tests the local-to-original index tracking that a constant input would not expose.

Then try xⱼ=(−1)^j. Root orthogonality predicts a single nonzero frequency bin, X₄=8, with all other bins zero. Inverting that spectrum divides by eight and recovers the alternating sequence. These exact cases can diagnose incorrect indexing, signs and scaling independently of another implementation.

For floating-point verification, compare against a separately coded direct DFT, test inverse reconstruction, and check Σ|Xₖ|²=nΣ|xⱼ|². Those checks are useful evidence, not an all-input error proof. Cancellation can make even a direct computation inaccurate: the exact zero-frequency sum of (10^16,1,−10^16,0) is 1, while left-to-right binary64 addition may lose the 1 and return zero. A different grouping can retain it. Use absolute as well as scale-dependent tolerances, especially at theoretically zero bins; never divide by a zero reference to form a relative error. A small disagreement requires investigating precision and accumulation order, not automatically declaring one implementation an exact oracle.

Continue

Record the sign convention, inverse normalization, supported lengths, recursive input order and output placement with any FFT implementation. A length-three DFT exists even though this radix-two algorithm does not handle it. Padding changes the transform length and sampled frequencies and must be an explicitly chosen operation.

The theorem here is exact algorithmic equality over complex numbers. The work recurrence describes one implementation convention. Finite numerical checks describe the tested inputs and arithmetic. Keeping these three statements separate makes both correctness proofs and debugging reports more useful.