Lecture 9 - Discrete Fourier Transform and Fast Fourier Transform (I) James Barnes ([email protected]) Spring 2009

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 1 / 16

Introduction ●

The next two lectures cover the Discrete Fourier Transform (DFT) and the Fast Fourier Transform technique for speeding up computation by reducing the number of multiplies and adds required.

●

To compute the DFT, we sample the Discrete Time Fourier Transform in the frequency domain, specifically at points spaced uniformly around the unit circle. The resulting DFT when transformed back to the time domain, contains images of the original time-domain sequence.

●

The DFT uses a factor WN which is the N-th root of unity. The FFT algorithm exploits symmetry in WN to greatly reduce the number of operations.

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 2 / 16

Continuous Time Fourier Pairs In continuous time, we represent the frequency spectra of periodic and aperiodic signals with Fourier Series and Fourier Transforms, respectively ●

Fourier Series expansion of continuous-time periodic signal with fundamental frequency Fo , period To . x(t) =

∞ X

j2πkFo t

ck e

,

k=−∞

1 ck = T

toZ+To

x(t)e−j2πkFo t dt.

(1)

to

where the {ck } give the spectral characteristics. ●

Fourier Transform for continuous-time aperiodic signal x(t) =

Z∞

X(F )ej2πF t dF,

−∞

Colorado State University Dept of Electrical and Computer Engineering

X(F ) =

Z∞

x(t)e−j2πF t dt.

(2)

−∞

ECE423 – 3 / 16

Discrete-time Fourier Pairs In discrete time, we replace t with nto in the previous, where to is the sampling period. ●

Fourier Series expansion of discrete-time periodic with period N sequence times, frequency 1/N: x[n] =

N −1 X

ck ej

k=0

●

2π N kn

,

N −1 1 X −j 2π ck = x[n]e N kn N n=0

k = 0, 1, . . . , N − 1.

(3)

Fourier Transform of discrete-time aperiodic signal (DTFT): 1 x[n] = 2π

Zπ

X(ω)ejωn dω,

−π

Colorado State University Dept of Electrical and Computer Engineering

X(ω) =

∞ X

x[n]e−jωn ,

ω = [−π, π]

(4)

n=−∞

ECE423 – 4 / 16

DTFT to DFT We can also write the DTFT as X(ω) = X(z)|z=ejω ,

(5)

where X(z) is the z-transform. But X(ω) is still a continuous function of ω. The only way we can use it in computation is in an analog computer (!), which is not practical. =⇒ we need a version of X(ω) which is defined for discrete frequency values. DFT: evaluate (sample) Z transform at N discrete points uniformly spaced around the unit circle 2πk , k = 0, 1, . . . , N − 1. (6) N Note that we define ωk to lie in the range [0, 2π] instead of [−π, π], but this is just re-mapping the negative frequency part to [π, 2π]. X[k] = X(z)|z=ejωk ,

Colorado State University Dept of Electrical and Computer Engineering

ωk =

ECE423 – 5 / 16

Artifact from Sampling Frequency How well does the sampled FT X[k] represent the original sequence in the frequency domain? In other words, can we go back X[k] → x[n] and recover x[n]?

X(ωk ) ≡ X[k] =

∞ X

x[n]e−j

2π N kn

= . . .+

n=−∞

−1 X

x[n]e−j

2π N kn

N −1 X

+

2π

x[n]e−j N kn +. . .

n=0

n=−N

(7) where we have split up the infinite sum into a number of segments encompassing N values of the sequence x[n]. We can rewrite the above: X[k] =

∞ X

m=−∞

X

(m+1)N −1 2π

x[n]e−j N kn .

(8)

n=mN

Now let n0 = n − mN and change order of summation. Then X[k] =

N −1 X

{

∞ X

0

0 −j 2π N k(n +mN )

x[n + mN ]e

}.

(9)

n0 =0 m=−∞ Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 6 / 16

Artifact from Sampling Frequency (2) 0 −j 2π N k(n +mN )

But in (9), e

X[k] =

N −1 X

0 −j 2π N kn

=e

xp [n0 ]e−j

n0 =0

. So we can write (9) as

0 2π N kn

,

xp [n] ≡

∞ X

x[n0 + mN ].

(10)

m=−∞

This is the ”inverse” Fourier Series expansion of a periodic signal xp [n0 ] , which is the ”periodically extended” version of x[n], with a fundamental period of N, the length of the original sequence x[n]. xp [n] thus contains an infinite number of replicas (”images”) of x[n] displaced by multiples of N. We will see how to determine xp [n] later. For now, we note that X[k] is related to the coefficients Fourier Series expansion of xp [n]. The images in xp would cause aliasing distortion if they overlapped x[n]. To avoid this distorition, the order N of the DFT (the number of sample points around the unit circle) must be greater than or equal to the length of x[n].

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 7 / 16

Artifact from Sampling Frequency (3)

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 8 / 16

Computing IDFT (recovered xp [n]) xp [n] is periodic in n with period N, so it can be expressed as a Fourier Series xp [n] =

N −1 X

2π

ck ej N kn ,

n = 0, . . . , N − 1.

(11)

k=0

Then the coefficients ck are given by: N −1 2π 1 X xp [n]e−j N kn , ck = N n=0

k = 0, . . . , N − 1.

(12)

Comparing the last with the left hand side of (10), we see 1 ck = X[k], N

N −1 2π 1 X xp [n] = X[k]ej N kn N

(13)

k=0

Then with NO TIME ALIASING, x[n] = xp [n] , 0 ≤ n ≤ N − 1.

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 9 / 16

Remarks on Sequence Length and Accuracy of Spectral Analysis What if x[n] is a segment of a sequence of arbitrary length? We need to window x[n]. Note: clearly, the frequency characteristics will only represent those of the signal during the window period. ●

The windows used are similar to those we saw with FIR: rectangular, Hanning,...

●

Windowing introduces an artifact called ”leakage” due to the finite length of the windowed sequence.

Consider rectangular window w[n] of length L. In the time domain, we compute the (inner) product x[n] w[n]. We know this is equivalent to convolution in the frequency domain: x[n] w[n] ↔ X(ω) ? W (ω).

(14)

The FTrans of a rectangular window is the sinc(.) function, convolved with the original spectrum with the sync(.) function. This causes the spectrum to smear out and ”leak”. This effect is the basis for the Heisenberg Uncertainty Principle of Physics. Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 10 / 16

FFT Notation Notation X[k] =

N −1 X

−j

x[n] e

2π N kn

n=0

=

N −1 X

x[n] WNkn ,

−j 2π N

WN = e

.

(15)

n=0

WN is the N-th root of unity; called ”twiddle” factor. This definition agrees with standard textbooks, but your lab handout defines WN without the minus sign. This difference can be compensated for easily. Inverse DFT can be written: N −1 N −1 X 2π 1 1 X X[k]E j N kn = X[k] WN−kn x[n] = N N k=0

Colorado State University Dept of Electrical and Computer Engineering

(16)

k=0

ECE423 – 11 / 16

The FFT Idea The DFT as given in eqn(15) requires 4N 2 real multiplications and N (4N − 2) real additions, so the comuputation is Order(N 2 ). There are two ways of reducing this computational load: 1. Break the problem down into smaller pieces which can be re-combined. Sort of like the ”quicksort” approach to sorting; 2. Exploit symmetries in the twiddle factors WN

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 12 / 16

Useful Symmetry Properties of Twiddle Factors 1.

m+ N WN 2

= −WNm

(17)

(When you go half way around the unit circle, you change sign). 2. WNm+N = WNm

(18)

(Adding N is equivalent to going one full revolution around the unit circle). 3. WNm = W N

m

(19)

(Easily shown from the definition of WN ).

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 13 / 16

FFT Approaches 1. Decimation in time 2. Decimation in frequency - this course will not cover Decimation in time (DIT) approach: ●

break the sequence down by repeatedly halving in length – all the way to sequence length 2, and compute FFTs of the shorter chunks.

●

DIT requires the sequence length N to be power of 2: N = 2ν , ν is a integer

●

The smaller sequence FFTs are then recombined in ν stages to produce the FFT for the whole sequence

●

Exploits symmetry in the twiddle factors.

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 14 / 16

Example: 8 pt FFT

X[k] =

N −1 X

x[n]WNkn ,

W = e−j

2∗π N kn

,

k, n = 0, 1, . . . N − 1,

N = 8.

(20)

n=0

Split the summation into two sums over even and odd values of n, as shown below: X

(N/2)−1

X[k] =

X

(N/2)−1 k·(2r) x[2r]WN

r=0

+

k·(2r+1)

x[2r + 1]WN

(21)

r=0

But using WN2 = WN/2 , easily shown from the definition of WN , we can rewrite the last expression as: X

X[k] =

G[k]+WNk H[k],

G[k] =

r=0

X

(N/2)−1

(N/2)−1 k·r x[2r]WN/2 ,

H[k] =

kr x[2r+1]WN/2

r=0

(22) Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 15 / 16

Example: 8 pt FFT - signal flow diagram of first stage of Decimation in Time

Colorado State University Dept of Electrical and Computer Engineering

ECE423 – 16 / 16