The Frequency Domain DS-GA 1013 / MATH-GA 2824 Mathematical Tools - PowerPoint PPT Presentation

Recovered Fourier coefficients ( N = 10) 0.5 0.4 0.3 Magnitude Signal Reconstruction 0.2 0.1 0.0 6 4 2 0 2 4 6 Frequency (Hz)

Recovered signal ( N = 10) 1.00 0.75 0.50 0.25 Signal Reconstruction 0.00 Samples 0.25 0.50 0.75 1.00 0.0 0.2 0.4 0.6 0.8 1.0 Time (s)

Sampling a real sinusoid x ( t ) := cos( 8 π t ) = 0 . 5 exp( − i 2 π 4 t ) + 0 . 5 exp( i 2 π 4 t ) N = 5 (as if k c = 2) � x rec [ k ] = ˆ x [ m ] ˆ { ( m − k ) mod 5 = 0 } x rec [ − 2 ] = 0 ˆ x rec [ − 1 ] = ˆ x rec [ 4 ] = 0 . 5 ˆ x rec [ 0 ] = 0 ˆ x rec [ 1 ] = ˆ x rec [ − 4 ] = 0 . 5 ˆ x rec [ 2 ] = 0 ˆ

Recovered Fourier coefficients ( N = 5) 0.5 0.4 0.3 Magnitude Signal Reconstruction 0.2 0.1 0.0 6 4 2 0 2 4 6 Frequency (Hz)

Recovered signal ( N = 5) 1.00 0.75 0.50 0.25 Signal Reconstruction 0.00 Samples 0.25 0.50 0.75 1.00 0.0 0.2 0.4 0.6 0.8 1.0 Time (s)

Aliasing Show videos

What happens if we sample too slowly? Let x be a signal that is with cut-off frequency k true / T We measure x [ N ] , N samples of x at 0, T / N , 2 T / N , . . . T − T / N What happens if we recover the signal assuming it is bandlimited with cut-off freq k samp / T , N = 2 k samp + 1, but actually k true > k samp ? x rec [ k ] := 1 N ( � F ∗ ˆ [ N ] x [ N ] )[ k ] � = x [ m ] ˆ { ( m − k ) mod N = 0 } This is called aliasing

Proof � � � i 2 π mj � N − 1 k true � � 1 [ N ] x [ N ] )[ k ] = 1 − i 2 π kj N ( � F ∗ exp x [ m ] exp ˆ N N N j = 0 m = − k true � � k true � = 1 ψ k , x [ m ] ψ m ˆ N m = − k true � = x [ m ] ˆ { ( m − k ) mod N = 0 }

Electrocardiogram: Fourier coefficients (magnitude) 10 1 2 10 Magnitude 10 3 10 4 60 40 20 0 20 40 60 Frequency (Hz)

Sampling an electrocardiogram Signal is approximately bandlimited at 50 Hz T = 8 s, so k c = 50 / ( 1 / T ) = 400 To avoid aliasing N ≥ 801

Recovered Fourier coefficients ( N =1,000) 10 1 10 2 Magnitude 10 3 10 4 Signal Reconstruction 10 5 80 60 40 20 0 20 40 60 80 Frequency (Hz)

Recovered signal ( N =1,000) Signal 1.5 Reconstruction Samples 1.0 Voltage (mV) 0.5 0.0 0.5 1.0 3.40 3.45 3.50 3.55 3.60 3.65 3.70 3.75 Time (s)

Sampling an electrocardiogram Signal is approximately bandlimited at 50 Hz T = 8 s, so k c = 50 / ( 1 / T ) = 400 N = 625 � x rec [ k ] = ˆ x [ m ] ˆ { ( m − k ) mod 625 = 0 } Component at m = ± 400 (50 Hz) shows up at ± 225 (28.1 Hz)

Recovered Fourier coefficients ( N = 625) 10 1 10 2 Magnitude 10 3 10 4 Signal Reconstruction 10 5 80 60 40 20 0 20 40 60 80 Frequency (Hz)

Recovered signal ( N = 625) Signal 1.5 Reconstruction Samples 1.0 Voltage (mV) 0.5 0.0 0.5 1.0 3.40 3.45 3.50 3.55 3.60 3.65 3.70 3.75 Time (s)

The frequency domain Sampling Discrete Fourier transform Frequency representations in multiple dimensions

Discrete complex sinusoids The discrete complex sinusoid ψ k ∈ C N with frequency k is � i 2 π kj � ψ k [ j ] := exp , 0 ≤ j , k ≤ N − 1 N √ N : orthonormal basis of C N Discrete complex sinusoids scaled by 1 /

ψ 2 (N=10)

ψ 3 (N=10)

Discrete Fourier transform The discrete Fourier transform (DFT) of x ∈ C N is   1 1 1 · · · 1 � � � � � �   − i 2 π ( N − 1 ) − i 2 π − i 2 π 2  1 exp exp · · · exp   N N N  � � � � � �   − i 2 π 2 ( N − 1 ) − i 2 π 2 − i 2 π 4 x := ˆ  1 exp exp · · · exp  x  N N N    · · · · · · · · · · · · · · ·  � � � � � �  − i 2 π ( N − 1 ) 2 − i 2 π ( N − 1 ) − i 2 π 2 ( N − 1 ) 1 exp exp · · · exp N N N = F [ N ] x x [ k ] = � x , ψ k � , ˆ 0 ≤ k ≤ N − 1

Inverse discrete Fourier transform y ∈ C N equals The inverse DFT of a vector ˆ y = 1 N F ∗ [ N ] ˆ y It inverts the DFT

Interpretation in terms of bandlimited signals If x ∈ C N contains samples of a bandlimited signal such that 2 k c + 1 ≤ N the DFT contains the Fourier series coefficients of the function x [ k c ] = 1 F ∗ � ˆ [ N ] x [ N ] N   · · · 1 1 1   � i 2 π ( − k c ) � � i 2 π ( − k c + 1 ) � � �   i 2 π k c exp exp · · · exp   N N N �   F [ N ] :=    · · · · · · · · · · · ·    � i 2 π ( − k c )( N − 1 ) � � i 2 π ( − k c + 1 )( N − 1 ) � � i 2 π k c ( N − 1 ) � exp exp · · · exp N N N Rows of � F [ N ] equal rows of F [ N ] in a different order!

Complexity of computing the DFT Complexity of multiplying N × N matrix with N -dim. vector is N 2 Very slow! We can exploit the structure of the matrix to do much better

Fast Fourier transform The most important numerical algorithm of our lifetime (G. Strang) Main insight: Action of N -order DFT matrix on vector can be decomposed into action of N / 2-order DFT submatrices on subvectors

Separation in even/odd columns and top/bottom rows x [0] ˆ x [0] � x [1] ˆ x [1] � x [2] ˆ � x [2] x [3] ˆ x [3] � = x [4] � x [4] ˆ x [5] � x [5] ˆ x [6] � x [6] ˆ x [7] � x [7] ˆ

Even columns can be scaled to yield odd columns e − 2 πi (0) / 8 e − 2 πi (1) / 8 e − 2 πi (2) / 8 e − 2 πi (3) / 8 = e − 2 πi (4) / 8 e − 2 πi (5) / 8 e − 2 πi (6) / 8 e − 2 πi (7) / 8

Top even submatrix and bottom even submatrix are both an N / 2-order DFT matrix =

FFT identity x [0] ˆ x [0] � e − 2 πi (0) / 8 � x [1] x [1] ˆ x [2] � e − 2 πi (1) / 8 � x [3] = + x [2] ˆ � x [4] � x [5] e − 2 πi (2) / 8 x [3] ˆ � x [6] x [7] � e − 2 πi (3) / 8 x [4] ˆ x [0] � � x [1] e − 2 πi (4) / 8 x [5] ˆ x [2] � � x [3] e − 2 πi (5) / 8 = + x [6] ˆ x [4] � e − 2 πi (6) / 8 x [5] � x [7] ˆ x [6] � e − 2 πi (7) / 8 x [7] �

Cooley-Tukey Fast Fourier transform 1. Compute F [ N / 2 ] x even . 2. Compute F [ N / 2 ] x odd . 3. For k = 0 , 1 , . . . , N / 2 − 1 set � � − i 2 π k F [ N ] x [ k ] := F [ N / 2 ] x even [ k ] + exp F [ N / 2 ] x odd [ k ] , N � � − i 2 π k F [ N ] x [ k + N / 2 ] := F [ N / 2 ] x even [ k ] − exp F [ N / 2 ] x odd [ k ] . N

Complexity DFT 2 L DFT 2 L − 1 DFT 2 L − 1 DFT 2 L − 2 DFT 2 L − 2 DFT 2 L − 2 DFT 2 L − 2 DFT 2 L − 3 DFT 2 L − 3 DFT 2 L − 3 DFT 2 L − 3 DFT 2 L − 3 DFT 2 L − 3 DFT 2 L − 3 DFT 2 L − 3

Complexity Assume N = 2 L L = log 2 N levels At level m ∈ { 1 , . . . , L } there are 2 m nodes At each node, scale a vector of dim 2 L − m and add to another vector Complexity at each node: 2 L − m Complexity at each level: 2 L − m 2 m = 2 L = N Complexity is O ( N log N ) !

In practice 10 1 DFT (matrix) FFT (recursive) 10 2 Running time (s) 10 3 10 4 10 5 10 2 10 3 10 4 N

The frequency domain Sampling Discrete Fourier transform Frequency representations in multiple dimensions

Multidimensional signals Square-integrable functions defined on a hyperrectangle I := [ a 1 , b 1 ] × . . . × [ a p , b p ] ⊂ R p Inner product: � � x , y � := x ( t ) y ( t ) d t . I Goal: Extension of frequency representations to multidimensional signals

Multidimensional sinusoid a cos ( 2 π � f , t � + θ ) . The frequency and time indices are now d -dimensional Periodic with period 1 / || f || 2 in direction of f For any integer m � � � � m f a cos 2 π f , t + + θ = a cos ( 2 π � f , t � + 2 π m + θ ) || f || 2 || f || 2 = a cos ( 2 π � f , t � + θ )

2D sinusoid 1.0 1.00 0.75 0.8 0.50 0.25 0.6 0.00 t 1 0.4 0.25 0.50 0.2 0.75 0.0 1.00 0.0 0.2 0.4 0.6 0.8 1.0 t 2

Multidimensional complex sinusoids Complex sinusoid with frequency f ∈ R d : exp( i 2 π � f , t � ) := cos( 2 π � f , t � ) + i sin( 2 π � f , t � ) . cos ( i 2 π � f , t � + θ ) = exp( i θ ) exp( i 2 π � f , t � ) + exp( − i θ ) exp( − i 2 π � f , t � ) 2 2

Multidimensional complex sinusoids Can be expressed as product of 1D complex sinusoids   d �  i 2 π  exp( i 2 π � f , t � ) := exp f [ j ] t [ j ] j = 1 d � = exp( i 2 π f [ j ] t [ j ]) j = 1 From now on d = 2: t [ 1 ] = t 1 , t [ 2 ] = t 2

Orthogonality of multidimensional complex sinusoids The family of complex sinusoids with integer frequencies � i 2 π k 1 t 1 � � i 2 π k 2 t 2 � φ 2D k 1 , k 2 ( t 1 , t 2 ) := exp exp , k 1 , k 2 ∈ Z , T T is an orthogonal set of functions on any interval of the form [ a , a + T ] × [ b , b + T ] , a , b , T ∈ R and T > 0

Proof We have φ 2D k 1 , k 2 ( t 1 , t 2 ) = φ k 1 ( t 1 ) φ k 2 ( t 2 ) , so that � a + T � b + T � � φ 2D k 1 , k 2 , φ 2D = φ k 1 ( t 1 ) φ k 2 ( t 2 ) φ j 1 ( t 1 ) φ j 2 ( t 2 ) d t 1 d t 2 j 1 , j 2 t 1 = a t 2 = b = � φ k 1 , φ j 1 � � φ k 2 , φ j 2 � = 0 as long as j 1 � = k 1 or j 2 � = k 2

φ 2D 0 , 5 + φ 2D 0 , − 5 10 9 8 7 6 5 4 3 2 1 0 k 1 1 2 3 4 5 6 7 8 9 10 10 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 k 2

φ 2D 0 , 5 + φ 2D 0 , − 5 1.0 1.00 0.75 0.8 0.50 0.25 0.6 0.00 t 1 0.4 0.25 0.50 0.2 0.75 0.0 1.00 0.0 0.2 0.4 0.6 0.8 1.0 t 2

φ 2D 10 , 0 + φ 2D − 10 , 0 10 9 8 7 6 5 4 3 2 1 0 k 1 1 2 3 4 5 6 7 8 9 10 10 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 k 2

φ 2D 10 , 0 + φ 2D − 10 , 0 1.0 1.00 0.75 0.8 0.50 0.25 0.6 0.00 t 1 0.4 0.25 0.50 0.2 0.75 0.0 1.00 0.0 0.2 0.4 0.6 0.8 1.0 t 2

φ 2D 3 , 4 + φ 2D − 3 , − 4 10 9 8 7 6 5 4 3 2 1 0 k 1 1 2 3 4 5 6 7 8 9 10 10 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 k 2

φ 2D 3 , 4 + φ 2D − 3 , − 4 1.0 1.00 0.75 0.8 0.50 0.25 0.6 0.00 t 1 0.4 0.25 0.50 0.2 0.75 0.0 1.00 0.0 0.2 0.4 0.6 0.8 1.0 t 2

φ 2D 8 , − 6 + φ 2D − 8 , 6 10 9 8 7 6 5 4 3 2 1 0 k 1 1 2 3 4 5 6 7 8 9 10 10 9 8 7 6 5 4 3 2 1 0 1 2 3 4 5 6 7 8 9 10 k 2

φ 2D 8 , − 6 + φ 2D − 8 , 6 1.0 1.00 0.75 0.8 0.50 0.25 0.6 0.00 t 1 0.4 0.25 0.50 0.2 0.75 0.0 1.00 0.0 0.2 0.4 0.6 0.8 1.0 t 2

2D Fourier series The Fourier series coefficients of a function x ∈ L 2 [ a , a + T ] for any a , T ∈ R , T > 0, are given by � � x , φ 2D x [ k 1 , k 2 ] := ˆ k 1 , k 2 � a + T � b + T � � � � − i 2 π k 1 t 1 − i 2 π k 2 t 2 = x ( t 1 , t 2 ) exp exp d t 1 d t 2 T T t 1 = a t 2 = b The Fourier series of order k c , 1 , k c , 2 is defined as k c , 1 k c , 2 � � F k c , 1 , k c , 2 { x } := 1 x [ k 1 , k 2 ] φ 2D ˆ k 1 , k 2 . T k 1 = − k c , 1 k 2 = − k c , 2

Magnetic resonance imaging Non-invasive medical-imaging technique Measures response of atomic nuclei in biological tissues to high-frequency radio waves when placed in a strong magnetic field Radio waves adjusted so that each measurement equals 2D Fourier coefficients of proton density of hydrogen atoms in a region of interest