FFT
A Fourier transform of N samples, written the way the definition is written, builds every one of N outputs from every one of N inputs. That is N-squared complex multiplications. Cooley and Tukey's 1965 paper is the one that made a divide-and-conquer transform standard for machine computation, and it is not the first time the factorisation was found: the standard survey of that prehistory is titled after Gauss, who worked out the same factorisation decades before Fourier's own work was published and whose version went unnoticed. The date on this page is a publication date, and it is worth reading as one. Cooley and Tukey's paper is not a new transform. It is the observation that the sums overlap: split the samples into the even ones and the odd ones, transform each half, and the two halves can be recombined into the whole with one multiplication per pair. Do that all the way down and the work falls to N/2 multiplications per halving. This page runs both routes over the same samples and counts what each one spends, and then asks the question a diagram cannot: do the two answers agree?
The same answer, and the bill for each way of getting it
1 A signal whose transform is written down in advance, so the answer owes nothing to this page
Pick a signal and how many samples. Every signal here has a transform anybody can write down without running either route, which is what makes it a check rather than a demonstration.
every sample the same, so all the energy is in the zero bin. Written down in advance, its transform is N in the zero bin and nothing anywhere else.
| n | sample |
|---|---|
| 0 | 1 |
| 1 | 1 |
| 2 | 1 |
| 3 | 1 |
| 4 | 1 |
| 5 | 1 |
| 6 | 1 |
| 7 | 1 |
2 The direct sum, every output built from every input, with the multiplications counted
MeasuredEvery output from every input
This is equation one of the paper, done as written. The count beside it was incremented once per complex multiplication while it ran.
64 complex multiplications
| bin | value | size |
|---|---|---|
| 0 | 8 | 8 |
| 1 | 0 | 0 |
| 2 | 0 | 0 |
| 3 | 0 | 0 |
| 4 | 0 | 0 |
| 5 | 0 | 0 |
| 6 | 0 | 0 |
| 7 | 0 | 0 |
Each of the 8 outputs took 8 complex multiplications, so this cost 64 of them and 64 additions.
3 The same transform split even from odd, one column of butterflies per halving, counted again
MeasuredThe same thing, halved and halved again
Each row is one column of butterflies. A butterfly takes a pair, multiplies one of them by a twiddle factor and produces their sum and their difference, so one multiplication serves two outputs. The values shown are the array after that column has run.
12 complex multiplications, a 5.3x saving
| column | span | butterflies | the array after it |
|---|---|---|---|
| 1 | 2 | 4 | 2, 0, 2, 0, 2, 0, 2, 0 |
| 2 | 4 | 4 | 4, 0, 0, 0, 4, 0, 0, 0 |
| 3 | 8 | 4 | 8, 0, 0, 0, 0, 0, 0, 0 |
3 columns of 4 butterflies each: 12 complex multiplications and 24 additions, against 64 multiplications for the direct sum.
4 The two answers side by side, and the largest place they disagree
Bin by bin, and where they part
The marked rows are the bins where the two routes returned different numbers.
| bin | direct sum | butterfly | apart by |
|---|---|---|---|
| 0 | 8 | 8 | none |
| 1 | 0 | 0 | 3.3e-16 |
| 2 | 0 | 0 | none |
| 3 | 0 | 0 | 3.3e-16 |
| 4 | 0 | 0 | none |
| 5 | 0 | 0 | 2.2e-16 |
| 6 | 0 | 0 | none |
| 7 | 0 | 0 | 1.1e-16 |
The two routes disagree. The worst is bin 1, by 3.33e-16, which is the cost of the twiddle factors not being exact. The closed form for this signal sits 0.00e+0 from the butterfly's answer.
Two routes to one answer is the sort of claim that gets rounded off to a claim that the two are identical. It is not identical past four samples, and the reason is worth more than the tidier sentence: the numbers a transform multiplies by are points on a circle, and a quarter of the way round is the only place where those points are whole numbers. Everywhere else the machine keeps an approximation, the two routes reach the same bin by different chains of them, and they land a few parts in a quadrillion apart. At four samples every factor is a quarter turn, and there the answers match exactly.
These ran in this browser when the page loaded. Each claim, whether it held, and the number behind it.
| claim | held | measured |
|---|---|---|
| eight points direct costs 64 complex multiplications, by butterfly 12 | yes | counted while running: 64 against 12, which is N squared against N/2 times log2 N |
| and at every size it stays under the 2N log2 N the paper promises | yes | checked at N = 4, 8, 16, 32, against 2N log2 N |
| the direct sum's cost is N squared at every size on the page | yes | 4 -> 16, 8 -> 64, 16 -> 256, 32 -> 1024 |
| every signal here has a transform written down in advance, and the butterfly finds it | yes | 16 signal-and-size pairs against their closed forms |
| at four points the two routes agree bit for bit | yes | all 4 signals, every twiddle a quarter turn and written down exactly |
| past four points they do not, and the whole disagreement is under a part in a trillion | yes | worst gap 2.89e-15 on constant at 32, bin 31 |
| the energy in the samples equals the energy in the bins, divided by N | yes | Parseval, in the unnormalised convention this page uses: the sum of the squared samples equals the sum of the squared bins over N. Checked for all 4 signals at sixteen points |
| transform it and transform it back and the samples return | yes | all 4 signals at sixteen points, within a part in a trillion |
What is real here, and what is not
The counts are counted, not worked out
Both routes increment a counter once per complex multiplication as they run, and the page prints that counter. Nothing here evaluates N squared or N/2 log N for display. The two are checked against each other in the table below, which is the point: if the implementation stopped matching the formula, the page would say so rather than keep printing the formula.
No twiddle is skipped, even the ones that are 1
The first butterfly in every column multiplies by 1, and a tuned implementation would skip it. This one does not, because skipping it would turn the count into a claim about how clever this file is rather than about the algorithm. The paper's bound of 2N log2 N is generous enough to hold either way, and the page checks against the bound rather than against a tidy formula.
Quarter turns are written down, not asked of Math.cos
Math.cos(Math.PI/2) returns 6.123233995736766e-17. Every twiddle factor at a multiple of a quarter turn is therefore a written-down 1, -1 or i here instead of a call. That is what makes the four-point claim true rather than nearly true, and it is also the reason the four-signal set is built from ones, zeroes and minus ones.
The two routes do not agree, and the page says so
Past four samples the direct sum and the butterfly return different numbers, worst case a few parts in a quadrillion. That is not a bug and it is not a rounding of the truth: they are different chains of approximate multiplications reaching the same bin. Claiming they match would be the one dishonest sentence available on this page, so the disagreement is measured and displayed instead.
The witness is a closed form, not the other route
A constant, a single spike, a cosine at a quarter of the sampling rate and those last two added together all have transforms that can be written down without running anything. Every claim about correctness here is against those, not against agreement between the two routes: two implementations that share a bug agree perfectly while both are wrong. Parseval is checked in the unnormalised convention this page uses, where the sum of the squared samples equals the sum of the squared bins DIVIDED BY N; the row that reports it says so rather than implying the two sums are equal.
Radix-2 only, and real inputs only
The sizes offered are powers of two and the samples are real. The 1965 paper is more general than that, covering any highly composite N, and real production transforms exploit the input being real to halve the work again. Neither is modelled. What is here is the case the paper singles out for a binary computer.
Sound: no
The temptation is obvious, since the subject is signals. What this page measures is an operation count and a rounding difference, and neither is a duration. Playing the signal would be a mood, not the measurement, and the studio's rule is that a sound has to carry the measurement or stay out.
1965 is a publication date, not the invention of the idea
The roster dates this page to the April 1965 issue of Mathematics of Computation, and that is the paper that made the algorithm general knowledge. The factorisation itself is older, and the recognised survey of that prehistory is Heideman, Johnson and Burrus, Gauss and the History of the Fast Fourier Transform, IEEE ASSP Magazine, October 1984. That survey is cited here for its existence and its title rather than quoted, because IEEE serves no readable text for it. An outside reader flagged the first draft of the lede for saying the direct sum is what a computer did until 1965, which is a claim about every computer before that year and is not supported by anything here.
Sources
- Cooley and Tukey, An Algorithm for the Machine Calculation of Complex Fourier Series, Mathematics of Computation 19(90), 1965 — the paper, including the operation bound this page checks itself against and the note that it was received in August 1964.
- Crossref's record for the paper's DOI, which is where the volume, the issue and the page range on this site come from rather than from memory.
- Crossref's record for Heideman, Johnson and Burrus, Gauss and the History of the Fast Fourier Transform, IEEE ASSP Magazine 1(4), October 1984 — the survey of what came before 1965. IEEE serves no text for it, so it is cited by its record rather than read here.
- Logical Art, the studio this belongs to.