Page 462 - A First Course In Stochastic Models

D. THE DISCRETE FAST FOURIER TRANSFORM 457

as the sum of two discrete Fourier transforms, each of length n/2. Suppose we
know the c k and wish to compute the f ℓ from (D.1). It holds that

1 1
n−1 2 n−1 2 n−1
2πikℓ/n 2πiℓ(2k)/n 2πiℓ(2k+1)/n
c k e = c 2k e + c 2k+1 e
k=0 k=0 k=0
1 n−1 1 n−1
2 2
2πikℓ/(n/2) ℓ 2πikℓ/(n/2)
= c 2k e + w c 2k+1 e . (D.4)
k=0 k=0
The discrete Fourier transform of length n can thus be written as the sum of two
discrete Fourier transforms each of length n/2. This beautiful trick can be applied
recursively. For the implementation of the recursive discrete FFT procedure it is
convenient to choose

n = 2 m

for some positive m (if necessary, zeros can be added to the sequence f 0 , . . . , f n−1
m
in order to achieve that n = 2 for some m). The discrete FFT method is numer-
ically very stable (it is a fast and accurate method even for values of n with an
order of magnitude of a hundred thousand). The discrete FFT method that calcu-
lates the original coefficients f j from the Fourier coefficients c k is usually called
the inverse discrete FFT method. Ready-to-use codes for the discrete FFT method
are widely available. The discrete FFT method is a basic tool that should be part
of the toolbox of any applied probabilist. It is noted that the discrete FFT method
can be extended to a complex function defined over a multidimensional grid.

Numerical inversion of the generating function
Suppose an explicit expression is available for the generating function

∞
ℓ
P (z) = p ℓ z , |z| ≤ 1.
ℓ=0
How do we obtain the unknown probabilities p ℓ ? Choose an integer n = 2 m
such that
∞

p j ≤ ǫ
j=n

for some prespecified accuracy number ǫ, say ǫ = 10 −12 (often one can find a

∞
∞
known distribution {a j } such that j=k p j ≤ j=k j for all k; otherwise, the
a
truncation integer n has to be found by trial and error). Then calculate the complex