Concepts and Algorithms of Scientific and Visual Computing Fast - PowerPoint PPT Presentation

Concepts and Algorithms of Scientific and Visual Computing –Fast Fourier Transform– CS448J, Autumn 2015, Stanford University Dominik L. Michels

Fast Fourier Transform (FFT) For practical application, given a signal x ∈ ℓ 2 ( Z ), we only take a finite number of elements into account. W.l.o.g. we describe x = ( x 0 ,..., x n − 1 ) with a tuple of length n , which is a power of two. We decompose x in elements with even respectively odd indexes: x e = ( x 0 , x 2 ,..., x n − 2 ) , x o = ( x 1 , x 3 ,..., x n − 1 ) . The corresponding polynomials are defined by n n 2 − 1 2 − 1 n − 1 � � � x j χ j , x 2 j χ j , x 2 j +1 χ j . p x ( χ ) = p x e ( χ ) = p x o ( χ ) = j =0 j =0 j =0

Fast Fourier Transform (FFT) The identity p x ( χ ) = p x e ( χ 2 ) + χ p x o ( χ 2 ) follows by simple calculation. Furthermore, for the primitive n -th root of unity Ω n := exp(2 π i / n ) holds Ω k − n = Ω k n n for all k ∈ Z , and therefore � p x e ( Ω 2 k n ) + Ω k n · p x o ( Ω 2 k if k ∈ { 0 ,..., n n ) 2 − 1 } p x ( Ω k n ) = 2 ,..., n − 1 } . p x g ( Ω 2 k − n n · p x o ( Ω 2 k − n if k ∈ { n ) + Ω k ) , n n

Fast Fourier Transform (FFT) According to this scheme, the Fourier transform of � p x (1) , p x ( Ω n ) ,..., p x ( Ω n − 1 � ˆ x = ) n can be computed from � n ) ,..., p x e ( Ω n − 2 � � n ) ,..., p x o ( Ω n − 2 � p x e (1) , p x e ( Ω 2 p x o (1) , p x o ( Ω 2 ˆ ˆ x e = ) , x o = ) n n leading to a recursive algorithm, the so-called Cooley-Tukey fast Fourier transform (FFT). This works for signals with lengths given by powers of two; see [Cooley, Tukey 1965]. It was extended later in [Bluestein 1970] to arbitrary n ∈ N .

Fast Fourier Transform (FFT) Cooley-Tukey fast Fourier transform for an input signal vector x = ( x 0 ,..., x n − 1 ). function FFT( x , Ω n ) begin if n = 1 then ζ 0 ← x 0 end else ( ϕ 0 ,...,ϕ n / 2 − 1 ) ← FFT( x e , Ω 2 x e ← ( x 0 , x 2 ,..., x n − 2 ) n ) ( ψ 0 ,...,ψ n / 2 − 1 ) ← FFT( x o , Ω 2 x o ← ( x 1 , x 3 ,..., x n − 1 ) n ) for k ← 0 to n / 2 − 1 do ζ k ← ϕ k + Ω k n · ψ k ζ n / 2+ k ← ϕ k + Ω n / 2+ k · ψ k n end end return ( ζ 0 ,...,ζ n − 1 ) end

Fast Fourier Transform (FFT) The runtime of the Cooley-Tukey fast Fourier transform can be described in dependence of the signal length n by T ( n ) = 2 · T ( n / 2) + O ( n ) with T (1) ∈ O (1). Iterative substitution leads to the time complexity T ( n ) ∈ O ( n log( n )) . Using the identity FFT − 1 ( ζ , Ω n ) = 1 n FFT( ζ , Ω − 1 n ) we also obtain the time complexity O ( n log( n )) for the inverse FFT.

Two-dimensional Discrete Fourier Transform For a given finite two-dimensional signal x : { 0 ,..., n 1 − 1 } × { 0 ,..., n 2 − 1 } → C , its Fourier transform ˆ x : { 0 ,..., n 1 − 1 } × { 0 ,..., n 2 − 1 } → C is given by n 1 − 1 n 2 − 1 � � �� 1 k 1 k 2 � � ˆ x ( k 1 , k 2 ) = x ( j 1 , j 2 )exp − 2 π i + j 2 √ n 1 n 2 j 1 n 1 n 2 j 1 =0 j 2 =0 as a natural extension of the one-dimensional case.

Two-dimensional Discrete Fourier Transform Moreover, the Fourier transform can be expressed using the equivalent representation n 1 − 1 n 2 − 1 � � �� 1 k 1 k 2 � � ˆ x ( k 1 , k 2 ) = x ( j 1 , j 2 )exp − 2 π i + j 2 √ n 1 n 2 j 1 n 1 n 2 j 1 =0 j 2 =0   n 2 − 1 n 1 − 1 � � 1 1 k 1 k 2 �  �    = x ( j 1 , j 2 )exp( − 2 π i j 1 )  exp − 2 π i j 2 √ n 2 √ n 1 ,      n 1  n 2    j 2 =0 j 1 =0 in which the inner term (in brackets) represents a column-wise Fourier transform and the outer term represents a row-wise Fourier transform. Each Fourier transform can be carried out using a FFT leading to the time complexity n 2 · O ( n 1 log( n 1 )) + n 1 · O ( n 2 log( n 2 )) = O ( n log( n )) with n := n 1 n 2 . A parallelization with c := max { n 1 , n 2 } processors leads to O ( c log( c )).

Multidimensional Discrete Fourier Transform More general, the Fourier transform of for a signal x : { 0 ,..., n 1 − 1 } × ··· × { 0 ,..., n d − 1 } → C of dimension d is defined by n 1 − 1 n d − 1  d  1 k l � � �   ˆ x ( k 1 ,..., k d ) = ... x ( j 1 ,..., j d )exp   − 2 π i j l  .      �� d n l   l =1 n l j 1 =0 j d =0 l =1 As in the one-dimensional case, the time complexity is given by O ( n log( n )) with n := � n l =1 n l using the FFT.