CS 294-73 Software Engineering for Scientific Computing Guest - - PowerPoint PPT Presentation

CS 294-73 Software Engineering for Scientific Computing Guest Lecture: The Discrete Fourier Transform

Doru Thom Popovici

x[n]

• Given a continuous function , its Fourier transform is defined as
• The inverse Fourier transform is defined as

The Discrete Fourier Transform

1

ˆ x[k] =

N−1

X

n=0

x[n]e−i 2π

N kn, ∀

0 ≤ k < N

x[n] = 1 N

N−1

X

k=0

ˆ x[k]ei 2π

N kn, ∀

0 ≤ n < N

The Discrete Fourier Transform – Change of Basis

2

• It maps a signal/function from the time/space domain to the frequency domain

Time domain Frequency domain

ˆ x[k] =

N−1

X

n=0

x[n]e−i 2π

N kn

ˆ x[k]

x[n]

x[n] = 1 N

N−1

X

k=0

ˆ x[k]ei 2π

N kn

The Circular Convolution

12

3

The Circular Convolution

12

4

ˆ x[k]

ˆ h[k]

ˆ y[k]

The Circular Convolution

Precomputed Values/Closed Form Function

12

5

ˆ x[k]

ˆ h[k]

ˆ y[k]

! 𝒊 [ 𝒍 ] c a n b e a n y t h i n g

Not Only Circular Convolution

• Time Shift
• Convolution
• Correlation
• Time Difference
• Differentiation

13

6

y[n] = x[n − n0] ↔ ˆ y[k] = ˆ x[k] · e−i 2π

N nok

y[n] = x[n] ~ h[n] ↔ ˆ y[k] = ˆ x[k] · ˆ h[k]

y[n] = x[n] ? h[n] ↔ ˆ y[k] = ˆ x[k] · ˆ h∗[k]

N−1

X

l=0

αl · y[n − l] =

M−1

X

l=0

βl · x[n − l] ↔ ˆ y[k] = ˆ x[k] · PM−1

l=0

βl · e−i 2π

N lk

PN−1

l=0 αl · e−i 2π

N lk

y(t) = d dtx(t)|t=kn ↔ ˆ y[k] = ˆ x[k] · ˆ h[k]

x[n]

• Given a continuous function , its Fourier transform is defined as
• Important mathematical kernel
• How to efficiently compute it?

The Discrete Fourier Transform

7

ˆ x[k] =

N−1

X

n=0

x[n]e−i 2π

N kn, ∀

0 ≤ k < N

• Apply a DFT of size 16 on the input vector 𝑦

Cooley-Tukey for The DFT of Size 16

8

x

ˆ x[k] =

15

X

n=0

x[n]e−i 2π

16 kn

ˆ e[k] =

7

X

n=0

x[2 · n]e−i 2π

8 kn

• Apply a DFT of size 16 on the input vector 𝑦

Cooley-Tukey for The DFT of Size 16

9

x

DFT8

• Apply a DFT of size 16 on the input vector 𝑦

Cooley-Tukey for The DFT of Size 16

10

x

DFT8

DFT8

ˆ e[k] =

7

X

n=0

x[2 · n]e−i 2π

8 kn

ˆ

• [k] =

7

X

n=0

x[2 · n + 1]e−i 2π

8 kn

ˆ

• 0[k] = e−i 2π

16 k · ˆ

• [k]

• Apply a DFT of size 16 on the input vector 𝑦

Cooley-Tukey for The DFT of Size 16

11

x

DFT8

DFT8

ˆ e[k] =

7

X

n=0

x[2 · n]e−i 2π

8 kn

e−i 2π

16 k·1, ∀

0 ≤ k < 8

• Apply a DFT of size 16 on the input vector 𝑦

Cooley-Tukey for The DFT of Size 16

12

x

DFT8

DFT8

y

e−i 2π

16 k·1, ∀

0 ≤ k < 8

ˆ e[k] + ˆ

• 0[k]

ˆ e[k] − ˆ

• 0[k]

x[n]

• Given a continuous function , its Fourier transform is defined as
• What is the summation?

The Discrete Fourier Transform

13

ˆ x[k] =

N−1

X

n=0

x[n]e−i 2π

N kn, ∀

0 ≤ k < N

x[n]

• Given a continuous function , its Fourier transform is defined as

The Discrete Fourier Transform

14

ˆ x[k] = ⇥ e−i 2π

N 0·k

e−i 2π

N 1·k

. . . e−i 2π

N (N−1)·k⇤

· 2 6 6 6 4 x[0] x[1] . . . x[N − 1] 3 7 7 7 5 , ∀ 0 ≤ k < N

x[n]

• Given a continuous function , its Fourier transform is defined as

The Discrete Fourier Transform

15

     ˆ x[0] ˆ x[1] . . . ˆ x[N − 1]      =      e−i 2π

N 0·0

e−i 2π

N 0·1

e−i 2π

N 0·2

. . . e−i 2π

N 0·(N−1)

e−i 2π

N 1·0

e−i 2π

N 1·1

e−i 2π

N 1·2

. . . e−i 2π

N 1·(N−1)

. . . . . . . . . . . . e−i 2π

N (N−1)·0

e−i 2π

N (N−1)·1

e−i 2π

N (N−1)·2

. . . e−i 2π

N (N−1)·(N−1)

     ·      x[0] x[1] . . . x[N − 1]     

Fourier matrix

x

y

DFTN

x[n]

• Given a continuous function , its Fourier transform is defined as
• It is a matrix vector multiplication, hence it is 𝑃(𝑂*)
• Cooley-Tukey algorithm is a matrix factorization of the Fourier matrix
• A symmetric matrix vector multiplication for prime problem sizes

The Matrix Representation of the Fourier Transform

16

y = DFTN · x

• First steps of the Cooley-Tukey Algorithm

Cooley-Tukey for The DFT of Size 16

17

x

y = DFT16 · x

• Do not treat the data as a flat 1D array

The DFT of Size 16

18

y = DFT16 · x

˜ x

˜ y = DFT2 · ⇣ T 16

8

(DFT8 · ˜ x)T ⌘

• Operations applied on high dimensional objects
• What is the memory footprint of this algorithm?

The DFT of Size 16

19

DFT8

˜ x

transposition pointwise

e−i 2π

N k·1, ∀

0 ≤ k < 8

DFT2

˜ y = DFT4 · ⇣ T 4

4 (DFT4 · ˜

x)T ⌘

The DFT of Size 16 – Different algorithm

20

˜ x

DFT4

transposition pointwise

DFT4

• Operations applied on high dimensional objects
• What is the memory footprint of this algorithm?

Cooley-Tukey – Divide and Conquer Algorithm

21

DFT4

DFT2

DFT32

DFT16

DFT8

DFT64

DFT8

DFT2

DFT4

DFT4

DFT2

DFT64

DFT4

DFT4

DFT4

DFT64

DFT8

DFT2

DFT4

• Multiple algorithms
• Which one performs better?

0.5 1 1.5 2 2.5 3 3.5 4

DFT64 DFT128 DFT256 DFT512 DFT64 batch 4 DFT128 batch 4 DFT256 batch 4 DFT512 batch 4

FFTW Flops/cycle Sizes Performance on Intel Haswell 4770K Custom Approach MKL

Higher is better

Peak Performance

Some Results and Ongoing Research

22

The DFT of Size 5

23

• Given a continuous function , its Fourier transform is defined as
• For prime problem sizes there are Rader and Bluestein algorithms
• Talk about a new approach – cast the computation as a special SYMV
• The DFT matrix has properties and symmetries that can be exploited

x[n]

y = DFT5 · x

Rader and Bluestein Algorithms

24

1 1 1 1 1 𝑦′- 𝑦′. 𝑦′* 𝑦′/ 𝑦′0 𝑧′- 𝑧′. 𝑧′* 𝑧′/ 𝑧′0

temp temp Bluestein Matrix

𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0

input

𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0

• utput

1 1

pre-process post-process

1 1 1 1 1 1 1 1 1 𝑦′- 𝑦′. 𝑦′/ 𝑦′0 𝑦′* 𝑧′- 𝑧′. 𝑧′/ 𝑧′0 𝑧′*

temp temp Rader Matrix

𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0

input

𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0

• utput

permute permute

1 1 1 1 1 1 1 1 1 𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0

input

• utput

𝐸𝐺𝑈5

Rader Algorithm Bluestein Algorithm

𝜕5

• 𝜕5
• 𝜕5
• 𝜕5
• 𝜕5
• 𝜕5
• 𝜕5

.

𝜕5

*

𝜕5

/

𝜕5 𝜕5

• 𝜕5

*

𝜕5 𝜕5

7

𝜕5

8

𝜕5

• 𝜕5

/

𝜕5

7

𝜕5

9

𝜕5

.*

𝜕5

• 𝜕5

𝜕5

8

𝜕5

.*

𝜕5

.7

𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0

𝜕5 = 𝑓<=*>

5

Exploit the Structure of the Matrix without Changing It

25

𝜕5

• 𝜕5
• 𝜕5
• 𝜕5
• 𝜕5
• 𝜕5
• 𝜕5

.

𝜕5

*

𝜕5

/

𝜕5 𝜕5

• 𝜕5

*

𝜕5 𝜕5

7

𝜕5

8

𝜕5

• 𝜕5

/

𝜕5

7

𝜕5

9

𝜕5

.*

𝜕5

• 𝜕5

𝜕5

8

𝜕5

.*

𝜕5

.7

𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0

𝜕5

• = 1

Fourier Matrix – Row and Column of 1s

26

1 1 1 1 1 1 1 1 1 𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0 𝜕5

.

𝜕5

*

𝜕5

/

𝜕5 𝜕5

*

𝜕5 𝜕5

7

𝜕5

8

𝜕5

/

𝜕5

7

𝜕5

9

𝜕5

.*

𝜕5 𝜕5

8

𝜕5

.*

𝜕5

.7

𝜕5

• = 1

Fourier Matrix – Row and Column of 1s

27

1 1 1 1 1 1 1 1 1 𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0 𝜕5

.

𝜕5

*

𝜕5

/

𝜕5 𝜕5

*

𝜕5 𝜕5

7

𝜕5

8

𝜕5

/

𝜕5

7

𝜕5

9

𝜕5

.*

𝜕5 𝜕5

8

𝜕5

.*

𝜕5

.7

𝜕5

@ = 𝜕5 @ % 5

Fourier Matrix – Symmetry and Persymmetry

28

1 1 1 1 1 1 1 1 1 𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0 𝜕5

.

𝜕5

*

𝜕5

/

𝜕5 𝜕5

*

𝜕5 𝜕5

.

𝜕5

/

𝜕5

/

𝜕5

.

𝜕5 𝜕5

*

𝜕5 𝜕5

/

𝜕5

*

𝜕5

.

𝜕5

@ = 𝜕5 @ % 5

Fourier Matrix – Symmetry and Persymmetry

29

1 1 1 1 1 1 1 1 1 𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0 𝜕5

.

𝜕5

*

𝜕5

/

𝜕5

*

𝜕5 𝜕5

.

𝜕5

/

𝜕5

/

𝜕5

.

𝜕5 𝜕5

*

𝜕5 𝜕5

/

𝜕5

*

𝜕5

.

𝜕5

𝜕5

@ = 𝜕5 <(5<@)

Fourier Matrix – Conjugate Symmetry

30

1 1 1 1 1 1 1 1 1 𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0 𝜕5

.

𝜕5

*

𝜕5

<*

𝜕5

*

𝜕5 𝜕5

<0

𝜕5

<*

𝜕5

<*

𝜕5

<0

𝜕5 𝜕5

*

𝜕5

<.

𝜕5

<*

𝜕5

*

𝜕5

.

𝜕5

<.

𝜕5

@ = 𝜕5 <(5<@)

Fourier Matrix – Conjugate Symmetry

31

1 1 1 1 1 1 1 1 1 𝑦- 𝑦. 𝑦* 𝑦/ 𝑦0 𝑧- 𝑧. 𝑧* 𝑧/ 𝑧0 𝜕5

.

𝜕5

*

𝜕5

<*

𝜕5

*

𝜕5 𝜕5

<0

𝜕5

<*

𝜕5

<*

𝜕5

<0

𝜕5 𝜕5

*

𝜕5

<.

𝜕5

<*

𝜕5

*

𝜕5

.

𝜕5

<.

Fourier Matrix – Conjugate Symmetry

32

34 10 20 40 80 160 320 640 1280 2560 5120 10240 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97

Prime-sized DFTs on Intel Haswell 4770K

SymDFT MKL FFTW

Problem Size CPU Cycles [LogScale] Lower is better

Some Results and Ongoing Research

33

• Prefer bottom-up approach – started with codelets that are typically execute on one core
• Extend to larger sizes and higher dimensions
• Exploit thread/core/socket/node parallelism
• Use systematic approaches to optimize things
• Generalize approach to GPU/FPGAs and whatever exotic architecture that is out there

What Next?

34