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.
WNkN  e j ( 2 / N ) Nk  e j 2k  1
kN
N
X [k ]  W
N 1
 x[l ]W
r 0
kr
N
N 1
  x[l ]WNk ( 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] nN
• It can be interpreted as a discrete convolution of the finite
duration sequence x[n],0  n  N  1 with the sequence WNknu[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
WNk
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  WNk 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(2k / 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 21
 x[2n]W
k ( 2n )
N
n 0

N 21
 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 m1[ p]  X m [ p]  WNr [q]
X m1[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 ]WNnk 
 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 
2r 
 
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 
2r 
 
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 
2k
2k
2
, k 
, N 
(fundamenta l frequency)
NT
N
N
 2k / N
29
8.5.2 Signal Analysis using DFT: Example-1
30
8.5.3 Signal Analysis using DFT: Example-2
31