Sampling and Fourier Transform Applications

420 5. Sampling and Discrete Fourier Transform Section 5.1 Sampling and Reconstruction Section 5.2 Discrete Fourier Transform [1] R. N. Bracewell, The Fourier Transform and Its Applications, 3rd ed., McGraw Hill, Boston, 2000. [2] A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, London: Prentice-Hall, 3rded., 2010.

421 Sampling and Discrete Fourier Transform Varying the Sampling Rate (Sec. 5-1-4) Sampling (Sec. 5-1-2) Reconstruction (Sec. 5-1-3) Sampling and Discrete Fourier Transform Implementing the Continuous FT (Sec. 5-2-2) Discrete Fourier Transform (Sec. 5-2-1) Properties (Sec. 5-2-3) Complexity (Sec. 5-2-4) 2D Version (Sec. 5-2-5)

422 5.1 Sampling 5.1.1 Impulse Train Impulse Train ( ) p x ( ) ( ) ( ) x ( ) ( ) = = + + + + + + 1 1 2 x n x x x n It is also called the comb function. .......... ..... x=2 x=0 x=1 x=-1 x=-2

423 Signal sampling can be express in terms of the impulse train ( ) ( ) ( ) x ( ) sampling = g x g g x x n s x n ( ) = g x n n x ( ) = n g g n where n x Since x x 1 1 ( ) = = x n n p x x x x x n n the sampled signal gs(x) can be expressed in terms of g x p x 1 ( ) x ( ) = g s x x

424 [Theorem 5.1.1] The impulse train is also an eigenfunction of the Fourier transform, i.e., ( ) P f ( ) p x ( ) = = f n n (Proof): Note that the impulse train is a periodic function ( ) ( ) 1 p x p x = + Therefore, it can be expanded by the Fourier series (page 321) of the complex form with T = 1 = ( ) p x ( ) exp 2 c j nx n n where 1 1 1/2 1/2 ( ) p x ( ) ( ) x ( ) 1/2 1/2 = nx dx = nx dx = exp 2 exp 2 1 nc j j (page 342(2))

425 Therefore, = ( ) p x ( ) exp 2 j nx n ( ) p x ( ) = F F exp 2 j nx n ( ) ( ) p f = = f n n ( ) ( ) = x n f n n n -1 0 1 2 -1 0 1 2

426 Varying the interval of the impulse train x x 1 1 ( ) x n = = n p x x x x x n n = x 1 ( ) ( ) (page 349(5)) x n = F F p p f x x x x n n 1 ( ) = = f n f x x x n n 1 x 1 2 1 2 0 0 x x x x x x

427 5.1.2 Sampling Theory [Theorem 5.1.2] Suppose that we perform sampling for a continuous signal with sampling interval x ( ) ( ) s g x g x = ( ) ( ) sampling g x x n x n ( ) = where g x n n x then ( ) = n g g n n x n 1 ( ) f = G G f s x x n p x ( )1 g x ( ) x (Proof): Since = g s x x x 1 ( ) x ( ) = F F F g g x p s x x n n 1 1 ( ) f ( ) ( ) ( ) = = = G G f p f G f f G f s x x x x x n n (page 387)

428 1 (fsis call the sampling frequency) f = ( ) s G If we set s x ( ) = then f f G f nf s s n ( ) G f ( ) 0 G f-axis f=0 ( ) = n f n f ( ) x ( ) g g n = = g g x x g x n x s n sampling s s n n ( ) f ( ) 0 ( ) = G f G f nf s s s ( ) f G n s sf G -fs f=0 2fs fs f-axis

429 [Sampling Theory] 1 The sampling frequency should be larger than twice of the f = s x bandwidth of the original continuous function: sf B B (Nyquist criterion) 2 sf B where ( ) G f = when f > B. 0 Otherwise, the original function cannot be reconstructed and the aliasing effect is led. B fs B B -fs f=0 fs+B 2fs fs f-axis

430 = 2 ? B sf Q: What will happen if Even component with frequency f = B preserved Odd component with frequency f = B destroyed

431 5.1.3 Reconstruction (Digital toAnalogous) When the Nyquist criterion is satisfied, one can apply the lowpass filter to reconstruct the original signal. ( ) ( ) 0 G f Frequency Domain ( ) G f G f-axis ( ) f G ( ) 0 sf G s ( ) ( ) f = G f G f nf s s s n B fs B fs f-axis -fs f=0 ( ) f G ( ) 0 f sf G 1 f ( ) ( ) f s = G f G 2 f s s c where B < fc< fs B B fs B fs -fs f=0 -fc fc

432 Time Domain Frequency Domain ( ) ( ) G f g x sampling ( ) = g g n ( ) ( ) f = G f G f nf n x s s s n f ( ) x = g g x n s n s n reconstruction f 1 f ( ) ( ) f = G f G 2 f s s c 2 f ( ) ( ) ( ) x where B < fc< fs B = sinc 2 g x g f x c f s c s

433 2 f ( ) ( ) ( ) = sinc 2 ( ) g x g f x d c f s c s f 2 n f ( ) = sinc 2 ( ) g f x d c f n c s s n 2 f n f ( ) = sinc 2 ( ) g x g f x c f n c s s n = Specially, when / 2 f f c s ( ) ( ) = sinc g x g f x n n s n Signal Reconstruction Formula: ( ) = g g n x where ( ) = sinc g x g n n x n x n constraint: 1 2 B x

434 [Example 1] Suppose that ( ) 2 = n = g g n = = n g = otherwise 1, 2, 0 g g g 1 1 0 ( ) for f 1 G f = 0 Try to reconstruct g(x). x= 1/2 (Solution): ( ) g x ( ) ( ) = + + + sinc 2 1 2sinc2 sinc 2 1 x x x

435 ( ) g x ( ) ( ) = + + + sinc 2 1 2sinc2 sinc 2 1 x x x ( ) x+ sinc 2 1 ( ) 2sinc 2x ( ) x sinc 2 1 ( ) g x ( ) ( ) ( ) + sinc 2 1 , 2sinc 2 , sinc 2 1 (Note): x x x do not interfere with one another at x = n/2.

436 5.1.4 Varying the Sampling Rate (i) D/A conversion x ( ) = sinc g x g n n x n (ii) Re-sampling n ( ) = = g sinc g n g m new n new m x m = k and k is an integer Note: When new x = n g g kn

437 From the view point of the spectrum, if ( ) ( ) x ( ) ( ) x = = g g x n g g x n s n new s n x n n then ( ) ( ) f ( ) ( ) f s G = G f G f nf = f G f nf s s s new new n n = = where 1/ , 1/ f f s x new new ( ) f ( ) 0 G sf G s B fs B ( ) 0 G fs f-axis -fs f=0 ( ) f G f s new fnew B fnew B f-axis -fnew f=0

438 5.2 Discrete Fourier Transform 5.2.1 Fourier Transform Derivation and Definitions of the Discrete To process discrete functions, the continuous Fourier transform should be converted into the discrete version. Continuous Fourier transform: ( ) ( ) = 2 j fx G f g x e dx If we set = = , f m x n f x then ( ) ( ) 2 j mn = G m g n e f x f x x n

439 ( ) ( ) 2 j mn = G m g n e f x f x x n Specially, if = 1/ N f x then ( ) ( ) 2 mn N j = G m g n e f x x n

440 Discrete Fourier Transform (DFT) 1 N = 2 mn N j = G m DFT g n g n e = 0 n m, n = 0, 1, 2, , N-1 Inverse Discrete Fourier Transform (IDFT) 1 N = 1 N 2 mn N j = g n IDFT G m G m e = 0 m Note that the output of g[n] is periodic = + G m G m N Also note that, on page 439, ( G m ) ( ) ( ) = DFT g n f x x

441 The DFT and the IDFT form a transform pair since 1 1 N N 1 1 1 N N N 1 N 1 N 1 N g k g k e 2 2 mk N mn N 2 2 2 mn N mk N mn N j j j j j = = e e G m e e = = = = = 0 0 k m 0 0 0 m m k Because 1 N 2 a j N 2 j a 1 1 1 e e e 2 a j m N = = = 0 e if a 0, N 2 2 a a j j 1 e N N = 0 m 1 N 2 a if a = 0, j m = e N N = unit impulse function (discrete Dirac delta function) 0 m 1 N 2 m n k ( ) j = n k e N N d = 0 m 1 1 N N 1 N 1 N g k N 2 mn N j = n k = G m e g n d = = 0 0 m k

442 Unit Impulse Function (discrete Dirac delta function) = 1 0 0 0 if n if n n = d 1 n -3 -2 -1 0 1 2 3 The unit impulse function has an explicit form. It does not a limitation of a distribution. n = Compared to page 341, 1 d n n = = 0 0 n if n d n d d

443 Other possible definitions of the DFT + 1 n N 0 = 2 mn N j G m g n e DFT n n = 0 + 1 m N 0 1 N 2 mn N j = IDFT g n G m e m m = 0 1 N 1 N 2 mn N j = G m g n e DFT = 0 n 0 1 N 1 N 2 mn N j = g n G m e IDFT = 0 m

444 5.2.2 Implementing Continuous FT by the DFT Suppose that we want to calculate the continuous FT of g(x) digitally and ( ) g x = for x [x1, x1+T] 0 (i) Shifting ( ) ( ) = + g x g x x 1 1 Note: ( ) g x = 0 for x [0, T] 1 (ii) Sampling n ( ) = g g n 1 d x (iii) DFT 1 N = n e 2 mn N j G m g d d = 0 n

445 (iv) Mapping to the true frequency ( 1 G m ) = G m (from pages 439 and 440) f d x Since f N 1 = = s = 1/ N N f f x x Therefore, f N = G m G m s if 0 m N/2, 1 d x f N = G m f G m s if N/2 m N-1 1 s d x (v) Using the modulation property ( ) ( ) f = 2 j x f G f e G 1 1

446 [Example 1] Suppose that for -1 x 1 g(x) = 0 otherwise g(x) = (1-|x|)2 Sampling interval x = 0.1 How do we obtain the FT of g(x) by the DFT? (Solution): ( ) ( ) = 1 g x g x (i) 1 ( ) ( ) g x 1g x x

n ( ) 447 = g g n (ii) 1 d x n ( ) n = G m DFT g (iii) d d blue: real part; red: imaginary part move to left m

448 (iv) Mapping blue: real part; red: imaginary part f change m to f f N = G m G m s if 0 m N/2, 1 d x f N = G m f G m s if N/2 m N-1 1 s d x

449 ( ) ( ) f = 2 j f G f e G (v) 1 blue: real part; red: imaginary part f

450 5.2.3 Transform Pairs and Properties [Duality Property] 1 N = G m DFT g n If = = (( ))N DFT G n Ng m (( ))N g m DFT G n then ((a))N : a modulo N the remainder of a after divided by N = 1,2, if m , 1 g N m if m N = (( )) g m 0 N = 0 g

451 1 N = (( ))N = g m DFT G n If G m DFT g n then 1 N 1 N 1 N 2 mn N j = (Proof): DFT G n e G n = 0 n 1 1 1 1 N N N N 1 N 1 N g k g k 2 2 2 2 mn N nk mn N nk j j j j = = e e e e N N = = = = 0 0 0 0 n k k n 1 N 2 n j a = (( )) a e N Since d N = (proved on the next page) 0 n 1 N 1 N 1 N g k = m k (( )) DFT G n N d N = 0 k 1 N g k = + = (( )) (( )) m k g m d N N = 0 k

452 2 a 1 N j n ((a))N : the remainder of a after divided by N = (( )) a e N N d N = 0 n When a bN where b is some integer 2 a 2 j N a 1 N 1 1 j n 1 e N = = = 0 e N 2 2 a a j j = 0 n 1 1 e e N N When a = bN where b is some integer (i.e., ((a))N= 0) 2 a 1 1 1 N N N j n = = = 2 j bn 1 e e N N = = = 0 0 0 n n n

453 [Determine the IDFT by the DFT] = g n = 1g m DFT G n 1 N g n 1 Note: (i) Computation loading of the IDFT = Computation loading of the DFT (ii) In industry, only the chip of the DFT is required.

454 [Transform Pair] g[n] G[m] (1) d[n] (see page 442) 1 m (2) 1 N d d[n k] exp[ j2 km/N] (3) exp[j2 kn/N] (4) m k N d N N cos[2 kn/N] (5) m k + ( ) m N k 2 2 d d N N ( sin sin[2 kn/N] (6) m k + ( ) j j m N k 2 2 d d ) (7) g[n] = 1 for 0 n W g[n] = 0 otherwise m W + W N sin ( = 1)/ N j m 0, e for m ( ) / m N 0 + 1 W for m (8) exp[-kn], k 0 Nk 1 e e k j m N 2 / 1

455 [Discrete Impulse Train] 1 0 if nis a multiple of c otherwise = c is a factor of N p n c N = 12 2 p n 3p n

456 [Example 2] Determine the DFT of pc[n] c is a factor of N n = ck 2 2 2 N c m m m 1 / 1 / 1 N N c N c j n j ck j k = = e p n e e / N N c = = = 0 0 0 n k k mN 2 N c c 2 N c j m / 1 N c j k 1 e / if m is not a multiple of N/c = = 0 e / 2 N c m j = 0 k 1 e / 2 N c m / 1 N c j k N c = if m is a multiple of N/c e / = 0 k N c m = DFT p n p Therefore, / c N c

457 DFT 6p m 2 p n 6 3p n DFT 4p m 4

458 Properties (1) Linear DFT{ax[n] + by[n]} = aX[m] + bY[m] 1 1 N N (2) DC Values 1 N 0 g n 0 = = , G g G m = = 0 0 n n (3) Shifting = km (( )) DFT g n k W G m N if 0 n N-1 where ((n))N= n ((n))N= n+N if -N n -1 ((n))N= n N if N n 2N-1 ( exp 2 / W j N = DFT W g n = (( )) N DFT g n ) G + k n (4) Modulation (( )) m k N (5) Time Reverse = (( )) G m N (6) Even /Odd Input If g[n] = g[((-n))N], then G[m] = G[((-m))N]; If g[n] = -g[((-n))N], then G[m] = -G[((-m))N];

459 (7) Conjugate DFT{g*[n]} = G*[((-m))N] (8) Real/Imaginary Input If g[n] is real, then G[m] = G*[((-m))N]; If g[n] is pure imaginary, then G[m] = G*[((-m))N]; 1 N h n g k h = (9) Circular Convolution If y n g n = = (( )) n k c N 0 k then Y m y n G m H m c g n 0 k = = = (( )) h n (10) Circular Correlation If N 1 N = + (( )) g k n h k N Y m m = G m H then (11) Parseval s Theorem (Energy Preservation) 1 1 N N g n 2 2 = N G m = = 1 0 0 n m 1 N N (12) Generalized Parseval s Theorem g n h n m = N G m H = = 0 0 n m

460 5.2.4 Discrete Circular Convolution 1 N [Discrete Circular Convolution] g k h h n = n k (( )) g n c N = 0 k (Proof of the convolution property) 2 2 2 2 m m m m 1 1 1 1 N N N N j n j k j s j n 1 N 1 N G m H m e g k e h s e = e N N N N = = = = 0 0 0 0 m m k s 2 ( ) n s k N 1 1 1 N N N j m 1 N g k h s = e = = = 0 1 0 1 0 k N s N m 1 N g k h s N = (( )) n s k d N = = 0 0 k s 1 1 N N g k h = n k (( )) N = = 0 0 k s 2 a 1 N j n ((a))N : the remainder of a after divided by N = (( )) a e N N Here we apply d N = 0 n

461 [Discrete Circular Convolution and Discrete Linear Convolution] A discrete linear time-invariant (LTI) system can always be expressed a discrete linear convolution: 1 0 k = However, the convolution implemented by the DFT is the discrete circular convolution: If 1 N g k h n k = = y n g n h n ( ) y n ( ) = = [ ] [ ] IDFT DFT g n DFT h n IDFT G m H m then 1 N g k h y n h n = = n k (( )) g n c N = 0 k ((a))N : the remainder of a after divided by N

462 1 N h n g k h n k = = y n g n linear convolution: 1 = 0 k N 1 g k h y n h n = = n k circular convolution: (( )) g n c N = 0 k For example, y = 0 h 1 h 2 h 3 4 + + + + + 12 2 1 0 1 2 g g g g h g h 2 0 h + 1 h 2 h 3 4 = + + + + + 2 1 1 3 0 1 2 y g g h g g h N g h N g N The condition where the circular convolution is equal to the linear convolution: 0 0 g n for n or n = 0 0 h n for n or n = (ii) M L (i) + 1 N M L (iii)

463 The condition where the circular convolution is equal to the linear convolution: 0 0 g n for n or n = 0 0 h n for n or n = (ii) M L (i) + 1 N M L (iii) 1 N g k h y n = n k (Proof): (( )) N = 0 k y n 0 2 g n + g N + 0 g h n 2 L h L 1 h n g N 1 + 1 + 1 g N g n h 1 h N + = + + + + + + + 1 1 0 1 g h n g h n + g n h + + + + g n h N h n + + L h L 0 1 h n 1 h n + + 1 1 n g N = + + 0 1 1 g g h n g n h + 1 1 n = + + = 0 g y n 1 (Since N+n+1-L N+1-L M)

464 5.2.5 Complexity 1 N = 2 mn N j G m g n e = 0 n Direct implementation: Complexity = O(N2) With the fast algorithm: Complexity = O(Nlog2N)

465 1 N = 2 mn N j G m g n e = 0 n 2-point DFT [0] [1] 1 1 1 [0] [1] g G G g = 1 g[0] G[0] g[1] G[1] -1

466 When N = 2k 2 mn N 1 N j g n e = G m = 0 n + 2 (2 ) N 2 (2 N 1) m n m n /2 1 /2 1 N N j j = + + 2 2 1 g n e g n e = = 0 0 n n 2 2 2 mn m mn /2 1 /2 1 N N j j j = + g n e e g n e /2 /2 N N N 1 2 = = 0 0 n n twiddle factors g1[n] = g[2n], g2[n] = g[2n+1] Therefore, one N-point DFT = two (N/2)-point DFT + twiddle factors

467 8-point DFT G[0] g[0] G[1] g[2] 4-point DFT G[2] g[4] G[3] g[6] 1 G[4] g[1] -1 w G[5] -1 g[3] 4-point DFT w2 G[6] -1 g[5] w3 G[7] g[7] -1 When N = 8 2 2 j jN = w e 8 = w e w = 4 1 2 2 2 2 N N N + + ( ) j m j m j j m ( ) m = = = = m w e e e e w w = 2 2 N N N N 2 8 1

468 8-point DFT g[0] G[0] g[4] G[1] -1 1 g[2] G[2] -1 w2 g[6] G[3] -1 -1 1 g[1] G[4] -1 w g[5] G[5] -1 -1 1 g[3] G[6] w2 -1 -1 w2 g[7] G[7] w3 -1 -1 -1 2 j = w e 8

469 complexity: log N N 2 Number of stages N/2 twiddle factors are required between stages two J. W. Cooley and J. W. Tukey, An algorithm for the machine computation of complex Fourier series, Mathematics of Computation, vol. 19, pp. 297-301, Apr. 1965. (Cooley-Tukey) C. S. Burrus, Index Mappings for multidimensional formulation of the DFT and convolution, IEEE Trans. Acoustics, Speech, and Signal Processing, vol. 25, pp. 1239-242, June 1977. (Prime factor)