B. Fast Algorithms for Discrete Fourier Transform

1
B. Fast Algorithms for Discrete Fourier
Transform B.1. Introduction (9.0, 9.1) B.2.
Radix-2 Decimation-in-Time FFT (9.3, 9.5) B.3.
Radix-2 Decimation-in-Frequency FFT (9.4,
9.5) B.4. FFT for a More General N (9.5) B.5.
Inverse Fast Fourier Transform
2
B.1. Introduction The discrete Fourier
transform (DFT) can be computed by using a class
of fast algorithms called the fast Fourier
transform (FFT). The principle of the FFT is to
decompose the DFT into shorter DFTs and butterfly
operations, which can be carried out more
efficiently. The FFT algorithms apply when the
length of the sequence is a composite integer
initially or after zero padding. B.2. Radix-2
Decimation-in-Time FFT B.2.1. Principle Let
us assume that x(n) is a finite-length sequence
over 0?n?N-1, where N2m (m?2) initially or after
zero-padding. The DFT of x(n) is defined as
(B.1)
3
x(n) is separated into two sequences. One
consists of even-numbered samples, and the other
consists of odd-numbered samples. Thus,
(B.2)
X(k) is halved. The first half is
(B.3)
4
The second half is
(B.4)
which is also written as
(B.5)
The radix-2 decimation-in-time FFT is based on
(B.3) and (B.5).
5
According to (B.3) and (B.5), the DFT can be
computed using the following algorithm (figure
B.1) (1) Find the DFT of N/2 the even-numbered
samples and the DFT of the N/2 odd-numbered
samples. (2) Compute the DFT of the original
sequence using (B.3) and (B.5). The above
decomposition is continued until each DFT has
only 1 point. The 1-point DFTs contain no
operation and do not have to be carried out.
The number of complex multiplications required in
this algorithm is N(log2N)/2. This number is less
than N2, the number of complex multiplications
required in the direct DFT. The following table
gives a few typical cases.
6
x(0)
X(0)
DFT (N/2)
x(2)
X(1)



x(N-2)
X(N/2-1)
x(1)
exp-j2?0/N
X(N/2)
DFT (N/2)
-
exp-j2?1/N
x(3)
X(N/21)
-



x(N-1)
exp-j2?(N/2-1)/N
X(N-1)
-
Figure B.1. Radix-2 Decimation-in-Time FFT.
7
B.2.2. Several Details This algorithm
consists of two steps the sequence reordering
and the butterfly operation.
The time sequence is reordered in the
bit-reversed way. That is, the sample at position
bm-1bm-2b0 (binary index) is transferred to
position b0b1bm-1 (binary index). The above
table shows the case of N8.
8
The butterfly operation can be carried out
using three nested loops or recursively. It is an
in-place operation. That is, in each stage of
computation, the output sequence and the input
sequence can use the same memory space. A few
variations of this algorithm can be obtained by
changing the orders of the input sequence and of
the output sequence. B.3. Radix-2
Decimation-in-Frequency FFT B.3.1. Principle
Let us assume that x(n) is a finite-length
sequence over 0?n?N-1, where N2m (m?2) initially
or after zero-padding. The DFT of x(n) is defined
as
(B.6)
X(k) is separated into two sequences. One
consists of even-numbered
9
samples, and the other consists of odd-numbered
samples. Thus,
(B.7)
and
(B.8)
x(n) is halved. (B.7) becomes
(B.9)
which is also written as
10
(B.10)
and further
(B.11)
(B.8) becomes
(B.12)
11
which is also written as
(B.13)
and further
(B.14)
The radix-2 decimation-in-frequency FFT is
developed according to (B.11) and (B.14).
12
According to (B.11) and (B.14), the DFT can
be computed using the following algorithm (figure
B.2) (1) Find sequences x(n)x(nN/2) and
x(n)-x(nN/2)exp(-j2?n/N), where 0?n?N/2-1. (2)
Find the DFTs of the two sequences. The above
decomposition is continued until each DFT has
only 1 point. The 1-point DFTs contain no
operation and do not have to be carried out.
The number of complex multiplications required in
this algorithm is N(log2N)/2. This number is less
than N2, the number of complex multiplications
required in the direct DFT. The following table
gives a few typical cases.
13
x(0)
X(0)
DFT (N/2)
x(1)
X(2)



x(N/2-1)
X(N-2)
exp-j2?0/N
x(N/2)
X(1)
-
DFT (N/2)
exp-j2?1/N
x(N/21)
X(3)
-



exp-j2?(N/2-1)/N
x(N-1)
X(N-1)
-
Figure B.2. Radix-2 Decimation-in-Frequency FFT.
14
B.3.2. Several Details This algorithm
consists of two steps the butterfly operation
and the sequence reordering. The butterfly
operation can be carried out by three nested
loops or recursively. It is an in-place
operation. That is, in each stage of computation,
the output sequence and the input sequence can
use the same memory space.
15
The frequency sequence is reordered in the
bit-reversed way. That is, the sample at position
bm-1bm-2b0 (binary index) is transferred to
position b0b1bm-1 (binary index). The above
table gives the case of N8. A few variations
of this algorithm can be obtained by changing the
orders of the input sequence and of the output
sequence. B.4. FFT for a More General N When
the length of the sequence is a composite integer
initially or after zero padding, the FFT
applies. Assume that x(n) is a finite-length
sequence over 0?n?N-1, where NN1N2 (N1 and N2
are two positive integers not equal to 1)
initially or after zero-padding. The DFT of x(n)
is defined as
(B.15)
16
Letting nn1N2n2, 0?n1?N1-1, 0?n2?N2-1 and
kk2N1k1, 0?k1?N1-1, 0?k2?N2-1, we obtain
(B.16)
which is further written as
(B.17)
17
Changing the order of the summations, we obtain
(B.18)
According to (B.18), the DFT can be computed
by the following algorithm (1) For each n2,
compute the DFT of x(n1N2n2) with respect to n1
to obtain p(k1,n2). (2) Multiply p(k1,n2) by
exp(-j2?k1n2/N) to obtain q(k1,n2). (3) For
each k1, compute the DFT of q(k1,n2) with respect
to n2 to obtain X(k2N1k1).
18
The number of complex multiplications
required in this algorithm is N(N1N21). This
number is less than N2, the number of complex
multiplications required in the direct DFT.
The above decomposition can be continued and the
computation can be further reduced until each DFT
has a prime number of points. B.5. Inverse Fast
Fourier Transform A class of fast algorithms
called the inverse fast Fourier transform (IFFT)
can be used to compute the inverse discrete
Fourier transform (IDFT). The DFT and the
IDFT are defined as follows
(B.19)
(B.20)
19
Comparing (B.19) and (B.20), we see that the
algorithms for the DFT can also be used to
compute the IDFT if we replace exp(-j2?kn/N) by
exp(j2?kn/N) and divide the resulting sequence by
N. x(n) can be written as
(B.21)
According to (B.21), the IDFT can be computed
using the following algorithm (1) Determine
X(k)/N. (2) Compute the DFT of X(k)/N. (3)
Conjugate the DFT of X(k)/N. x(n) can be
written as
(B.22)
where X?(k) is X(k) extended with period N.
According to (B.22), the
20
IDFT can be computed using the following
algorithm (1) Determine X?(-k)/N, 0?k?N-1. (2)
Compute the DFT of X?(-k)/N, 0?k?N-1.