Transcript Document
Chapter 8
Fast Fourier Transform (FFT)
Section 9.1-9.4, 10.1
• 8.1 Introduction
• 8.2 Goertzel Algorithm
• 8.3 Fast Fourier Transform (FFT)
• 8.4 Inverse FFT (IFFT)
• 8.5 Signal Analysis using DFT
1
8.1.1 DFT Expression
• Consider how complicated this is when x[n] is complex.
• Represent x[n]=Re{x[n]}+j Im{x[n]}.
• We have the expression:
N 1
X [k ] x[n ]WNkn
(0 k N 1)
n 0
N 1
X [k ] x[n]WNkn
n 0
N 1
1
x[n ]
N
WN e
j
N 1
X [k ]W
k 0
kn
N
(0 n N 1)
2
N
X [k ] Re( x[n]) Re(WNkn ) Im(x[n]) Im(WNkn )
n 0
j Re( x[n]) Im(WNkn ) Im(x[n]) Re(WNkn )
2
8.1.2 DFT Computational Complexity
• For each K we have 4N real multiplications, 4N-2 real adds.
• For all N of X[k] we have
– N(4N) = 4N2 real multiplication
– N(4N-2) = 4N2 -2N real adds
• Example. N=1000
– # multiplications
= 4,000,000 !
– # additions
4,000,000
• The computational complexity is O(N2). As N increase, the
order becomes huge!
• Storage Requirements (memory): x[n] -N values (complex),
WN - N values (complex), which can be reduced due to
periodicity and symmetry.
3
8.1.3 DFT Symmetry and Periodicity Properties
• Symmetry and Periodicity
k ( N n )
N
kn
N
W
W
W
kn *
N
• Periodic in both k and n with period N
k ( n N )
N
W W
kn
N
(k N )n
N
W
W86 ,W814
Example N=8
W84 ,W812
W87 ,W815
W81,W89
4
8.1.4 Symmetry and Periodicity Properties
•
Based on the symmetry, we can write for real part:
Re{x[n]}Re{WNkn }+ Re{x[N-n]}Re{WNk(N-n) }
Group terms, drop one multiplication we get:
(Re{x[n]}+ Re{x[N-n]})Re{WNkn}
Since Re{WNkn }= Re{WNk(N-n) }
•
Similarly for other three terms, and this results in a reduction
by 50% of real multiplications without really doing anything.
•
However, the overall complexity is still O(N2)
•
The periodicity property is better that leads to the Fast
Fourier Transform (FFT) with O(Nlog(N))
5
8.1.5 FFT History
• Known for many years
– Gauss 1805
– Runge 1905
– Danielson and Lanczos (1942)
• Hand calculation of the DFT
• Relatively short sequence
• Cooly & Tukey – 1965
– Modern algorithm that could decompose a DFT into a
sum of shorter DTFs
– When N is a composite number, i.e., the product of two or
more integers.
6
8.1.6 FFT Principles
• Decompose a sequence into successively smaller sequence to
reduce overall computation .
• Example:
– Consider N=100
– Brute Force O(N2)=10000
– If we can break this into 2 (50 point) DFTs, then: O(502) + O(502)
=5000<10000 !
– Decompose further, gaining each time.
• Two basic classes
– Decimation in time (DIT)
– Decimation in frequency (DIF)
7
8.2.1 Goertel Algorithm: Periodicity
• The Goertel algorithm is an example of how the periodicity of
the sequence WNkn can be used to reduce computation.
WNkN e j ( 2 / N ) Nk e j 2k 1
kN
N
X [k ] W
N 1
x[l ]W
r 0
kr
N
N 1
x[l ]WNk ( N r )
r 0
• To suggest the final result, let us define the sequence
yk [ n ]
x[r]W
r
k ( N r )
N
u[n r ].
X [k ] yk [n] nN
• It can be interpreted as a discrete convolution of the finite
duration sequence x[n],0 n N 1 with the sequence WNknu[n].
• Consequentially, yk [n] can be viewed as the response of a system
kn
with impulse response WN u[n]. to a finite-length input x[n]. In
particular, X[k] is the value of the output when n=N.
8
8.2.2 Goertel Algorithm: Complexity
xe [n]
yk [n ]
xe [ N ] 0
yk [1] 0
WNk
Z-1
• For a particular value of k, we need 4N real multiplications and
4N real additions to compute X[k].
• This procedure is slightly less efficient than the direct method.
• It avoids the computation or storage of the coefficient WNkn,
since these quantities are implicitly computed by the recursion.
9
8.2.3 Goertel Algorithm: Complexity Reduction
1
H k ( z ) (W z )
1 WNk z 1
n 0
k 1 n
N
(1 WNk z 1 )
(1 WNk z 1 )
Hk ( z)
k 1
k 1
(1 WN z )(1 WN z ) 1 2 cos(2k / N ) z 1 z 2
• 2(N+2) real multiplications
• 4(N+1) real additions
10
8.2.4 Goertel Algorithm: Comments
• In either the direct method or the Goertzel method, we do not
need to evaluate X[k] at all N values.
• We can evaluate X[k] for any M values of k, with each DFT
value being computed by a recursive system.
• This is slightly less efficient than DFT performed directly.
• But, WNkn does not need to be computed and stored in
advance due to recursive nature of the algorithm.
• Advantage: Goertzel’s algorithm can be used to compute
efficiently a small number of DFT coefficients.
• Next – the real FFT
11
8.3.1 DIT-FFT: Concepts
• Decimation in Time (DIT) FFT – Classical approach:
– Based on decomposing x[n] into small sequences
– Exploits both periodicity and symmetry of W
– Usually consider N=2v so that we can successively subdivide
by factor of 2 (v times) all the way down to length – 2
subsequences.
• Divide into even and odd indexed points
Let N=2v
x[n ] xeven[n ] xodd [n ]
N 1
X [k ] x[n ]WNkn
n 0
X [k ]
kn
x
[
n
]
W
N
n odd
kn
x
[
n
]
W
N
n even
12
8.3.2 DIT-FFT: Formulas
So we can write
X [k ]
N 21
x[2n]W
k ( 2n )
N
n 0
N 21
x[2n 1]W
k ( 2 n 1)
N
n 0
Want to write this as a sum of N/2 points DFT
X [k ]
N 2 1
x[2n](W )
2 kn
N
n 0
X [k ]
N 2 1
x[2n](W
2 kn
N
n 0
WN e
j 2N
)
W e
2
N
N 2 1
k 2n
k
x
[
2
n
1
]
(
W
)
W
N
N
n 0
WNk
j 2N 2
e
N 2 1
x[2n 1](W
2 kn
N
n 0
j 4N
e
j N2 / 2
)
WN
2
13
8.3.2 DIT-FFT: Formulas (Cont’d)
• Rewrite
X [k ]
N 2 1
x[2n]WN 2 W
kn
n 0
k
N
N 2 1
kn
x
[
2
n
1
]
W
N
2
n 0
X [k ] G[k ] H [k ]
X [k ] G[k ] WNk H [k ]
0 k N 1
• Where G[k] and H[k] are N/2 point DFT
• G[k] and H[k] are both periodic with period N/2
• So, must be computed over only N/2 points (not N)
• X(k) is still N points long.
• Lets check the flow diagram
14
8.3.3 DIT-FFT: Flow Graph
X [k ] G[k ] WNk H [k ]
X [0] G[0] W80 H [0]
X [1] G[1] W81H [1]
X [4] G[0] W84 H [0]
X [5] G[1] W85 H [1]
X [7] G[3] W87 H [3]
Because of periodicity
G[4]=G[0]
H[4]=H[0]
G[5]=G[1]
H[5]=H[1]
G[6]=G[2]
H[6]=H[2] ….
15
8.3.4 DIT-FFT: Reduced Complexity
• Consider computation:
– Recall that DFT is O(N2)
– Now length is N/2, not N, but there are two N/2 point DFTs
– At first (Complex)
– N2 multiplications and adds (each)
– Now, 2*(N/2)2 = N2/2 multiplications and additions (each)
– Plus N complex multiplications for W
– Plus N complex additions
• Total = N+N2/2 complex multiplications and adds.
• Since N+N2/2 < N2 ! A savings for N=2v.
• We can continue this decomposition process and get more
savings, if N/2 is even and greater than 2.
16
8.3.5 DIT-FFT: Procedure
• Assume N=8
– Stage I : Form two 4 point DFTs, plus combining algebra
– Stage II: Split 4-point DFTs into four 2-point DFTs, plus combining
algebra
– Would continue if N>8
• Consider the two 4-point DFTs
G[ k ]
N 2 1
kn
g
[
n
]
W
N 2
wit h g[n ] x[2n ]
n 0
G[ k ]
N 4 1
g[2l ]W
l 0
G[ k ]
N 4 1
g[2l ]W
l 0
H [k ]
2 lk
N 2
lk
N 4
N 4 1
h[2l ]W
l 0
lk
N 4
N 4 1
g[2l 1]W
( 2 l 1) k
N 2
l 0
N 4 1
W
k
N 2
W
k
N 2
g[2l 1]W
l 0
lk
N 4
N 4 1
lk
h
[
2
l
1
]
W
N 4
l 0
17
8.3.5 DIT-FFT: Procedure (Cont’d)
• Where h[n] x[2n 1]
• This represents the 4-point DFTs as four 2-point DFTs plus 2
stages of combining algebra
• If N=2v, it will take v stages to do this completely.
• log2N decomposition to get down to 2-point DFTs, now require
4(N/4)2+2(N/2)+N=N2/4+N+N<N2
• A 2-point DFT – what is the simplest stage like?
1
P[k ] p[n ]W2nk
0 k 1
p[0]
P[0]=p[0]+p[1]
p[1]
P[1]=p[0]+W21p[1]
n 0
P[0] p[0]W20 p[1]W21 p[0] p[1]
P[1] p[0]W20 p[1]W21 p[0] W21 p[1]
W21 1
18
8.3.5 DIT-FFT: Procedure (Cont’d)
19
8.3.6 DIT-FFT: Complexity
• With this decomposition, each stage requires N complex mults
and ~N complex adds
• There are log2N stages
• Thus Nlog2N complex mults and adds, not N2
• Further reduction results in each stage of Butterflies of form as
xm[p]
xm[q]
In general
r
N
W
WN( r N 2)
Xm+1[p]
Xm+1[q]
Xm+1 [p]=Xm [p]+WrNXm [q]
Xm+1 [q]=Xm [p]+WN(r+N/2)Xm [q]
20
8.3.6 DIT-FFT: Complexity (Cont’d)
But
WNr N 2 WNrWNN 2 WNr ( 1) WNr
WNN 2 e
j
2 N
N 2
1
Then the general Butterfly can be written as
xm[p]
xm[q]
Xm+1[p]
r
N
W
So
-1
Xm+1[q]
X m1[ p] X m [ p] WNr [q]
X m1[q] X m [ p] WNr [q]
21
8.3.6 DIT-FFT: Complexity (Cont’d)
• This signal flow graph saves one complex multiply per
butterfly – which cuts the total by half
# complex mults = N/2 log2N
# complex adds = Nlog2N
• As opposed to N2 of the DFT direct computation.
• Consider reductions:
N
DFT
FFT
Ratio
256
2562=65536
256*8=2048
32:1
512
5122=262144
512*9=4608
57:1
1024
10242=1048576
1024*10=10240
102:1
2048
22
8.3.7 DIT-FFT: In-Place Computations
• Due to “flow of calculation” from left to right, 2-point
calculations in each Butterfly – only a single complex
storage array is needed in the computation
• Bit reward ordering: data is stored in “bit-reversed” order
due to the even/odd splits on the data.
Index
Binary
Bit Rev
Order
0
000
000
0
1
001
100
4
2
010
010
2
3
011
110
6
4
100
001
1
5
101
101
5
6
110
011
3
7
111
111
7
Even
Odd
23
8.3.8 Decimation-In-Frequency (DIF)
• Subdivide X[k] into successively smaller subsequences
• Let N=2v Even
X [2 r ]
N 2 1
( x[n] (1)
n 0
Odd
X [2r 1]
2r
x[n N 2 ])WN2 rn
N 2 1
N ])W nr W n
(
x
[
n
]
x
[
n
2
N 2 N
n 0
• These are 2 N/2 point DFTs
Let
g [n ] x[n ] x[n N 2 ]
and
h[n ] x[n ] x[n N 2 ]
X [2 r ]
N 2 1
(g[n])W
rn
N 2
n 0
X [2 r 1]
N 2 1
n
nr
(
h
[
n
]
)
W
W
N
N 2
n 0
24
8.3.9 DIF-FFT: Signal Flow Graph
25
8.4.1 Inverse FFT: Ideas
1
x[n]
N
Want to compute IDFT using FFT
N 1
X [k ]W
k 0
nk
N
IDFT
1. Multiply IDFT by N
2. Take the complex conjugate
N 1
*
*
( Nx[n ]) Nx [n ] X [k ]WNnk
k 0
*
N 1
X *[k ]WNnk
k 0
3. Take the conjugate again
which is the DFT of
4. Divided by N
X*[k]
1
Nx[n] N
N
N 1
X [k ]W
k 0
nk
N
*
N 1 *
Nx[n ] X [k ]WNnk
k 0
*
1 N 1 *
x[n] X [k ]WNnk
N k 0
26
8.4.2 Inverse FFT: Procedure and Example
1.
2.
3.
4.
Get X*[k];
Compute FFT(X*[k]);
Compute conjugate to get Nx[n];
Divide result by N.
Original Image A
FFT(A)
Sine Wave B
A+B
FFT(A+B)
FFT(A+B)-FFT(B)
FFT(B)
IFFT[FFT(A+B)-FFT(B)] 27
8.5.1 Signal Analysis using DFT
• Anti-aliasing filter is to incorporated to eliminate or minimize the effect
of aliasing when the continuous-time is converted to a sequence.
• The need for multiplication of x[n] by w[n] is a consequence of the
finite-length requirement.
• A finite-duration window w[n] is applied to x[n] prior to computation of
the DFT in order to smooth sharp peaks and discontinuities in X (e j ) .
• The C/D conversion is represented in the frequency-domain as
1
2r
X (e ) X c j j
T r T
T
j
28
Illustration of the Fourier transforms:
(a) Fourier transform of continuous-time input signal.
S ( j)
(b) Frequency response of anti-aliasing filter.
H aa ( j)
(c) Fourier transform of output of anti-aliasing filter.
X c ( j) Sc ( j) H aa ( j)
(d) Fourier transform of sampled signal.
1
2r
X (e ) X c j j
T r T
T
j
(e) Fourier transform of window sequence.
V (e j ) W (e j ) * X (e j )
(f) Fourier transform of windowed segment and
frequency samples obtained using DFT samples.
N 1
V [k ] v[n]e j ( 2 / N ) kn V (e j )
n 0
k
2k
2k
2
, k
, N
(fundamenta l frequency)
NT
N
N
2k / N
29
8.5.2 Signal Analysis using DFT: Example-1
30
8.5.3 Signal Analysis using DFT: Example-2
31