05 DFT

Download as pdf or txt
Download as pdf or txt
You are on page 1of 20

The Discrete Fourier Transform

The discrete-time Fourier transform (DTFT) of a sequence is a continuous


function of ω, and repeats with period 2π. In practice we usually want to
obtain the Fourier components using digital computation, and can only
evaluate them for a discrete set of frequencies. The discrete Fourier transform
(DFT) provides a means for achieving this.
The DFT is itself a sequence, and it corresponds roughly to samples, equally
spaced in frequency, of the Fourier transform of the signal. The discrete
Fourier transform of a length N signal x[n], n = 0, 1, . . . , N − 1 is given by
N
X −1
X[k] = x[n]e−j(2π/N )kn .
n=0

This is the analysis equation. The corresponding synthesis equation is


N −1
1 X
x[n] = X[k]ej(2π/N )kn .
N
k=0

When dealing with the DFT, it is common to define the complex quantity

WN = e−j(2π/N ) .

With this notation the DFT analysis-synthesis pair becomes


N
X −1
X[k] = x[n]WNkn
n=0
N −1
1 X
x[n] = X[k]WN−kn .
N
k=0

An important property of the DFT is that it is cyclic, with period N , both in the

1
discrete-time and discrete-frequency domains. For example, for any integer r,
N −1 N −1
(k+rN )n
X X
X[k + rN ] = x[n]WN = x[n]WNkn (WNN )rn
n=0 n=0
N
X −1
= x[n]WNkn = X[k],
n=0

since WNN = e−j(2π/N )N = e−j2π = 1. Similarly, it is easy to show that


x[n + rN ] = x[n], implying periodicity of the synthesis equation. This is
important — even though the DFT only depends on samples in the interval 0 to
N − 1, it is implicitly assumed that the signals repeat with period N in both
the time and frequency domains.
To this end, it is sometimes useful to define the periodic extension of the signal
x[n] to be
x̃[n] = x[n mod N ] = x[((n))N ].
Here n mod N and ((n))N are taken to mean n modulo N , which has the
value of the remainder after n is divided by N . Alternatively, if n is written in
the form n = kN + l for 0 ≤ l < N , then

n mod N = ((n))N = l.

x[n]

n
0 N
x̃[n]

n
0 N
Similarly, the periodic extension of X[k] is defined to be

X̃[k] = X[k mod N ] = X[((k))N ].

2
It is sometimes better to reason in terms of these periodic extensions when
dealing with the DFT. Specifically, if X[k] is the DFT of x[n], then the inverse
DFT of X[k] is x̃[n]. The signals x[n] and x̃[n] are identical over the interval 0
to N − 1, but may differ outside of this range. Similar statements can be made
regarding the transform X[k].

1 Properties of the DFT


Many of the properties of the DFT are analogous to those of the discrete-time
Fourier transform, with the notable exception that all shifts involved must be
considered to be circular, or modulo N .
D D D
Defining the DFT pairs x[n]←→X[k], x1 [n]←→X1 [k], and x2 [n]←→X[k],
the following are properties of the DFT:
• Symmetry:

X[k] = X ∗ [((−k))N ]
Re{X[k]} = Re{X[((−k))N ]}
Im{X[k]} = −Im{X[((−k))N ]}
|X[k]| = |X[((−k))N ]|
∢X[k] = −∢X[((−k))N ]

D
• Linearity: ax1 [n] + bx2 [n]←→aX1 [k] + bX2 [k].
D
• Circular time shift: x[((n − m))N ]←→WNkm X[k].
• Circular convolution:
N −1
D
X
x1 [m]x2 [((n − m))N ]←→X1 [k]X2 [k].
m=0

Circular convolution between two N-point signals is sometimes denoted


by x1 [n] N x[n].

3
• Modulation:
N −1
D 1 X
x1 [n]x2 [n]←→ X1 [l]X2 [((k − l))N ].
N
l=0

Some of these properties, such as linearity, are easy to prove. The properties
involving time shifts can be quite confusing notationally, but are otherwise
quite simple. For example, consider the 4-point DFT
3
X
X[k] = x[n]W4kn
n=0

of the length 4 signal x[n]. This can be written as

X[k] = x[0]W40k + x[1]W41k + x[2]W42k + x[3]W43k

The product W41k X[k] can therefore be written as

W41k X[k] = x[0]W41k + x[1]W42k + x[2]W43k + x[3]W44k


= x[3]W40k + x[0]W41k + x[1]W42k + x[2]W43k

since W44k = W40k . This can be seen to be the DFT of the sequence
x[3], x[0], x[1], x[2], which is precisely the sequence x[n] circularly shifted to
the right by one sample. This proves the time-shift property for a shift of
length 1. In general, multiplying the DFT of a sequence by WNkm results in an
N-point circular shift of the sequence by m samples. The convolution
properties can be similarly demonstrated.
It is useful to note that the circularly shifted signal x[((n − m))N ] is the same
as the linearly shifted signal x̃[n − m], where x̃[n] is the N-point periodic
extension of x[n].

4
x[n]

n
0 N
x̃[n]

n
0 N
x̃[n − m]

n
0 N
x[((n − m))N ]

n
0 N

On the interval 0 to N − 1, the circular convolution


N
X −1
x3 [n] = x1 [n] N x2 [n] = x1 [m]x2 [((n − m))N ]
m=0

can therefore be calculated using the linear convolution product


N
X −1
x3 [n] = x1 [m]x̃2 [n − m].
m=0

Circular convolution is really just periodic convolution.


Example: Circular convolution with a delayed impulse sequence
Given the sequences

x1 [n] x2 [n]
n n
0 N 0 N

5
the circular convolution x3 [n] = x1 [n] N x2 [n] is the signal x̃[n] delayed by
two samples, evaluated over the range 0 to N − 1:

x3 [n]
n
0 N

Example: Circular convolution of two rectangular pulses


Let 
1 0≤n≤L−1
x1 [n] = x2 [n] =
0 otherwise.
If N = L, then the N-point DFTs are

N −1 N
X k=0
X1 [k] = X2 [k] = WNkn =
0 otherwise.
n=0

Since the product is



N 2 k=0
X3 [k] = X1 [k]X2 [k] =
0 otherwise,

it follows that the N-point circular convolution of x1 [n] and x2 [n] is

x3 [n] = x1 [n] N x2 [n] = N, 0 ≤ n ≤ N − 1.

Suppose now that x1 [n] and x2 [n] are considered to be length 2L sequences by
augmenting with zeros. The N = 2L-point circular convolution is then seen to
be the same as the linear convolution of the finite-duration sequences x1 [n] and
x2 [n]:

6
x1 [n] = x2 [n]

n
0 L N

L
x1 [n] N x2 [n]

n
0 L N

2 Linear convolution using the DFT


Using the DFT we can compute the circular convolution as follows
• Compute the N -point DFTs X1 [k] and X2 [k] of the two sequences x1 [n]
and x2 [n].
• Compute the product X3 [k] = X1 [k]X2 [k] for 0 ≤ k ≤ N − 1.
• Compute the sequence x3 [n] = x1 [n] N x2 [n] as the inverse DFT of
X3 [k].
This is computationally useful due to efficient algorithms for calculating the
DFT. The question that now arises is this: how do we get the linear convolution
(required in speech, radar, sonar, image processing) from this procedure?

2.1 Linear convolution of two finite-length sequences

Consider a sequence x1 [n] with length L points, and x2 [n] with length P
points. The linear convolution of the sequences,

X
x3 [n] = x1 [m]x2 [n − m],
m=−∞

is nonzero over a maximum length of L + P − 1 points:

7
2

x [n]
1
0
0 L

1
x [n]
2

0
0 P
8
x [n]
3

0
0 L+P−1
n

Therefore L + P − 1 is the maximum length of x3 [n] resulting from the linear


convolution.
The N-point circular convolution of x1 [n] and x2 [n] is
N
X −1 N
X −1
x1 [n] N x2 [n] = x1 [m]x2 [((n − m))N ] = x1 [m]x̃2 [n − m] :
m=0 m=0

It is easy to see that the circular convolution product will be equal to the linear
convolution product on the interval 0 to N − 1 as long as we choose
N ≥ L + P − 1. The process of augmenting a sequence with zeros to make it
of a required length is called zero padding.

2.2 Convolution by sectioning

Suppose that for computational efficiency we want to implement a FIR system


using DFTs. It cannot in general be assumed that the input signal has a finite
duration, so the methods described up to now cannot be applied directly:

8
h[n]
0 P
x[n]

0 L 2L 3L
n
The solution is to use block convolution, where the signal to be filtered is
segmented into sections of length L. The input signal x[n], here assumed to be
causal, can be decomposed into blocks of length L as follows:

X
x[n] = xr [n − rL],
r=0

where 
x[n + rL] 0≤n≤L−1
xr [n] =
0 otherwise.
x0 [n]

0 L
n
x1 [n]

0 L
n
x2 [n]

0 L
n

9
The convolution product can therefore be written as

X
y[n] = x[n] ∗ h[n] = yr [n − rL],
r=0

where yr [n] is the response

yr [n] = xr [n] ∗ h[n].


y0 [n]

0 L+P-1
n
y1 [n]

0
n
y2 [n]

0
n

Since the sequences xr [n] have only L nonzero points and h[n] is of length P ,
each response term yr [n] has length L + P − 1. Thus linear convolution can
be obtained using N-point DFTs with N ≥ L + P − 1. Since the final result is
obtained by summing the overlapping output regions, this is called the
overlap-add method.
y[n]

0 L 2L 3L

An alternative block convolution procedure, called the overlap-save method,

10
corresponds to implementing an L-point circular convolution of a P-point
impulse response h[n] with an L-point segment xr [n]. The portion of the
output that corresponds to linear convolution is then identified (consisting of
L − (P − 1) points), and the resulting segments patched together to form the
output.

3 Spectrum estimation using the DFT

Spectrum estimation is the task of estimating the DTFT of a signal x[n]. The
DTFT of a discrete-time signal x[n] is

X

X(e ) = x[n]e−jωn .
n=−∞

The signal x[n] is generally of infinite duration, and X(ejω ) is a continuous


function of ω. The DTFT can therefore not be calculated using a computer.
Consider now that we truncate the signal x[n] by multiplying with the
rectangular window wr [n]:

1
w [n]

0.5
r

0
0 N−1
n
The windowed signal is then xw [n] = x[n]wr [n]. The DTFT of this windowed
signal is given by

X N
X −1
jω −jωn
Xw (e ) = xw [n]e = xw [n]e−jωn .
n=−∞ n=0

11
Noting that the DFT of xw [n] is
N −1
2πkn
X
Xw [k] = xw [n]e−j N ,
n=0

it is evident that
Xw [k] = Xw (ejω ) ω=2πk/N
.
The values of the DFT Xw [k] of the signal xw [n] are therefore periodic
samples of the DTFT Xw (ejω ), where the spacing between the samples is
2π/N . Since the relationship between the discrete-time frequency variable and
the continuous-time frequency variable is ω = ΩT , the DFT frequencies
correspond to continuous-time frequencies
2πk
Ωk = .
NT
The DFT can therefore only be used to find points on the DTFT of the
windowed signal xw [n] of x[n].
The operation of windowing involves multiplication in the discrete time
domain, which corresponds to continuous-time periodic convolution in the
DTFT frequency domain. The DTFT of the windowed signal is therefore
Z π
1
Xw (ejω ) = X(ejθ )W (ej(ω−θ) )dθ,
2π −π

where W (ejω ) is the frequency response of the window function. For a simple
rectangular window, the frequency response is as follows:

12
1

wr [n]
0.5
0
0 8
|Wr (ejω )|

0
−π 0 π
ω

The DFT therefore effectively samples the DTFT of the signal convolved with
the frequency response of the window.
Example: Spectrum analysis of sinusoidal signals Suppose we have the
sinusoidal signal combination

x[n] = cos(π/3n) + 0.75 cos(2π/3n), −∞ < n < ∞.

Since the signal is infinite in duration, the DTFT cannot be computed


numerically. We therefore window the signal in order to make the duration
finite:

13
1

x[n]
0
−1
0 8

1
wr [n]

0.5

0
0 8

1
xw [n]

0
−1
0 8
n

The operation of windowing modifies the signal. This is reflected in the


discrete-time Fourier transform domain by a spreading of the frequency
components:
|X(ejω )|

0.5

0
−π − 2π
3
− π3 0
π
3

3
π
|Wr (e )|

8

0
−π − 2π
3
− π3 0
π
3

3
π
|Xw (e )|

5

0
−π − 2π
3
− π3 0
π
3

3
π
ω

14
The operation of windowing therefore limits the ability of the Fourier
transform to resolve closely-spaced frequency components. When the DFT is
used for spectrum estimation, it effectively samples the spectrum of this
modified signal at the locations of the crosses indicated:
6
|X[k]|

0
0 1 2 3 4 5 6 7
k
Note that since k = 0 corresponds to ω = 0, there is a corresponding shift in
the sampled values.
In general, the elements of the N -point DFT of xw [n] contain N evenly-spaced
samples from the DTFT Xw (ejω ). These samples span an entire period of the
DTFT, and therefore correspond to frequencies at spacings of 2π/N . We can
obtain samples with a closer spacing by performing more computation.
Suppose we form the zero-padded length M signal xM [n] as follows:

x[n] 0≤n≤N −1
xM [n] =
0 N ≤ n ≤ M − 1.

The M -point DFT of this signal is


M −1 N −1
−j 2π 2π
X X
XM [k] = xM [n]e M kn = xw [n]e−j M kn
n=0 n=0


X
= xw [n]e−j M kn
n=−∞

The sample Xp [k] can therefore be seen to correspond to the DTFT of the the
windowed signal xw [n] at frequency ωk = 2πk/M . Since M is chosen to be
larger than N , the transform values correspond to regular samples of Xw (ejω )

15
with a closer spacing of 2π/M . The following figure shows the magnitude of
the DFT transform values for the 8-point signal shown previously, but
zero-padded to use a 32-point DFT:

1
xM [n]

0
−1

0 5 10 15 20 25 30
n
5
|X[k]|

0
0 5 10 15 20 25 30
k

Note that this process increases the density of the samples, but has no effect on
the resolution of the spectrum.
If W (ejω ) is sharply peaked, and approximates a Dirac delta function at the
origin, then Xw (ejω ) ≈ X(ejω ). The values of the DFT then correspond quite
accurately to samples of the DTFT of x[n]. For a rectangular window, the
approximation improves as N increases:

16
1

wr [n]
0.5

0
0 32
|Wr (ejω )|

32

0
−π 0 π
|Xw (ejω )|

0
−π 0 π
ω

The magnitude of the DFT of the windowed signal is

15
|X[k]|

10
5
0
0 5 10 15 20 25 30
k

which is clearly easier to interpret than for the case of the shorter signal. As
the window length tends to ∞, the relationship becomes exact.
The rectangular window inherent in the DFT has the disadvantage that the peak
sidelobe of Wr (ejω ) is high relative to the mainlobe. This limits the ability of
the DFT to resolve frequencies. Alternative windows may be used which have
preferred behaviour — the only requirement is that in the time domain the
window function is of finite duration. For example, the triangular window

17
1

wr [n]
0.5

0
0 32
|Wr (ejω )|

0
−π 0 π
|Xw (ejω )|

0
−π 0 π
ω

leads to DFT samples with magnitude


10
|X[k]|

0
0 5 10 15 20 25 30
k

Here the sidelobes have been reduced at the cost of diminished resolution —
the mainlobe has become wider.
The method just described forms the basis for the periodogram spectrum
estimate. It is often used in practice on account of its perceived simplicity.
However, it has a poor statistical properties — model-based spectrum
estimates generally have higher resolution and more predictable performance.
Finally, note that the discrete samples of the spectrum are only a complete
representation if the sampling criterion is met. The samples therefore have to

18
be sufficiently closely spaced.

4 Fast Fourier transforms


The widespread application of the DFT to convolution and spectrum analysis
is due to the existence of fast algorithms for its implementation. The class of
methods are referred to as fast Fourier transforms (FFTs).
Consider a direct implementation of an 8-point DFT:
7
X
X[k] = x[n]W8kn , k = 0, . . . , 7.
n=0

If the factors W8kn have been calculated in advance (and perhaps stored in a
lookup table), then the calculation of X[k] for each value of k requires 8
complex multiplications and 7 complex additions. The 8-point DFT therefore
requires 8 × 8 multiplications and 8 × 7 additions. For an N-point DFT these
become N 2 and N (N − 1) respectively. If N = 1024, then approximately one
million complex multiplications and one million complex additions are
required.
The key to reducing the computational complexity lies in the observation that
the same values of x[n]W8kn are effectively calculated many times as the
computation proceeds — particularly if the transform is long.
The conventional decomposition involves decimation-in-time, where at each
stage a N-point transform is decomposed into two N/2-point transforms. That
is, X[k] can be written as
N/2−1 N/2−1
(2r+1)k
X X
X[k] = x[2r]WN2rk + x[2r + 1]WN
r=0 r=0
N/2−1 N/2−1
X X
= x[2r](WN2 )rk + WNk x[2r + 1](WN2 )rk .
r=0 r=0

19
Noting that WN2 = WN/2 this becomes
N/2−1 N/2−1
X X
rk
X[k] = x[2r](WN/2 ) + WNk x[2r + 1](WN/2 )rk
r=0 r=0

= G[k] + WNk H[k].

The original N-point DFT can therefore be expressed in terms of two


N/2-point DFTs.
The N/2-point transforms can again be decomposed, and the process repeated
until only 2-point transforms remain. In general this requires log2 N stages of
decomposition. Since each stage requires approximately N complex
multiplications, the complexity of the resulting algorithm is of the order of
N log2 N .
The difference between N 2 and N log2 N complex multiplications can
become considerable for large values of N . For example, if N = 2048 then
N 2 /(N log2 N ) ≈ 200.
There are numerous variations of FFT algorithms, and all exploit the basic
redundancy in the computation of the DFT. In almost all cases an off-the-shelf
implementation of the FFT will be sufficient — there is seldom any reason to
implement a FFT yourself.

20

You might also like