Transcript 9.3 decimation-in-time FFT Algorithms

Biomedical Signal processing
Chapter 9 Computation of the Discrete
Fourier Transform
Zhongguo Liu
Biomedical Engineering
School of Control Science and Engineering, Shandong
University
2015/4/13
1
Zhongguo Liu_Biomedical Engineering_Shandong Univ.
Chapter 9 Computation of the
Discrete Fourier Transform
9.0 Introduction
9.1 Efficient Computation of Discrete Fourier Transform
9.2 The Goertzel Algorithm
9.3 decimation-in-time FFT Algorithms
9.4 decimation-in-frequency FFT Algorithms
9.5 practical considerations（software realization)
2
9.0 Introduction
Implement a convolution of two sequences
by the following procedure:
1. Compute the N-point DFT X 1 k  and X 2 k 
of the two sequence x1 n and x2 n
2. Compute X 3 k   X1k X 2 k  for 0  k  N  1
3. Compute x3 n  x1n N x2 n as the inverse
DFT of X 3 k 
Why not convolve the two sequences directly?
There are efficient algorithms called Fast
Fourier Transform (FFT) that can be orders of
magnitude
more
efficient
than
others.
3
9.1 Efficient Computation of Discrete
Fourier Transform
The DFT pair was given as
N 1
X  k    x[n]e
n 0
 j  2 / N kn
1
x[n] 
N
N 1
 X k  e
k 0
DFTcomputational complexity:
Each DFT coefficient requires
N complex multiplications;
N-1 complex additions
All N DFT coefficients require
N2 complex multiplications;
N(N-1) complex additions
4
j  2 / N kn
4
9.1 Efficient Computation of Discrete
Fourier Transform
N 1
X  k    x[n]e
 j  2 / N kn
n 0
Complexity in terms of real operations
4N2 real multiplications
2N(N-1) real additions (approximate 2N2)
5
5
9.1 Efficient Computation of Discrete
Fourier Transform
Most fast methods are based on Periodicity properties
Periodicity in n and k; Conjugate symmetry
e
 j  2 / N k  N n 
e
Re
6
 j  2 / N kn

e
e
 j  2 / N kN  j  2 / N k   n 
e
 j  2 / N k  n  N 
e
e
j  2 / N kn
j  2 / N  k  N n
]
6
9.2 The Goertzel Algorithm
Makes use of the periodicity e j  2 / N  Nk  e j 2 k  1
Multiply DFT equation with this factor
X k   e
Define
j  2 / N kN
N 1
 x[r ]e
 j  2 / N rk
r 0
yk  n 

 x[r ]e
N 1
  x[r ]e
r 0
j  2 / N k  n r 
r 
using x[n]=0 for n<0 and n>N-1
j  2 / N k  N r 
u n  r 
X  k   yk  n  n  N
X[k] can be viewed as the output of a filter to the input
x[n]
j  2 / N kn
u  n
Impulse response of filter: h[n]  e
X[k] is the output of the filter at time n=N
7
7
9.2 The Goertzel Algorithm
Goertzel h[n]  e j  2 / N knu[n]  W  knu[n]
Filter:
N
1
Hk  z  
1  WN k z 1
yk [n]  yk [n 1]WNk  x[n], n  0,1,..., N ,
X  k   yk  n  n  N , k  0,1,..., N
yk [1]  0
N 1
X  k    x[n]WNkn
n 0
Computational complexity
4N real multiplications; 4N real additions
Slightly less efficient than the direct method
But it avoids computation and storage of WNkn
8
Second Order Goertzel Filter
Goertzel Filter
1
Hk  z 
2
1 e
j
N
k 1
z
Multiply both numerator and denominator
 j 2 k
2
k
N
1 e
z 1
Hk  z  

2
2
2 k 1 2
j
k
 j k 1 




1
N
N
1

2
cos
z

z
1

e
z
1

e
z



N



2 k
y[n]   y[n  2]  2 cos
y[n  1]  x[n], n  0,1,..., N
N
yk [ N ]  y[ N ] WNk y[ N 1]  X k  , k  0,1, ..., N
1 e
9
N
z 1
j
9
Second Order Goertzel Filter
2 k
y[n]   y[n  2]  2 cos
y[n  1]  x[n], n  0,1,..., N
N
yk [ N ]  y[ N ] WNk y[ N 1]  X k  , k  0,1, ..., N
Complexity for one DFT coefficient ( x(n) is complex
sequence).
Poles: 2N real multiplications and 4N real additions
Zeros: Need to be implement only once:
4 real multiplications and 4 real additions
Complexity for all DFT coefficients
Each pole is used for two DFT coefficients
Approximately N2 real multiplications and 2N2 real
additions
10
10
Second Order Goertzel Filter
2 k
y[n]   y[n  2]  2 cos
y[n  1]  x[n], n  0,1,..., N
N
yk [ N ]  y[ N ] WNk y[ N 1]
 X k , k  0,1, ..., N
If do not need to evaluate all N DFT coefficients
Goertzel Algorithm is more efficient than FFT if
less than M DFT coefficients are needed,M < log2N
11
11
9.3 decimation-in-time FFT Algorithms
Makes use of both periodicity and symmetry
Consider special case of N an integer power of 2
Separate x[n] into two sequence of length N/2
Even indexed samples in the first sequence
Odd indexed samples in the other sequence
N 1
X  k    x[n]e
 j  2 / N kn
n 0

 x[n]e
n even
12
 j  2 / N kn

 x[n]e
 j  2 / N kn
n odd
12
9.3 decimation-in-time FFT Algorithms
X k  
 x[n]e
 j  2 / N kn
n even

 x[n]e
 j  2 / N kn
n odd
Substitute variables n=2r for n even and n=2r+1 for odd
X k  

N /21
 x[2r ]W
2 rk
N
r 0
N /21

r 0

N /21
 x[2r  1]W
 G k   W H k 
k
N
N
r 0
x[2r ]WNrk/2  WNk
 2 r 1k
N /21

x[2r  1]WNrk/2
WN2
 j 2 2
e N
r 0
 j 2
 e N /2
G[k] and H[k] are the N/2-point DFT’s of each
13 subsequence
 WN /2
13
9.3 decimation-in-time FFT Algorithms
X k  
N /2 1
 x[2r ]W
rk
N /2
r 0
W
 G k   W H k 
k
N
N 1
k  0,1,...,
2
N

G k    G k 
2

k
N
N /2 1
 x[2r  1]W
rk
N /2
r 0
 j 2 2rk
e N
 j 2 rk
 e N /2
k  0,1,..., N
N

H k    H k 
2

G[k] and H[k] are the N/2-point DFT’s of each
subsequence
14
W
rk
N /2
14
8-point DFT using decimation-in-time
15
Figure 9.3
computational complexity
Two N/2-point DFTs
2(N/2)2 complex multiplications
2(N/2)2 complex additions
Combining the DFT outputs
N complex multiplications
N complex additions
Total complexity
N2/2+N complex multiplications
N2/2+N complex additions
More efficient than direct DFT
16
16
9.3 decimation-in-time FFT Algorithms
Repeat same process , Divide
N/2-point DFTs into
Two N/4-point DFTs
Combine outputs
N=8
17
17
9.3 decimation-in-time FFT Algorithms
After two steps of decimation in
time
Repeat until we’re left with two-point DFT’s
18
18
9.3 decimation-in-time FFT Algorithms
flow graph for 8-point decimation in time
N=2m
m=log2N
Complexity:
19
Nlog2N complex multiplications and additions
19
Butterfly Computation
Flow graph constitutes of butterflies
We can implement each butterfly with one multiplication
Final complexity for decimation-in-time FFT
(N/2)log2N complex multiplications and Nlog2N additions
20
20
9.3 decimation-in-time FFT Algorithms
Final flow graph for 8-point decimation in time
N=2m
m=log2N
Nlog2N
complex
multiplications
and additions
Complexity:
(Nlog2N)/2 complex multiplications and Nlog2N additions
21
21
9.3.1 In-Place Computation同址运算
Decimation-in-time flow graphs require two sets of registers
Input and output for each stage one set of registers
X 0  0  x  0
X 0  1  x  4
X 0  2  x  2
X 0 3  x  6
X 0  4  x 1
X 0  5  x  5
X 0  6  x 3
22X 0
7  x 7
x  0
x  4
x  2
x  6
x 1
x  5
x  3
x 7
X 2  0
X 2  1
X 2  2
X 2 3
X 2  4
X 0
X 1
X  2
X  3
X  4
X 2  5
X 5
X 2  6
X 6
X 2 7
X 7
22
9.3.1 In-Place Computation同址运算
Note the arrangement of the input indices
Bit reversed indexing（码位倒置）
X 0  0  x  0  X 0  000  x  000 x  0
X 0
X 0  1  x  4  X 0  001  x 100 x  4
X 1
X 0  2  x  2  X 0  010  x  010 x  2
X 0 3  x  6  X 0  011  x 110
X 0  4  x 1  X 0 100  x  001
X 0 5  x 5  X 0 101  x 101
X 0  6  x 3  X 0 110  x  011
X 0  7   x  7   X 0 111  x 111
X  2
x  6
X  3
x 1
X  4
x  5
X 5
x  3
X 6
x 7
X 7
23
23
normal order
binary coding for position：
x  000  000
x  0
x n2n1n0 
x  001
001
x  4
x  010 
010
x  2
x  011
011
x  6
x 100 
100
x 1
x 101
101
x 110 
110
x 111
111
must padding 0 to
24
N  2M
x  5
x  3
x 7
Figure 9.13
cause of bit-reversed order
binary coding for position：
x  0
x  000  000
x n2n1n0 
must padding 0 to
25
N  2M
x 100 
001
x  4
x  010 
010
x  2
x 110 
011
x  001
100
x 101
101
x  011
110
x 111
111
x  6
x 1
x  5
x  3
x 7
Figure 9.13
9.3.2 Alternative forms
Note the arrangement of the input indices
Bit reversed indexing（码位倒置）
X 0  0  x  0  X 0  000  x  000 x  0
X 0
X 0  1  x  4  X 0  001  x 100 x  4
X 1
X 0  2  x  2  X 0  010  x  010 x  2
X 0 3  x  6  X 0  011  x 110
X 0  4  x 1  X 0 100  x  001
X 0 5  x 5  X 0 101  x 101
X 0  6  x 3  X 0 110  x  011
X 0  7   x  7   X 0 111  x 111
26
X  2
x  6
X  3
x 1
X  4
x  5
X 5
x  3
X 6
x 7
X 7
26
9.3.2 Alternative forms
strongpoint：in-place computations
shortcoming：non-sequential access of data
27
Figure 9.14
Figure 9.15
28
Input, output same normal order
shortcoming：not in-place computation
non-sequential access of data
Figure 9.16
29
Same structure, Input, output same normal order
shortcoming：not in-place computation
strongpoint: sequential access of data
9.4 Decimation-In-Frequency FFT Algorithm
The DFT equation
N 1
X  k    x[n]WNnk
n 0
Split the DFT equation into even and odd frequency indexes
N 1
N /2 1
X  2r    x[n]WNn 2 r 

n 0
Substitute
variables



30
n 0
N /21
 x[n]W
n 0
N /2 1
n 2r
N
x[n]WNn 2 r 

N 1

n N / 2
N /21
x[n]WNn 2 r
 x[n  N / 2]W
 n  N / 22 r
N
n 0
  x[n]  x[n  N / 2]W
nr
N /2
n 0
N /2 1

n 0
g (n)WNrn/2
30
9.4 Decimation-In-Frequency FFT Algorithm
N 1
X  k    x[n]WNnk
The DFT equation
N 1
X  2r  1   x[n]W
n ( 2 r 1)
N
n 0



N /21
 x[n]W
n 0
N /2 1
n (2 r 1)
N


N /21
 x[n]W
N /21
N 1


n ( 2 r 1)
N
x[n]W
n  N /2
 x[n  N / 2]W
 n  N /2(2 r 1)
N
n 0
  x[n]  x[n  N / 2]W
n (2 r 1)
N
n 0
N /21
  x[n]  x[n  N / 2]W W
n 2r 1
N
31
n ( 2 r 1)
N
n 0
n
N
n 0
W
n 0
W W W W
2rn
N
n
N
rn
N /2
n
N
rn
N /2

N /2 1

n 0
N
(2 r 1)
2
N
W
h(n)WNnWNrn/2
W W
Nr
N
N /2
N
 1
31
decimation-in-frequency decomposition of an Npoint DFT to N/2-point DFT
X  2r  
N /2 1

n 0N /2 1
X  2r  1 
32
nr
x
[
n
]

x
[
n

N
/
2]
W

 N /2

n 0
N /2 1

n 0
rn
N /2
n
x
[
n
]

x
[
n

N
/
2]
W

 NW

g (n)WNrn/2
N /2 1

n 0
h(n)WNnWNrn/2
32
decimation-in-frequency decomposition of an 8point DFT to four 2-point DFT
X  2*2s  
N / 41
sn

[
g
(
n
)

g
(
n

N
/
4)]
W

N /4
n 0
X  2*(2s  1) 
33
N / 41
 [ g (n)  g (n  N / 4)]W W
n 0
2n
N
sn
N /4
N /4 1

n 0

p(n)WNsn/4
N /4 1

n 0
q(n)WN2 n33
WNsn/4
2-point DFT
X v ( p)  X v1( p)  X v1(q)
X v (q)   X v 1 ( p)  X v 1 (q)  W
0
8
34
when N  8
34
X  2r  
N /2 1
  x[n]  x[n  N / 2]
n 0
X  2* 2 s  

N /4 1

n 0


n 0
g (n)WN2/2sn 
g (n)WN2/2sn 
sn
N /4
g (n)W

N /4 1

n 0
N /4 1


N /2 1

n 0
g (n)WNrn/2
N /2 1

n  N /4
g (n)WN2/2sn
g ( n  N / 4)WN2/2s ( n  N /4)
g ( n  N / 4)W
sn
N /4
n 0
N /4 1
 [ g (n)  g (n  N / 4)]W
n 0
35

N /4 1
n 0

N /4 1
WNnr/2
sn
N /4

N /4 1

n 0
p(n)WNsn/4
X  2r  
N /2 1
  x[n]  x[n  N / 2]
n 0
X  2*(2 s  1)  

N /4 1

n 0

N /4 1

n 0


n 0
g (n)WN(2/2s 1) n 

sn
N /4
2n
N
g (n)W W
N /2 1
36

n 0
g (n)WNrn/2

n N / 4

g (n)WN(2/2s 1) n
N /4 1

n 0
N /4 1

g (n  N / 4)WN(2/2s 1)( n  N /4)
g (n  N / 4)W W W
sn
N /4
n 0
N /4 1
 [ g (n)  g (n  N / 4)]W
n 0

g ( n)WN(2/2s 1) n
g (n)WNsn/4WNn/2 
N /4 1
n 0

N /2 1
nr
WN /2
N /2 1
2n
N
sn
N /4
W

2n
N
(2 s 1) N /4
N /2
N /4 1

n 0
q(n)WN2 nWNsn/4
WN(2/2s1) N /4  WNsN/2/2WNN/2/4  1
X  2*2s  
N /4 1

n 0
X  2* 2* 2t  

2 tn
N /4
p(n)W
n 0

N /4 1
N /81

N /81

2tn
N /4
p(n)W
n 0

p(n)WNsn/4

2 tn
N /4
p(n)W
n 0

N /4 1

2 tn
N /4
p(n)W
n  N /8

N /81

2t ( n  N /8)
N /4
p(n  N / 8)W
n 0
N /81
tn
[
p
(
n
)

p
(
n

N
/
8)]
W

N /8
n 0
 p(n)  p(n  1)
37
when N  8
X  2*2s  
N /4 1

n 0
p(n)WNsn/4
X  2*2*(2t  1)  

N /81

n 0


p(n)WN(2/4t 1) n 
n 0


p(n)WN(2/4t 1) n 
n 0
N /81

N /4 1
N /81

n 0
N /81
p(n)WN(2/4t 1) n
N /4 1

n N /8
N /81

n 0
tn
p(n)WN2/4
WNn/4 
p(n)WN(2/ 4t 1) n
p(n  N / 8)WN(2/4t 1)( n N /8)
N /81

n 0
tn
p(n  N / 8)WN2/4
WNn/4WN(2/4t 1) N /8
tn
4n
[
p
(
n
)

p
(
n

N
/
8)]
W
W

N /8 N
n 0
 [ p(n)  p(n  1)]W80 when N  8
38
WN(2/4t 1) N /8  WNtN/4/4WNN/4/8  1
Final flow graph for 8-point DFT decimation in
frequency
N=8
39
39
9.4.1 In-Place Computation同址运算
DIF FFT
DIT FFT
40
40
9.4.1 In-Place Computation同址运算
DIF FFT
transpose
DIT FFT
41
41
9.4.2 Alternative forms
decimation-in-frequecy Butterfly Computation
transpose
decimation-in-time Butterfly Computation
42
42
The DIF FFT is the transpose of the DIT FFT
DIF FFT
transpose
DIT FFT
43
43
9.4.2 Alternative forms
DIF FFT
transpose
DIT FFT
44
9.4.2 Alternative forms
DIF FFT
transpose
DIT FFT
45
Figure 9.24 erratum
x  0
x  4
x  2
x  6
x 1
x  5
x  3
x 7
46
9.4.2 Alternative forms
DIF FFT
transpose
DIT FFT
47
Chapter 8 HW
9.1, 9.2, 9.3,
48
返
回2015/4/13
上一页
Zhongguo Liu_Biomedical Engineering_Shandong Univ.
下一页