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 O : L 2 → L 2 satisfying the two properties • Linearity O { α 1 x 1 + x 2 } = α 1 O { x 1 } + O { x 2 } • Time invariance � � T a { u } ( t ) = u ( t − a ) T a { x } = T a O { x } O ⇒ • Sound can be discretized…

Focus on the discrete, finite case • Linearity means linear operator L : R n → R n • Time invariance could mean different things • Most powerful definition uses circular shifts , with • � � � � T ( u ) = T L ( u ) L • T ( u ) = v , v i mod n = u i − 1 mod n     u 0 u n − 1   1 u 1 u 0     1       u 2 u 1 T T =       LT = TL ... =       . .       . . . .     1 u n − 1 u n − 2

What can an LTI system do? • Look at the effect of L on vector u n − 1 ! n − 1 n − 1 n − 1 X X X X T i ( e 0 ) � � u i T i � � u i L ( e i ) = v = L ( u ) = L = = L ( e 0 ) u i e i u i L i =0 i =0 i =0 i =0 • It is completely defined by ! h = L ( e 0 ) n − 1 n − 1 Impulse Impulse response X ⇥ T i ( h ) ⇤ X v j = u i u i h j − i mod n j = i =0 Circular convolution between u and h i =0 • What about matrix L ? Circulant matrix   h 0 h n − 1 h 1 · · · h 1 h 0 h 2 · · ·   ⇥ h T n − 1 ( h ) ⇤ L = T ( h )   = · · · . . . ...   . . . . . .   h n − 1 h n − 2 h 0 · · ·

Summary of LTI systems • Every discrete linear time-invariant system L can be expressed as a circular convolution • The circular convolution between , u , h ∈ R n denoted or , is v ∈ R n u ∗ h h ∗ u n − 1 n − 1 X X v j = u i h j − i mod n u j − k mod n h k = i =0 k =0 is the impulse response of L • h = L ( δ ) ( δ = e 0 ) • The corresponding matrix L is circulant • It commutes with the unit circular shift T

Computational cost too high • Same as a matrix-vector product! O ( n 2 ) • 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?

Diagonalize T • First notice that T has n distinct eigenvalues � � 1 − λ � � � � 1 − λ = ( − λ ) n + ( − 1) 1+ n = 0 � � det( T − λ I ) = � ... ... � � � � � � � 1 − λ � � T n = I ( λ k ) n = 1 2 π i λ k = ω k , k ∈ { 0 , 1 , . . . , n − 1 } , ω = e n ⇒ • Find the eigenvector basis   1 1 1 1 · · · Tv k = λ k v k ω 2 ω n − 1 1 1 ω · · ·   P =  . . . .  ... . . . . [ Tv k ] j = λ k [ v k ] j = [ v k ] j − 1 mod n √ n   . . . .   ω n − 1 ( ω 2 ) n − 1 ( ω n − 1 ) n − 1 1 [ v k ] j = ( λ k ) j = ( ω k ) j · · · n − 1 n − 1 ( 1 − ( ω k − j ) n = 0 , k 6 = j ( ω k − j ) p = X X ( ω k ) p ( ω j ) p 1 − ω k − j k v j = = v ∗ P ∗ P = PP ∗ = I n, k = j p =0 p =0

Simultaneously diagonalize L • Eigenvectors of L and T are the same! � Tv = λ v T ( Lv ) = LTv = λ ( Lv ) ⇒ LT = TL • Lv and v are eigenvectors of T associated to λ • But � � dim N ( T − λ I ) = 1 • So we must have Lv = β v • And therefore L is diagonalizable by same P !

The Discrete Fourier Transform • Given a vector , the Discrete Fourier v ∈ R n Transform of v , denoted , is the vector F ( v ) n − 1 1 v k e − 2 π ijk X v = P ∗ v ˆ v j = ˆ n √ n k =0 • The Inverse Discrete Fourier Transform, denoted , restores v F − 1 (ˆ v ) n − 1 1 2 π ijk X v = PP ∗ v v = P ˆ v j = ˆ n v k e √ n k =0

The Convolution Theorem • If u , v are two vectors in R n , then F ( u ⇤ h ) = F ( u ) � F ( h ) [ a � b ] j = a j b j (Hadamard product, element-wise product) • Proof F ( u ∗ h ) = P ∗ ( u ∗ h ) = P ∗ Lu = P ∗ LPP ∗ u = ( P ∗ LP ) F ( u ) ⇤ T � F ( u ) ⇥ β 0 = D F ( u ) = β n − 1 · · · = P ∗ h � F ( u ) = F ( h ) � F ( u ) 1 ⇤ T ⇥ 1 P ∗ h = P ∗ Le 0 = P ∗ LP ( P ∗ e 0 ) = D 1 · · · ⇤ T ⇥ β 0 β n − 1 = · · ·

Alternative to convolution O ( n ) u , ˆ u � ˆ h ˆ h ˆ P P ∗ O ( n 2 ) O ( n 2 ) F − 1 F u , h u ∗ h O ( n 2 )

Fast Fourier Transform (FFT) • Divide and conquer!   h 0 h 1     n − 1 . h =   . X ˆ h k ω − jk .   h j = n   h n − 2   k =0 h n − 1 n/ 2 − 1 n/ 2 − 1 X X h 2 k ω − j 2 k h 2 k +1 ω − j (2 k +1)   h 0 = + n n h 2 k =0 k =0     . h e =   . m − 1 m − 1 .   n ) − jk + ω − j X X h e k ( ω 2 h o k ( ω 2 n ) − jk   = h n − 4 n   h n − 2 k =0 k =0 m − 1 m − 1   h 1 X X h e k ω − jk + ω − j h o k ω − jk = m n m h 3   k =0 k =0   . h o =   . .   = ˆ h e j mod m + ω − j n ˆ h o   j mod n h n − 3   h n − 1

Alternative to convolution O ( n ) u , ˆ u � ˆ h ˆ h ˆ P P ∗ O ( n log n ) O ( n log n ) F − 1 F u , h u ∗ h O ( n 2 )

Demo!