The Fast Fourier Transform

Press Gauss, 3 periods of 4. These twelve numbers are from Gauss's manuscript: the declination of the northern edge of the band of sky the asteroid Pallas can wander in, every 30° of right ascension, in arcminutes, copied from von Zach's journal. He wanted the trigonometric series through them. Instead of one twelve-term sum he took every third value (0°, 90°, 180°, 270°; then 30°, 120°, …; then 60°, …), fitted each set of four on its own, and combined the three short fits. The dashed curves are those fits, and they disagree: four points cannot see the faster wiggles, so each set gets a different answer, which is what sent him looking for the combination. The table under the results is his calculation done again: the first period gives 779.5, −415.0, −726.0 and 43.5, the figures in the manuscript, and the combined series is the direct twelve-term one. He did it a second way, four periods of three, and got the same series. The index arithmetic is the one Cooley and Tukey published in 1965, as Heideman, Johnson and Burrus showed in 1984. What the manuscript leaves out matters as much. It never counts the work saved: where Gauss says the split is quicker, he leaves the proof to practice, praxis tentantem docebit, whoever tries it will see. Its worked examples have two factors, 12 = 3 × 4 and 36 = 6 × 6, and article 27 only says the split extends to more. It was not published in his lifetime, appearing in Latin in his collected works in 1866, and it is undated: the Werke's editor calls it a reworking of a draft apparently begun in October 1805, and Heideman, Johnson and Burrus place it in October or November 1805 from a set of clues, among them an eccentricity of Juno, 0.254236, that the manuscript uses and that was printed in May 1805.

Press 1,024 points. The direct sum multiplies every sample by every twiddle factor: N² = 1,048,576 multiplications. The FFT runs log₂N = 10 stages of N/2 = 512 butterflies with one multiplication each: 5,120. That is 204.8 times fewer, and the ratio is 2N/log₂N, so it keeps growing. Press A million points: 1,099,511,627,776 multiplications against 10,485,760, a ratio of 104,857.6. Double N and the direct sum's work goes up four times; the FFT's goes up 2.2 times at a thousand points, and closer to twice the further you go. At 8 points the saving is only 5.3 times, which is why the diagram looks like a lot of wiring for little gain. Both counts include the multiplications by W⁰ = 1, as Cooley and Tukey's own N² does; they counted multiply-and-add pairs and got fewer than 2N log₂N, and either way the ratio grows like N/log N. No computer changes either count, or their ratio: a faster one changes the timings, and not necessarily both by the same factor. The button under the tables times both in your browser; the measured ratio climbs with N too, though not to these figures, because memory traffic and additions cost time and the count ignores them.

The accuracy table runs both methods with every operation rounded to single precision, about seven significant digits, against a double-precision answer. At 1,024 points the direct sum is off by 5.75×10⁻⁷ of the answer's size and the FFT by 1.23×10⁻⁷, about 4.7 times less. Going from 16 points to 1,024, 64 times as many, made the direct sum's error 7.3 times larger and the FFT's 1.4 times. Each output of the direct sum is a running total of N rounded products, so its rounding errors pile up over 1,024 additions; each output of the FFT passes through only log₂N = 10 butterflies. At 16 points there is nothing in it, and on this input the FFT is even slightly worse, 8.74×10⁻⁸ against 7.91×10⁻⁸. Both methods get each twiddle factor the same way, computed in double precision and rounded once, so the table compares the two algorithms and not two ways of getting a sine.