Computational Physics

600 Appendix A

From now on, we shall use the notationf(xj)=fjandfˆ(xn)=fˆn. We assume
that the number of pointsNis equal to a power of 2:N = 2 m. If this is not the
case, we put the function forj-values fromNupward to the nearest power of 2 to
zero and do the transform on this larger interval. The Danielson–Lanczos lemma is
simple. The pointsxjfall into two subsets: one withjeven and one withjodd. The
Fourier transforms of these subsets are calledˆf 2 | 0 andfˆ 2 | 1 This notation indicates
that we take the indicesjmodulo 2 and check if the result is equal to 0 (jeven) or
1(jodd). We obtain

fˆn=

N∑− 1

j= 0

fje^2 πinj/N

=

N∑/ 2 − 1

j= 0

f 2 je^4 πinj/N+e^2 πin/N

N∑/ 2 − 1

j= 0

f 2 j+ 1 e^4 πinj/N

=fˆ 2 | 0 +e^2 πn/Nfˆ 2 | 1. (A.139)

The Fourier transform on the full set ofNpoints is thus split into two transforms
on sets containing half the number of points. If we were to calculate these sub-
transforms using the direct method, each would take(N/ 2 )^2 steps. Adding them as
in the last line of (A.139) takes anotherNsteps, so in total we haveN^2 / 2 +Nas
opposed toN^2 operations if we had applied the direct method to the original series.
Where does the gain in speed come from? Compare the transforms fork 0 andkN/ 2
respectively. For any pointxjwithjeven, we have the exponential exp( 2 πinj/N)
which are equal for these twok-values, but in the direct method, these exponentials
are calculated twice! The same holds for oddj-values and for other pairs ofkn-values
spaced bykN/ 2. Indeed, the two sub-transformsfˆ 2 | 0 andfˆ 2 | 1 are periodic with period
N/2 instead ofN– therefore their evaluation for allj=0,...,N−1 takes only
(N/ 2 )^2 steps, and the total work involved in finding the Fourier transform is reduced
approximately by a factor of two. That is not yet theNlog 2 Npromised above. But
this can be found straightforwardly by applying the transform of (A.139) to the
sub-transforms of lengthN/2. One then splits the sum into four sub-transforms of
lengthN/4. The reader can verify that the new form becomes

fˆn=ˆf 4 | 0 +e^4 πikn/Nfˆ 4 | 2 +e^2 πikn/Nfˆ 4 | 1 +e^6 πikn/Nˆf 4 | 3. (A.140)

The same transformation can be performed recursivelymtimes (rememberN= 2 m)
to arrive atNlog 2 Nsub-transforms of length 1, which are trivial.
It will be clear that the FFT can be coded most easily using recursive program-
ming. When recursion is avoided for reasons of efficiency, a bit more bookkeeping
isrequired:seeRefs.[ 1 , 23 , 24 ].