74 CHAPTER 1. ANALYSIS OF DISCRETE-TIME LINEAR TIME-INVARIANT SYSTEMS
1.4 Fast Fourier Transform (FFT) Algorithm
Fast Fourier Transform, or FFT, is any algorithm for computing the N-point DFT with
a computational complexity of O(log N). Tt is not a new transform, but simply an
efficient method of ealeulating the DFT of x(n)
If we assume that N is even, wo can write the N-point DFT of x(n) as
XO) = y a(nje a"
fee Medes SE
+ y a(n)eF"
mode neBE LEO Ad
= r(QmnoF 2 4S a(2t 4 1H CHY (1.31)
mmo ay
We make the following substitutions:
ra(m) = (2m), where m=0,
ait)
(21-41), where f=
Rewriting Eq. (1.31), we get
q- Behn gh
XM) = SY ayime FP" 40 ae
im
xD ay + FF x! (1.32)
where X(7) (k) is the ¥-point DP’ of the even-numbered samples of x(n) and X{?) (k)
is the ¥-point DFT of the odd-numbered samples of 2(n). Note that both of them are
3 periodic discrete-time fictions.
We have the following algorithm to compute X(®)(k) for k
(N-)e
1. Compute Xj (A) for k = 0,
2. Compute X{7)(k) for k=0,--- Y= 1,
Perform the computation (1.32) with NV’ complex multiplications and N’ complex
ditions.
Actually, it is possible to use fewer than N complex multiplications. LetSec 14, Fast Fourier Transform (FFT) Algorithm 75
20) xo xPo+wex{Po
2)
DFr
2-2)
Total smbor of complex operations
Ca] Swen reson
Figure 1.97, The recursive implementation of the FPT supposing that N = 24, ‘There is a to-
tal of M = loz, N’ stages of computation, each requiring 3NV complex operations. Hence, the total
‘computational complexity is O(WV Tox N).
1, For large N, the FFT is much faster than the direct application of the definition
of DFT, which is of complexity O(N?),
2. ‘The particular implementation of the FFT deseribed above is called decimation.
in-time radia-2 F!
2. The number of operations required by an FFT algorithm ean be approximated
as CN log N, where C is a constant. There are many variations of FFT aimed at
reducing this constant eg., if N= 3", it may be better to uso a radix-3 FFT,
4, Note that
{yprreen} =
which is the IDF of x(n). ‘Thus, the FFT can also be used to compute the IDFT.Sec 14, Fast Fourier Transform (FFT) Algorithm 17
“9
sow vera( 0) pan orn ( 3 }
=i)
20)
(0) X®(0)
2(4) X®(1)
(2) X)(2)
(6) mer X®(3)
ini 2
pout DET of ( 21)
(1) x(4)
(5) X@(5)
(3) X®(6)
a(7) X®(7)
” (0)
ovine er ot (2) apie om 2 }
Figure 1.38, The & point PFT.
Example 1.26. The &-point FFT is depicted in Fig. 1.38. The values of the twiddle
factors are:
&
ion78 CHAPTER 1. ANALYSIS OF DISCRETE-TIME LINEAR TIME-INVARIANT SYSTEMS
&] = A) x
Nx1 Nx1
Ey
202)
(v2)
=)
23)
2(V=1)
Figure 1.39. ‘The FFT reduces the number of operations required to ealeulate the DFT by reducing
A®) to two AU!) that is only half the size of A). This operation is repeated with every recursion
‘until we reach the I-poine DPT.Sec 14, Fast Fourier Transform (FFT) Algorithm 79
Recall that the DF'T is a matrix multiplication (Pig. 1.35). One stage of the FFT
ly the multiplication by an N x A’ matrix to two multiplications by
matrices. This reduces the number of operations required to calculate the DFT
by almost a factor of two (Fig. 1.39).
Another interpretation of FFT involves analyzing the matrix
hype lh (1 o%
SEEM eo J
where and L are nonnegative integers such that & <2! Note that
(Aazy)* (Aix)
(Aix Arty)
nultiplication by Ag, preserves distances and angles — ronghly speaking, it is a
1 or reflection. Continuing the matrix decomposition of Fig. 1.39 further until
wwe get the full FFT, it can be shown that FFT consists of log NV multiplications by
2x2 matrices of the form y2Aq,., each operating on a pair of coordinates.2 Therefore,
FFT breaks down the multiplication by the DPT matrix A into elementary planar
transformations,
1.4.1 Fast Computation of Convolution
Consider a Tinear system described by
Sx, (1.33)
where x is the Nx 1 input veetor, representing an N-periodie input signal; $ is an Nx
matrix; aud y is the N x 1 output vector, representing an N-periodie output signal,
What conditions must the matrix S satisfy in order for the system to be time-invariant,
i.0., invariant to cireular shifts of the input vector?
Note that a cireular shift by one sample is
(0) r(-1) =2(N —1)
a(1) x(0)
(2) = a(1)
a(N—1) a(N
)
"The same conclusion can be reached by examining an FFT diagram such as Fig. 1.38.80 CHAPTER 1. ANALYSIS OF DISCRETE-TIME LINEAR TIME-INVARIANT SYSTEMS
Let the first column of S be
n(0)
n(1)
h(2)
(N -1)
Note that when
1
0
0], then y =n,
0
‘and when
0
1
x-| 0
0
then y is the second column of S, which therefore, in order for S to be
circular shifts, must be equal to:
MN ~1)
h(O)
AC)
Similarly, when
6 |. then y is the third column of S, ete.
0,
‘Thus, the matrix S must have the following structure:
0) -A(N 1) K(N-2) KL)
A(1) 10) ACN 1) (2)
n(2) Al) (0) H(3)
AN =1) M2) AN=3) = W(0)Sec 14, Fast Fourier Transform (FFT) Algorithm 81
‘This is called a eéreulant matréx, We ean then write Eq. (1.33) as
wa
y(n) = > x(m)h(n—m)
mo
xa
= ¥ x(m)n((n—m) moa N) (1.34)
mo
= r@b(n)=2@h (1.35)
Eq, (1.35) is called a circular convolution or a periodic convolution. Note that formula,
(1.34) works even when + or fare non-periodie. Observe the following:
«© For y(0), the sum of the indices of 2x and f is always 0’ mod N for every term,
y(0) = 2(0)h(O) + 2(1)A(N — 1) + 2(2)K(N — 2) +--+ (NW — 1)A(1)
© For y(1), the sum of the indices of zr and h is always I mod N for every term,
y(l) = 2(O)A(1) +-2(1)R(O) + x2)A(N = 1) ++ +2(W — AQ)
isis true for all y(k), k= 0,1,-+-,N— 1.
‘What are the eigenvectors of 5? Let us try
e=| xe™ » Where k= 0,1,
Wo have:
Ain) ® we
N=
= SY hlon)ge(n — m)
mb
Naa
-»> hm) heron)
{x yew
= HW, pe
DET of a82 CHAPTER 1. ANALYSIS OF DISCRETE-TIME LINEAR TIME-INVARIANT SYSTEMS
have that
SB = Hk) Be
where ge is the k-th eigenvector and H(k) gives the corresponding eigenvalue. ‘There-
fore,
(0) o
si A)
S( go at gv-1)=(g0 Bt gy )
o
Torr main B H(N-1)
‘Then $ can be written as
H(0)
0
na
A
o
H(N~1)
where the DPT matrix A is:
si
A=NBH *
8X1
Complex exponentials are the cigenveetors of cirenlant matrices. They diagonalize
circulant matrices. ‘Thus, for any x € CN,
H(0)
H(1)
H(N 1)
‘Let ws compare two algorithms for computing the circular convolution of x and hh.
Algorithm 1 Directly perform the multiplication Sx. ‘This has computational eom-
plexity © (2).
Algorithm 2 — 1. Represent x in the eigenbasis of S, ie., the Fourier basis,
X=Ax.
This stop can be done with FFT whose complexity is O(N log N),Sec 14, Fast Fourier Transform (FFT) Algorithm 83
Step 1 Step 2 Step 3
N-poiat DPT
a(n) > X(k) potas ET
Y(k) = X(k)A(k) | ¥(K) — y(n)
N-point DPT:
h(n) > H()
Figure 1.40, An illustration of the FFT implementation of the circular convolution,
2. Compute the representation of y in the eigenbasis of S:
(0)
WN 1)
‘This computation has complexity O(N).
3. Reconstruct y from its Fourier coefficients:
y= BY.
‘This has complexity O(N log N), if done using the FET.
‘This algorithm is summarized in Fig, 1.40, Its total complexity is O(V loz N).
(Note that the second algorithm does not necessarily perform better for any matrix.)
Example 1.27. This example explores the relationship between the convolution and the
cireular convolution. Let x and h be N-periodie siguals, and let
_ f[2t), 0