lgebra Linear e Aplicaes DISCRETE FOURIER TRANSFORM (AND THE FFT) - - PowerPoint PPT Presentation

Álgebra Linear e Aplicações

DISCRETE FOURIER TRANSFORM

(AND THE FFT)

Motivation

• Consider the propagation of sound from fixed

source to fixed destination in an environment

• Sound is pressure variation through time
• It is a vector space
• Each environment can be seen as an operator
• Different environments, different operators
• What are the properties we can expect?

Linear, time-invariant systems

• An LTI system is an operator

satisfying the two properties

• Linearity
• Time invariance
• Sound can be discretized…

O : L2 → L2

O{α1x1 + x2} = α1O{x1} + O{x2} Ta{u}(t) = u(t − a) ⇒ O

• Ta{x}

= Ta

• O{x}

Focus on the discrete, finite case

• Linearity means linear operator
• Time invariance could mean different things
• Most powerful definition uses circular shifts
• , with
• T(u) = v,

vi mod n = ui−1 mod n

T =      1 1 ... 1      T        u0 u1 u2 . . . un−1        =        un−1 u0 u1 . . . un−2       

L

• T(u)
• = T
• L(u)
• L : Rn → Rn

LT = TL

What can an LTI system do?

• Look at the effect of L on vector u
• It is completely defined by !
• What about matrix L?

vj =

n−1

X

i=0

ui ⇥ T i(h) ⇤

j = n−1

X

i=0

ui hj−i mod n =      h0 hn−1 · · · h1 h1 h0 · · · h2 . . . . . . ... . . . hn−1 hn−2 · · · h0     

Circulant matrix Impulse Circular convolution between u and h Impulse response

v = L(u) = L n−1 X

i=0

uiei ! =

n−1

X

i=0

ui L(ei) =

n−1

X

i=0

ui L

• T i(e0)
• =

n−1

X

i=0

ui T i L(e0)

• h = L(e0)

L = ⇥h T(h) · · · T n−1(h)⇤

Summary of LTI systems

• Every discrete linear time-invariant system L

can be expressed as a circular convolution

• The circular convolution between ,

denoted or , is

• is the impulse response of L
• The corresponding matrix L is circulant
• It commutes with the unit circular shift T

h = L(δ) (δ = e0)

vj =

n−1

X

i=0

ui hj−i mod n

u, h ∈ Rn u ∗ h v ∈ Rn h ∗ u

=

n−1

X

k=0

uj−k mod n hk

Computational cost too high

• Same as a matrix-vector product!
• But matrix has only n degrees of freedom
• Doesn’t make sense
• Vectors can be very large
• Millions of entries
• Need a faster alternative
• What if we found a better basis?

O(n2)

Diagonalize T

• First notice that T has n distinct eigenvalues
• Find the eigenvector basis

det(T − λI) =

• −λ

1 1 −λ ... ... 1 −λ

• = (−λ)n + (−1)1+n = 0

Tn = I ⇒ (λk)n = 1 Tvk = λkvk λk = ωk, k ∈ {0, 1, . . . , n − 1}, ω = e

2πi n

[Tvk]j = λk[vk]j = [vk]j−1 mod n [vk]j = (λk)j = (ωk)j v∗

kvj = n−1

X

p=0

(ωk)p(ωj)p P∗P = PP∗ = I =

n−1

X

p=0

(ωk−j)p = (

1−(ωk−j)n 1−ωk−j

= 0, k 6= j n, k = j P = 1 √n      1 1 1 · · · 1 1 ω ω2 · · · ωn−1 . . . . . . . . . ... . . . 1 ωn−1 (ω2)n−1 · · · (ωn−1)n−1     

Simultaneously diagonalize L

• Eigenvectors of L and T are the same!
• Lv and v are eigenvectors of T associated to
• But
• So we must have
• And therefore L is diagonalizable by same P!

λ dim

• N(T − λI)
• = 1

Lv = βv

Tv = λv LT = TL

• ⇒

T(Lv) = LTv = λ(Lv)

The Discrete Fourier Transform

• Given a vector , the Discrete Fourier

Transform of v, denoted , is the vector

• The Inverse Discrete Fourier Transform,

denoted , restores v

v ∈ Rn F(v) ˆ v = P∗v

ˆ vj = 1 √n

n−1

X

k=0

vke− 2πijk

n

F −1(ˆ v) v = Pˆ v = PP∗v

vj = 1 √n

n−1

X

k=0

ˆ vke

2πijk n

The Convolution Theorem

• If u,v are two vectors in Rn, then
• Proof

F(u ⇤ h) = F(u) F(h) [a b]j = ajbj

(Hadamard product, element-wise product)

F(u ∗ h) = P∗(u ∗ h) = P∗Lu = P∗LPP∗u = (P∗LP)F(u) = D F(u) P∗h = P∗Le0 = P∗LP(P∗e0) = ⇥β0 · · · βn−1 ⇤T F(u) = P∗h F(u) = F(h) F(u) = ⇥β0 · · · βn−1 ⇤T = D ⇥1 1 · · · 1⇤T

Alternative to convolution

u, h u ∗ h O(n2) ˆ u ˆ h O(n) ˆ u, ˆ h P∗ F O(n2) P F −1 O(n2)

Fast Fourier Transform (FFT)

• Divide and conquer!

ˆ hj =

n−1

X

k=0

hkω−jk

n

=

n/2−1

X

k=0

h2k ω−j2k

n

+

n/2−1

X

k=0

h2k+1 ω−j(2k+1)

n

=

m−1

X

k=0

he

k (ω2 n)−jk + ω−j n m−1

X

k=0

ho

k (ω2 n)−jk

=

m−1

X

k=0

he

k ω−jk m

+ ω−j

n m−1

X

k=0

ho

k ω−jk m

= ˆ he

j mod m + ω−j n ˆ

ho

j mod n

h =        h0 h1 . . . hn−2 hn−1        he =        h0 h2 . . . hn−4 hn−2        ho =        h1 h3 . . . hn−3 hn−1       

Alternative to convolution

u, h u ∗ h ˆ u, ˆ h ˆ u ˆ h P P∗ F F −1 O(n2) O(n) O(n log n) O(n log n)

Demo!