Your Federal Quarterly Tax Payments are due April 15th

# MAGIC015 Lecture 9 The Discrete Fourier Transform by qxc16070

VIEWS: 13 PAGES: 23

• pg 1
```									         MAGIC015 Lecture 9: The Discrete Fourier Transform

Jeremy Levesley jl1@mcs.le.ac.uk

University of Leicester

November 5, 2009

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   1 / 23
Introduction

We receive some sort of signal

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   2 / 23
Fourier approximation

In signal processing the need to rapidly compute approximate Fourier series:

SN f (x) =                    αk (f )exp(ikx),
|k|≤N

with
2π
1
αk (f ) =                         f (x) exp(−ikx)dx,
2π      0
is fundamental.
Approximate using all Newton-Cotes rules

2π                                             2N−1
2π                    2πj                      2πijk
f (x) exp(−ikx)dx ≈                                f                exp −                 ,    |k| ≤ N.
0                                           2N                    2N                        2N
j=0

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                             November 5, 2009   3 / 23
Discrete Fourier transform

Thus it is important to compute the
ˆ
Deﬁnition 9.1 Discrete Fourier transform (DFT) v of a vector v of length N:
N−1
2πijk
ˆ
vk =             vj exp −                        ,     k = 0, 1, · · · , N − 1.
N
j=0

Remark 9.2 We will view the input v as being periodic (vk = vN+k ).

History: Famous fast Fourier transform (FFT) algorithm by Cooley and Tukey
1965. First used by Gauss (1805).

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                        November 5, 2009   4 / 23
In this lecture we will look at:
1  Properties of the DFT;
2  FFT algorithm.
3  Smoothing and compression.
4  What is wrong with Fourier series.
5  The wavelet paradigm.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   5 / 23
Fourier basis
The Fourier basis,
2πijk
fj (k) = exp                          ,     k = 0, · · · , N − 1,
N
j = 0, 1, · · · , N − 1.
Lemma 9.3
The sum of the roots of unity is 0.
Proposition 9.4
The Fourier basis is orthogonal.
Proof
N−1
2πijk                     2πilk
f j , fl     =              exp                     exp −
N                         N
k=0
N−1
2πi(j − l)k
=              exp
N
k=0
N/d−1
2πik(j − l)/d
= d                   exp                                = 0,
N/d
k=0
MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                        November 5, 2009   6 / 23
Expanding in the Fourier basis
We make an orthonormal basis noting that

fj , fj = N,             j = 0, 1, · · · , N − 1.

Remark 9.5
f
√j ,          j = 0, 1, · · · , N − 1,
N
is an orthonormal basis for RN .
We can write any v ∈ RN in the form
N−1
1
v=                    v, fj fj ,
N
j=0

where
N−1
2πijk
v, fj =               vj exp −                          ˆ
= vj .
N
k=0

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                    November 5, 2009   7 / 23
The Inverse Fourier transform

Remark 9.6
The Fourier transform computes the coeﬃcients of a vector with respect to the
Fourier basis.

N−1
1
v=                  ˆ
vj fj ,
N
j=0
means
N−1
1                          2πijk
vk =                   vj exp
ˆ                   ,            k = 0, 1, · · · , N − 1.
N                           N
j=0

So we recover the original coeﬃcients via the inverse transform above.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                        November 5, 2009   8 / 23
Parseval and Plancherel

Inner product is independent of orthonormal basis:
Theorem 9.7 (Plancherel’s theorem)
1
v, u =          ˆ ˆ
v, u ,            v, u ∈ RN .
N

Corollary 9.8 (Parseval’s theorem)
1
v =√ v ,
ˆ                       v ∈ RN .
N

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   9 / 23
Matrix notation

We can write the DFT in the form
                                                                                                       
1    1        1                                               ···                           1
 1
      w        w2                                               ···                     · · · w N−1    

v= 1
ˆ
     w2        w4                                               ···                      w 2(N−1)       v,

 .     .        .                                               ..                             .
 .     .        .                                                                              .       
.    .        .                                                  .                           .       
1 w N−1 w 2(N−1)                                       w 3(N−1) · · ·           w (N−1)(N−1)

where w = exp − 2πi .
N

Write v = WN v, where WN (i, j) = w (i−1)(j−1) .
ˆ

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                               November 5, 2009   10 / 23
For example

For N = 4, w = exp −iπ/4 = −i. Then

                          
1 1                      1  1
 1 w                      w2 w3 
W4        = 
 1 w2

w4 w6 
1 w3                     w6 w9
                                
1 1                       1  1
 1 −i                     −1 i 
= 
 1 −1
.
1 −1 
1 i                      −1 −i

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   11 / 23
Filtering and convolution
Fourier transform gives frequency information. In audio engineering we wish to
manipulate the frequencies.
Can we arrange:
vj → uj = wj vj ,
ˆ    ˆ    ˆˆ                     j = 0, 1, · · · , N − 1.

Solved by
Deﬁnition 9.
The convolution of v, w ∈ RN is
N−1
[v ∗ w]j          =              vk wj−k ,          j = 0, 1, · · · , N − 1.
k=0

Theorem 9.

[v ∗ w]j = vj wj ,
ˆˆ                    j = 0, 1, · · · , N − 1.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                       November 5, 2009   12 / 23
Noise reduction
Noisy signal has high frequency. Filter the high frequency.
DFT v;
1

ˆ          ˆ
Make w by setting vj = 0, j > M;
2

3        ˆ
IDFT w .
Example 9.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   13 / 23
Compression
Smooth signal has few low frequencies. Filter the small frequencies.
1   DFT v;
2   Make w by setting vj = 0, if |ˆj | < ;
ˆ            ˆ           v
3        ˆ
IDFT w .

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   14 / 23
FFT idea

First reorder rows and columns:
                                
1 1      1   1       1 1  1  1
 1 −i −1 i 
 →  1 1 −1 −1  .
            

 1 −1 1 −1           1 −1 −i  i 
1 i     −1 −i        1 −1 i  −i

We need to keep track of changes in inputs and outputs - Binary Bit Reverse order
                                                                                               
1 1  1  1                                         1 1             0 0     1 1  0                0
 1 1 −1 −1                                       1 −1            0 0  0 0   1                1 
            =                                                                                  ,
 1 −1 −i  i                                      0 0             1 1   1 −1 0                0 
1 −1 i  −i                                        0 0             1 −1    0 0 −i                i
which is a sparse factorisation.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   15 / 23
Decimation formula
Theorem 9.10 Let N = 2M. Then
M−1
2πijk
v2j
ˆ        =              (vk + vk+M ) exp −
M
k=0
M−1
πik                                           2πijk
v2j+1
ˆ            =               exp                  (vk − vk+M ) exp −                            ,
M                                              M
k=0

j = 0, 1, · · · , M − 1.
Proof

N−1
2πi(2j + 1)k
v2j+1
ˆ             =              vk exp −
N
k=0
M−1
2πi(2j + 1)k                                            2πi(2j + 1)k + M
=                  vk exp −                                      + vk+M exp −
2M                                                         2M
k=0
M−1
πik                                            2πijk
=              exp −                    (vk − vk+M ) exp −                               .
M                                               M
k=0

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                              November 5, 2009   16 / 23
Complexity

Remark 9.11
We have replaced a length N DFT with 2 N/2 DFTs.
Number of complex multiplications to compute DFT is O(N 2 ).
After decimation we have 2 times (N/2)2 = N 2 /2, so a saving of 1/2. Now we
iterate.
The cost has been N/2 complex multiplications.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   17 / 23
The FFT

Suppose N = 2m . We can do m = log2 N decimations. Each takes N/2 complex
operations.
2 times N/2 to N/4 is 2 times N/4 complex multiplications = N/2
multiplications.
Total complexity
N
log2 N,
2
which brings large saving.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   18 / 23
Translation operator

Deﬁnition 9.9
The translation operator Tk : RN → RN is deﬁned by

[Tl v]j = v (j − l),                j = 0, · · · , N − 1.

Periodic extention of v.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   19 / 23
Lemma 9.10

2πijl
[Tˆv]j = exp −
l                                        ˆ
vj .
N
Proof

N−1
2πijk
[Tˆv]j
l            =             [Tl v]k exp −
N
k=0
N−1
2πijk
=              v (k − l) exp −
N
k=0
N−1
2πij(k + l)
=              v (k) exp −
N
k=0
2πijl
=      exp −                      ˆ
vj .
N

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   20 / 23
Why wavelets?
Theorem 9.
Suppose {Tk v : k = 0, 1, · · · , N − 1} is an orthogonal basis for RN . Then
|ˆk | = constant,
v                                 ∀k.

Proof
We have, using Plancherel

Tk v, Tj v = Tk v, Tj v = 0,                            k = j.

Using Lemma 9. we see that
N−1
2πikl                      2πijl
Tk v, Tj v            =              exp −                      ˆ
vl exp                ˆ
vl
N                          N
l=0
N−1
2πi(k − j)l
=              exp −                                  |vl |2 = 0.
N
l=0

|ˆl |2 is orthogonal to fj (Fourier basis) for j = 1, 2, · · · , N − 1.
v
v parallel to f0 = [1, 1, · · · , 1]T .
ˆ
MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                               November 5, 2009   21 / 23
Uncertainty principle

Localisation in space means no localisation in frequency.
Next step, use translates of two vectors:

RN = span({T2k u} ∪ {T2k v} : k = 0, 1, · · · , N/2 − 1).

This is the Wavelet Paradigm.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   22 / 23
Summary

In this lecture we have:
1   Properties of the DFT;
2   FFT algorithm.
3   Smoothing and compression.
4   What is wrong with Fourier series.
5   The wavelet paradigm.

MAGIC015 Lecture 9: The Discrete Fourier Transform
Jeremy Levesley jl1@mcs.le.ac.uk (University of Leicester)                                                 November 5, 2009   23 / 23

```
To top