0.10 Implementing ffts in practice  (Page 20/21)

Numerical accuracy in ffts

An important consideration in the implementation of any practical numerical algorithm is numerical accuracy: how quickly dofloating-point roundoff errors accumulate in the course of the computation? Fortunately, FFT algorithms for the most part haveremarkably good accuracy characteristics. In particular, for a DFT of length $n$ computed by a Cooley-Tukey algorithm with finite-precision floating-point arithmetic, the worst-case error growth is $O(logn)$ [link] , [link] and the mean error growth for random inputs is only $O\left(\sqrt{logn}\right)$ [link] , [link] . This is so good that, in practical applications, a properly implemented FFT will rarely be asignificant contributor to the numerical error.

The amazingly small roundoff errors of FFT algorithms are sometimes explained incorrectly as simply a consequence of the reduced number of operations: since there are fewer operations compared to a naive $O\left(n^2\right)$ algorithm, the argument goes, there is less accumulation of roundoff error. The real reason, however, is more subtle than that, and has to do with the ordering of the operations rather than their number. For example, consider the computation of only the output $Y\left[0\right]$ in the radix-2 algorithm of [link] , ignoring all of the other outputs of the FFT. $Y\left[0\right]$ is the sum of all of the inputs, requiring $n-1$ additions. The FFT does not change this requirement, it merely changes the order of the additions so as to re-use some of them for other outputs. In particular, this radix-2 DIT FFT computes $Y\left[0\right]$ as follows: it first sums the even-indexed inputs, then sums the odd-indexed inputs, then adds the two sums; the even- and odd-indexed inputs are summed recursively by the same procedure. This process is sometimes called cascade summation , and even though it still requires $n-1$ total additions to compute $Y\left[0\right]$ by itself, its roundoff error grows much more slowly than simply adding $X\left[0\right]$ , $X\left[1\right]$ , $X\left[2\right]$ and so on in sequence. Specifically, the roundoff error when adding up $n$ floating-point numbers in sequence grows as $O\left(n\right)$ in the worst case, or as $O\left(\sqrt{n}\right)$ on average for random inputs (where the errors grow according to a random walk), but simply reordering these n-1 additions into a cascade summation yields $O\left(\log n\right)$ worst-case and $O\left(\sqrt{\log n}\right)$ average-case error growth [link] .

However, these encouraging error-growth rates only apply if the trigonometric “twiddle” factors in the FFT algorithm are computedvery accurately. Many FFT implementations, including FFTW and common manufacturer-optimized libraries, therefore useprecomputed tables of twiddle factors calculated by means of standard library functions (which compute trigonometric constants to roughlymachine precision). The other common method to compute twiddle factors is to use a trigonometric recurrence formula—this savesmemory (and cache), but almost all recurrences have errors that grow as $O\left(\sqrt{n}\right)$ , $O\left(n\right)$ , or even $O\left(n^2\right)$ [link] , which lead to corresponding errors in the FFT. For example, one simplerecurrence is $e^{i(k+1)\theta }=e^{ik\theta }{e}^{i\theta }$ , multiplying repeatedly by $e^{i\theta }$ to obtain a sequence of equally spaced angles, but the errors when using this process grow as $O\left(n\right)$ [link] . A common improved recurrence is $e^{i(k+1)\theta }=e^{ik\theta }+e^{ik\theta }(e^{i\theta }-1)$ , where the small quantity In an FFT, the twiddle factors are powers of ${\omega }_n$ , so $\theta$ is a small angle proportional to $1/n$ and $e^{i\theta }$ is close to 1. $e^{i\theta }-1=cos\left(\theta \right)-1+isin\left(\theta \right)$ is computed using $cos\left(\theta \right)-1=-2{sin}^2(\theta /2)$ [link] ; unfortunately, the error using this method still grows as $O\left(\sqrt{n}\right)$ [link] , far worse than logarithmic.