Transcript The Fast Fourier Transform

The Fast Fourier Transform
(and DCT too…)
Nimrod Peleg
Oct. 2002
1
Outline
• Introduce Fourier series and transforms
• Introduce Discrete Time Fourier Transforms,
(DTFT)
• Introduce Discrete Fourier Transforms (DFT)
• Consider operational complexity of DFT
• Deduce a radix-2 FFT algorithm
• Consider some implementation issues of FFTs
with DSPs
• Introduce the sliding FFT (SFFT) algorithm
2
The Frequency Domain
• The frequency domain does not carry any
information that is not in the time domain.
• The power in the frequency domain is that it is simply
another way of looking at signal information.
• Any operation or inspection done in one domain is
equally applicable to the other domain, except that
usually one domain makes a particular operation or
inspection much easier than in the other domain.
• frequency domain information is extremely important
and useful in signal processing.
3
3 Basic Representations for FT
• 1. An Exponential Form
• 2. A Combined Trigonometric Form
• A Simple Trigonometric Form
4
The Fourier Series: Exponential Form
Periodic signal expressed as infinite
sum of sinusoids.
Complex Numbers !
xp (t ) 
ck 
1
Tp

c e
k  

k
jk  0 t
,
w here
x p ( t )e  jk  0 t dt
Tp
 Ck’s are frequency domain amplitude and phase representation
 For the given value xp(t) (a square value), the sum of the first four
terms of trigonometric Fourier series are: xp(t)  1.0 + sin(t) + C2
sin(3t) + C3sin(5t)
5
The Combined Trigonometric Form
• Periodic signal: xp (t) = xp(t+T) for all t
and cycle time (period) is: T  1  2
f0
0
f0 is the fundamental frequency in Hz
w0 is the fundamental frequency in radians:   2 f
xp(t) can be expressed as an infinite sum of orthogonal functions.
When these functions are the cosine and sine, the sum is called the
Fourier Series. The frequency of each of the sinusoidal functions in
the Fourier series is an integer multiple of the fundamental frequency.
Basic frequency + Harmonies
6
Fourier Series Coefficients
jk t
C
k
e
• Each individual term of the series,
,is the
frequency domain representation and is generally
complex (frequency and phase), but the sum is
real.
• The second common form is the combined

trigonometric form:
0
xp ( t )  C 0   2 Ck sin( k 0t  k )
k 1
k  tan 1
Im(Ck )
Re(Ck )
Again: Ck are Complex Numbers !
7
The Trigonometric Form
Reminder:

xp (t )=A0 +  AkCos(k 0t )+BkSin(k 0t )
e j  e-j
Cos =
2
e j -e-j
Sin =
2j
k =1
A0 
1
xp (t )dt = DC = Average value of xp (t )

T p Tp
Ak =  xp (t )Cos(k 0t )dt
Tp
Bk =  xp (t )Sin(k 0t )dt
Tp
All three forms are identical and are related using Euler’s identity:
e  j  Cos  jSin
Thus, the coefficients of the different forms are related by:
2Ck  Ak  jBk
k  tan 1
;
C 0 = A0
Bk
Im(Ck )
 tan 1
Ak
Re(Ck )
8
The Fourier Transform
1/3
The Fourier series is only valid for periodic signals.
For non-periodic signals, the Fourier transform is used.
Most natural signals are not periodic (speech).
We treat it as a periodic waveform with an infinite period.
If we assume that TP tends towards infinity, then we can produce
equations (“model”) for non-periodic signals.
If Tp tends towards infinity, then w0 tends towards 0.
Because of this, we can replace w0 with dw, and it leads us to:
1  fp  
Tp
2
1 d

t  Tp
2
lim
9
The Fourier Transform
2/3
d 
C ( ) 
x ( t )e  jt dt
2 
z
 Increase TP = Period Increases : No Repetition:
 Discrete frequency variable becomes continuous:
 Discrete coefficients Ck become continuous:
1 
d


Tp 2
2
k 0  
Ck  C( )

C ( )
2
  x(t )e jt dt
d 
10
The Fourier Transform
We define: 2

C ( )
d
3/3
X ( )  F{ x(t )

1
X ( )  x ( t )e  jt dt  x(t ) 
x ( )e jt d
2 

z
z
11
Signal Representation by Delta Function
Instead of a continuous signal we have a “collection of samples
This is equivalent to sampling the signal with one Delta
Function each time, moving it along X-axis, and summing
all the results:

xs( t )   x ( t ) ( t  nTs)

Note that the Delta is “1” only
If its index is zero !
12
Discrete Time Fourier Transform
1/3
• Consider a sampled version, xs(t) , of a continuous

signal, x(t) :
xs( t )   x ( t ) ( t  nTs)

Ts is the sample period. We wish to take the Fourier transform of this
sampled signal. Using the definition of Fourier transform of xs(t)
and some mathematical properties of it we get:

 Replace continuous time t with (nTs)
 Continuous x(t) becomes discrete x(n)
 Sum rather than integrate all discrete samples
xs( )   x ( nTs) e  jnTs
n 
13
Discrete Time Fourier Transform
Fourier
Transform x( ) 


x(t )e
 j t
dt

 x(n)e
x() 
Discrete Time
Fourier Transform
 jn
n 

Inverse

1
j
x
(

)
e
d
Fourier x(t ) 

2 
Transform
t

x(n) 
1
j( n )
x
(

)
e
2 
2/3
Inverse
Discrete
d
Time Fourier
Transform
 Limits of integration need not go beyond ± because the spectrum
repeats itself outside ± (every 2): X ( )  X (   2 )
 Keep integration because X ( ) is continuous:   Ts
means that  is periodic every Ts
14
Discrete Time Fourier Transform
3/3
• Now we have a transform from the time
domain to the frequency domain that is
discrete, but ...
DTFT is not applicable to DSP because it
requires an infinite number of samples and
the frequency domain representation is a
continuous function – impossible to represent
exactly in digital hardware.
15
st
1
result: Nyquist Sampling Rate
1/2
• The Spectrum of a sampled signal is
periodic, with 2*Pi Period: X ()  X (  2 )
Easy to see: X (  2 ) 


x(n)e jn (  2 ) 
n 





x(n)e jn e j 2 n
n 
x(n)e jn  X ()
n 
e j 2 n  cos(2 n)  j sin(2 n)  1
16
1st result: Nyquist Sampling Rate
2/2
• For maximum frequency wH :
|   
  Ts
H
   HTs

Ts 
H
 s 

BUT : Ts 
s
2
=2 H
Ts
17
Practical DTFT
 Take only N time domain samples
x() 

 x ( n )e
 j n
n  
N 1
x()   x(n)e
 j n
n 0
 Sample the frequency domain, i.e. only evaluate x() at N
discrete points. The equal spacing between points is  = 2/N
N 1
2k
 j 2 kn / N
x(
)   x ( n) e
N
n0
k  0,1,2,..., N  1
18
The DFT
Since the only variable in 2k / N is k , the DTFT is written:
N 1
x ( k )   x ( n) e
 j 2 kn / N
k  0,1,2,..., N  1
n0
 j 2 / N
W
N

e
Using the shorthand notation:
(Twiddle Factor)
The result is called Discrete Fourier Transform (DFT):
N 1
X N (k )  
n0
x(n)WNkn
and
1 N 1
 kn
x ( n) 
X
(
k
)
W
 N
N
N k 0
19
Usage of DFT
• The DFT pair allows us to move between
the time and frequency domains while using
the DSP.
• The time domain sequence x[n] is discrete
and has spacing Ts, while the frequency
domain sequence X[k] is discrete and has
spacing 1/NT [Hz].
20
DFT Relationships
Time Domain
Frequency Domain
|x(k)|
N Samples
N Samples


0
Ts
0
1
2Ts 3Ts (N-1)Ts
2
3
N-1

X(n)
t
0
1
2
N/2
n
0
Fs
N
2Fs
N
Fs
2
N-2 N-1
k
 2Fs
N
f
 Fs
N
21
Practical Considerations
N 1
 Standard DFT:
X N (k )   xn (k )WNkn
0  k  N 1
n 0
 An example of an 8 point DFT:
7
X N (k )   xN (k )W7kn
k  0,1, 2,..., 7
n 0
 Writing this out for each value of n :
Xn (k)  x(0)W7k 0  x(1)W7k1  .......  x(7)W7k7,k  0,1,...,7
k0
 Each term such as x(0)W7
requires 8 multiplications
 Total number of (complex !) multiplications required: 8 * 8 = 64
 1000-point DFT requres 10002 = 106 complex multiplications
And all of these need to be summed….
22
Fast Fourier Transform
Symmetry Property
WNk N / 2  WNk
Periodicity Property
WNk N  WNk
N
1
2
N
1
2
r 0
r 0
XN (k )   x(2r ).WN2rk  x(2r  1).WN( 2r 1)k
Splitting the DFT in two
(odd and even)
or
N
1
2
N
1
2
r 0
r 0
XN (k )   x(2r ).(WN2 )rk  WNk  x(2r  1).(WN2 )rk
Manipulating the twiddle factor
W e
2
N
j( 
2
2)
N
e
j( 
2
)
N/ 2
 WN / 2
THE FAST FOURIER TRANSFORM
N
1
2
N
1
2
r 0
r 0
Xn (k )   x(2r )WNrk 2  WNk  x(2r  1)WNrk 2
23
FFT complexity
xN (k ) 
N 1
2
N 1
2
r 0
r 0
rk
k
x
(
2
r
)
W

W

N /2
N
(N/2)2 multiplications
rk
x
(
2
r

1
)
W

N /2
(N/2)2 multiplications
N/2 Multiplications
 For an 8-point FFT, 42 + 42 + 4 = 36 multiplications, saving 64 - 36 = 28
multiplications
 For 1000 point FFT, 5002 + 5002 + 500 = 50,500 multiplications, saving
1,000,000 - 50,500 = 945,000 multiplications
24
Time Decimation
 Splitting the original series into two is called decimation in time
 Let us take a short series where N = 8
 Decimate once Called Radix-2 since we divided by 2
n = {0, 1, 2, 3, 4, 5, 6, 7}n = { 0, 2, 4, 6 } and { 1, 3, 5, 7 }
 Decimate again
n = { 0, 4 } { 2, 6 } { 1, 5 } and { 3, 7 }
 The result is a savings of N2 – (N/2)log2N multiplications:
1024 point DFT = 1,048,576 multiplications
1024 point FFT = 5120 multiplication
 Decimation simplifies mathematics but there are more twiddle
factors to calculate, and a practical FFT incorporates these extra factors
into the algorithm
25
Simple example: 4-Point FFT
 Let us consider an example where N=4:
3
X 4 ( k )   x (n)W4kn
0
 Decimate in time into 2 series:
n = {0 , 2} and {1, 3}
1
1
X 4( k )   x (2r )W  W  x (2r  1)W2rk
rk
2
r 0
k
4
r 0
 {x (0)  x (2)W2k }  W4k {x (1)  x (3)W2k }
 We have two twiddle factors.
Can we relate them?
WNk  e
W2k  e
j
2
k
N
j
2
k
2
e
j
2
2k
4
 W22 k
 Now our FFT becomes:
X 4 ( k )  {x (0)  x (2)W42 k }  W4k {x (1)  x (3)W42 k }
26
4-Point FFT Flow Diagram
X 4 ( k )  {x (0)  x (2)W42 k }  W4k {x (1)  x (3)W42 k }
The 2 DFT’s:
for k=0,1,2,3
X 4 ( 0)  {x (0)  x (2)W40 }  W40 {x (1)  x (3)W40 }
For k=0 only:
A ‘flow-diagram’ of it:
x(0)
x(2)
W4
0
x(1)
+
W40
+
W40
x(3)
+
This is for only 1/4 of the
whole diagram !
27
A Complete Diagram
X 4( 0)  {x(0)  x(2)W40}  W40 {x(1)  x(3)W40}
X 4(1)  {x(0)  x(2)W42 }  W41{x(1)  x(3)W42 }
X 4( 2)  {x(0)  x(2)W40}  W42 {x(1)  x(3)W40}
X 4( 3)  {x(0)  x(2)W42 }  W43 {x(1)  x(3)W42 }
X4(0)
x(0)
0
x(2)
1
2
X4(1)
X4(2)
0
x(3)
2
2
k
N
0
2
x(1)
WNk  e
Note:
j
3
X4(3)
W44  e
j
2
4
4
j
2
6
4
W46  e
 1  W40
 1  W42
28
The Butterfly
Twiddle Conversions
A Typical Butterfly
x1
X1
0
W4 = 1
k
X1 = x1 + WN x2
1
W4 = -j
2
x2
k
WN
X2
W4 = -1
k
X2 = x1 – WN x2
3
W4 = j
4 Point FFT Equations
0
X0 = (x0 + x2) + W4 (x1+x3)
4 Point FFT Butterfly
x0
X0
1
X1 = (x0 – x2) + W4 (x1–x3)
X2 = (x0 + x2) –
0
W4
(x1+x3)
X3 = (x0 – x2) –
1
W4
(x1–x3)
x2
X1
0
W4
x1
X2
1
x3
W4
X3
29
Summary
 Frequency domain information for a signal is important for processing
 Sinusoids can be represented by phasors
 Fourier series can be used to represent any periodic signal
 Fourier transforms are used to transform signals
 From time to frequency domain
 From frequency to time domain
 DFT allows transform operations on sampled signals
 DFT computations can be sped up by splitting the original series into
two or more series
 FFT offers considerable savings in computation time
 DSPs can implement FFT efficiently
30
Bit-Reversal
• If we look at the inputs to the butterfly FFT, we can see
that the inputs are not in the same order as the output.
• To perform an FFT quickly, we need a method of
shuffling these input data addresses around to the
correct order.
• This can be done either by reversing the order of the
bits that make up the address of the data, or by pointer
manipulation (bit reversed addition).
• Many DSPs have special addressing modes that allow
them to automatically shuffle the data in the
background.
31
Bit-Reversal example
X4(0)
x(0)
0
x(2)
0
2
1
x(1)
X4(1)
2
X4(2)
3
X4(3)
0
x(3)
2
• To obtain the output
in ascending order
the input values must
be loaded in the
order: {0,2,1,3}
• for 512 or 1024 it is
much more
complicated...
32
8-point Bit-Reversal
• Consider a 3-bit address (8 possible locations).
• After starting at zero, we add half of the FFT
length at each address access with carrying from
left to right (!)
Start at 0 = 000
000+100 = 100
100+100 = 010
010+100 = 110
110+100 = 001
001+100 = 101
101+100 = 011
011+100 = 111
=x(0)
=x(4)
=x(2)
=x(6)
=x(1)
=x(5)
=x(3)
=x(7)
Note that reversing the order
of the address bits gives same
result !
33
And what about DCT ???
  ( n  12 )k 
b(k)x ( n ) cos 


N

n 0
N 1
  ( n  12 )k 
II
b(k)X (k) cos 


N

k 0
N 1
X II (k)=
2
N
DCT Type II*:
x(n) 
2
N
 1 2 if k=0
b(k)= 
 1 if k=1,...,L-1
The rest of the math is quite similar…..
Note: at least 4 types of DCT !!!
* after: A course in Digital Signal Processing, Boaz Porat34
DCT Type II
Type I
Type II
Most used for compression:
JPEG, MPEG etc.
Type III
DCT Basis Vectors for N=8
Type IV
35
DCT Features
•
•
•
•
•
Real transformation
Reversible transformation
2D Transformation exists and separable
Better than the DFT as a “de-correlator”
Fast algorithm exists (NlogN Complexity)
36