Fast Fourier Transform Algorithm

Fast Fourier Transform
Jongwook Choi
2012 SNUPS Seminar
Seoul National University
August 4, 2012
Jongwook Choi Fast Fourier Transform

View Slide

Fourier Transform ?

View Slide

!?!?!?!?!???

View Slide

!?!?!?!?!???
F[x(t)] =
∞
−∞
x(t)e−jωtdt = 2
1
0
(1 − t) cos(ωt)dt
= 2
1
0
cos(ωt)dt −
1
0
t cos(ωt)dt
=
2
ω
sin(ω) −
1
0
t d sin(ωt)
=
2
ω
sin(ω) − t sin(ωt)
1
0
−
1
0
sin(ωt)dt
=
2
ω2
cos(ωt)
0
1
=
2
ω2
(1 − cos(ω))
=
2
ω
2
sin2
ω
2
=
sin2 (πf)
(πf)2
= · · ·

View Slide

Jean Baptiste Joseph Fourier (1768 – 1830)

View Slide

다항식 곱
두 다항식 A(z), B(z) 의 곱 C(z) = A(z)B(z)을 구하고 싶다.
A(z) = a0 + a1z + a2z2 + · · · + an−1zn−1 =
n−1
j=0
ajzj
B(z) = b0 + b1z + b2z2 + · · · + bn−1zn−1 =
n−1
j=0
bjzj

View Slide

다항식 곱
두 다항식 A(z), B(z) 의 곱 C(z) = A(z)B(z)을 구하고 싶다.
A(z) = a0 + a1z + a2z2 + · · · + an−1zn−1 =
n−1
j=0
ajzj
B(z) = b0 + b1z + b2z2 + · · · + bn−1zn−1 =
n−1
j=0
bjzj
=⇒ C(z) = A(z)B(z) = c0 + c1z + · · · + c2n−2z2n−2,
where ci = a0bi + a1bi−1 + · · · + aib0

View Slide

평범한 방법
O(n2)
C(z) = A(z)B(z) = c0 + c1z + · · · + c2n−2x2n−2,

View Slide

평범한 방법
O(n2)
C(z) = A(z)B(z) = c0 + c1z + · · · + c2n−2x2n−2,
좀 더 빠르게 하고 싶습니다.
FFT! (Fast Fourier Transform)

View Slide

신묘한 방법
Idea. m − 1차 다항식은 m개의 점에서의 함수값이 결정되면
유일하게 결정된다 (Lagrange’s Interpolation).

View Slide

신묘한 방법
C(z) 의 degree 는 (at most) 2n − 2 이므로,
서로 다른 2n − 1 개의 점에 대해 C 의 함수값만 알면 되겠다.

View Slide

신묘한 방법
C(z) 의 degree 는 (at most) 2n − 2 이므로,
서로 다른 2n − 1 개의 점에 대해 C 의 함수값만 알면 되겠다.
서로 다른 2n − 1 개의 점 z0, z1, · · · , z2n−2
에 대해,
A(z0), A(z1), · · · , A(z2n−2) 와
B(z0), B(z1), · · · , B(z2n−2) 를 안다면,
C(z0), C(z1), · · · , C(z2n−2) 를 빠르게 계산할 수 있지 않을까?

View Slide

Multiplication by FFT
서로 다른 2n − 1 개의 점 z0, z1, · · · , z2n−2
에 대해,
A(z0), A(z1), · · · , A(z2n−2) 와
B(z0), B(z1), · · · , B(z2n−2) 를 안다면,

View Slide

Multiplication by FFT
서로 다른 2n − 1 개의 점 z0, z1, · · · , z2n−2
에 대해,
A(z0), A(z1), · · · , A(z2n−2) 와
B(z0), B(z1), · · · , B(z2n−2) 를 안다면,
C(zk) = A(zk) · B(zk) k = 0, 1, · · · , 2n − 2

View Slide

Discrete Fourier Transform
다항식 A(x) = a0 + a1x + · · · + an−1xn−1 을 생각하자.
ω ∈ C 을 primitive root of unity (ωn = 1 라는 말) 라고 할 때,
n개의 점 ω0, ω1, ω2, · · · , ωn−1 에서의 함수값
A(ωi) =
n−1
j=0
aj(ωi)j
들을 모으면 n 개의 수가 얻어진다.

View Slide

Discrete Fourier Transform
다항식 A(x) = a0 + a1x + · · · + an−1xn−1 을 생각하자.
ω ∈ C 을 primitive root of unity (ωn = 1 라는 말) 라고 할 때,
n개의 점 ω0, ω1, ω2, · · · , ωn−1 에서의 함수값
A(ωi) =
n−1
j=0
aj(ωi)j
들을 모으면 n 개의 수가 얻어진다. 이 때 transform Cn → Cn
(a0, a1, · · · , an−1) → (A(a0), A(a1), · · · , A(an−1))
를 DFT (이산 푸리에 변환) 이라고 한다.

View Slide

ω : primitive n-th root of unity
ωn = 1

View Slide

ω : primitive n-th root of unity
ωn = 1
ω = e2πi
n
= cos(2π/n) + i sin(2π/n)

View Slide

Multiplication by DFT
(a0, a1, · · · , an−1) →


n−1
j=0
aj(ω0)j, · · · ,
n−1
j=0
aj(ωn−1)j


(a0, · · · , an−1
을 계수로하는 다항식에 ω0, · · · , ωn−1을 넣음)

View Slide

(a0, a1, · · · , an−1) →


n−1
j=0
aj(ω0)j, · · · ,
n−1
j=0
aj(ωn−1)j


(a0, · · · , an−1
그러면, C(z) = A(z)B(z)일 때,
A(ωj) · B(ωj) = C(ωj)
인가?

View Slide

(a0, a1, · · · , an−1) →


n−1
j=0
aj(ω0)j, · · · ,
n−1
j=0
aj(ωn−1)j


(a0, · · · , an−1
그러면, C(z) = A(z)B(z)일 때,
A(ωj) · B(ωj) = C(ωj)
인가?
· · · · · · · · · 너무 당연하다.

View Slide

그래서 우리의 관심은
DFT 구하는 것, 즉 ω0, · · · , ωn−1 에서의 함수값
A(ω0), · · · , A(ωn−1) 을 모두 알아낼 수 있으면 OK.
얼마나 빨리?

View Slide

그래서 우리의 관심은
DFT 구하는 것, 즉 ω0, · · · , ωn−1 에서의 함수값
A(ω0), · · · , A(ωn−1) 을 모두 알아낼 수 있으면 OK.
얼마나 빨리? O(n log n)
naive 하게 계산하면 O(n2) 인데 그러면 하나마나.

View Slide

Divide-And-Conquer
다항식 A(x)가 있을 때, A(ω0), · · · , A(ωn−1)을 구하고 싶다.
A(x) = a0 + a1x + a2x2 + · · · + an−1xn−1

View Slide

Divide-And-Conquer
A(x) = a0 + a1x + a2x2 + · · · + an−1xn−1
= (a0 + a2x2 + a4x4 + · · · ) + (a1x + a3x3 + a5x5 + · · · )
= (a0 + a2x2 + a4x4 + · · · ) + x(a1 + a3x2 + a5x4 + · · · )

View Slide

Divide-And-Conquer
A(x) = a0 + a1x + a2x2 + · · · + an−1xn−1
= (a0 + a2x2 + a4x4 + · · · ) + (a1x + a3x3 + a5x5 + · · · )
= (a0 + a2x2 + a4x4 + · · · ) + x(a1 + a3x2 + a5x4 + · · · )
= Aeven(x2) + xAodd(x2)
n이 짝수면..
Aeven(x) = a0 + a2x + a4x2 + · · · + an−2xn/2−1
Aodd(x) = a1 + a3x + a5x2 + · · · + an−1xn/2−1

View Slide

Divide-And-Conquer
A(x) = Aeven(x2) + xAodd(x2)

View Slide

Divide-And-Conquer
1 n/2개의 점 (ω0, ω2, · · · , ωn−2) 에서 두 다항식 Aeven, Aodd
에서의 함수값을 각각 구한다 (즉, DFT 한다).
원래와 동일한 모양의 문제: ω를 ω2로 두고 재귀ㄱㄱㄱ1
Aeven
(ω0), Aeven
(ω2), · · · , Aeven
(ωn−2)
Aodd
(ω0), Aodd
(ω2), · · · , Aodd
(ωn−2)
1n/2 차 다항식에 대하여 (ω2)n/2 = 1 이므로 ω2가 primitive root of unity.

View Slide

Divide-And-Conquer
2 그러면 O(n) time 에 n개의 함수값 A(ωj)를 모조리 구할 수
있다.1
A(ωj) = Aeven(ω2j) + ωj · Aodd(ω2j)
A(ωj+n/2) = Aeven(ω2j+n) + ωjωn/2 · Aodd(ω2j+n)
= Aeven(ω2j) − ωj · Aodd(ω2j)
for 0 ≤ j < n/2.
1Note that ωn = 1, ωn/2 = −1.

View Slide

Divide-And-Conquer
2 그러면 O(n) time 에 n개의 함수값 A(ωj)를 모조리 구할 수
있다.
T(n) = 2T(n/2) + O(n) =⇒ T(n) = O(n log n)

View Slide

역변환
거꾸로, DFT의 결과를 알고 있을 때 원래 다항식을 구하는
방법은? 즉, 최대 2n차 다항식 C에 대해
C(ω0), C(ω1), · · · , C(ω2n−1)
를 알고 있을 때 (ω2n = 1) C의 계수 c0, c1, · · · , c2n−1
를 구하려면
어떻게 하는가?

View Slide

DFT 역변환
n − 1차 다항식 A(x) 에 대해 wn = 1이고, yk = A(wk) 를 모두
알고 있다고 하자. (k = 0, 1, 2, · · · , n − 1)
yk = A(wk) =
n−1
j=0
aj(wk)j

View Slide

DFT 역변환
n − 1차 다항식 A(x) 에 대해 wn = 1이고, yk = A(wk) 를 모두
알고 있다고 하자. (k = 0, 1, 2, · · · , n − 1)
yk = A(wk) =
n−1
j=0
aj(wk)j







y0
y1
y2
.
.
.
yn−1







=







1 1 1 · · · 1
1 ω ω2 · · · ωn−1
1 ω2 ω4 · · · ω2(n−1)
.
.
.
.
.
.
.
.
.
...
.
.
.
1 ωn−1 ω2(n−1) · · · ω(n−1)(n−1)














a0
a1
a2
.
.
.
an−1








View Slide

DFT 역변환
역행렬을 구해서 풀면..







y0
y1
y2
.
.
.
yn−1







=







1 1 1 · · · 1
1 ω ω2 · · · ωn−1
1 ω2 ω4 · · · ω2(n−1)
.
.
.
.
.
.
.
.
.
...
.
.
.
1 ωn−1 ω2(n−1) · · · ω(n−1)(n−1)














a0
a1
a2
.
.
.
an−1














a0
a1
a2
.
.
.
an−1







=
1
n







1 1 1 · · · 1
1 ω−1 ω−2 · · · ω−(n−1)
1 ω−2 ω−4 · · · ω−2(n−1)
.
.
.
.
.
.
.
.
.
...
.
.
.
1 ω−(n−1) ω−2(n−1) · · · ω−(n−1)(n−1)














y0
y1
y2
.
.
.
yn−1








View Slide

DFT 역변환
즉, (y0, y1, · · · , yn−1) 를 DFT 하는데

View Slide

DFT 역변환
즉, (y0, y1, · · · , yn−1) 를 DFT 하는데
ω를 ω−1 로 replace 하고
최종 결과를 n 으로 나누면

View Slide

DFT 역변환
즉, (y0, y1, · · · , yn−1) 를 DFT 하는데
ω를 ω−1 로 replace 하고
최종 결과를 n 으로 나누면
(a0, a1, a2, · · · , an−1) 이 얻어진다.
즉 FFT 함수가 있으면 역변환 하는것도 바로 된다. O(n log n)

View Slide

FFT로 O(n log n) 곱셈하는 알고리즘 정리
A(z) = a0 +a1z+· · ·+an−1zn−1, B(z) = b0 +b1z+· · ·+bn−1zn−1
(image from Introduction to Algorithms; Corment et al.)

View Slide

FFT로 O(n log n) 곱셈하는 알고리즘 정리
(Cooley–Tukey FFT algorithm)
A(z) = a0 +a1z+· · ·+an−1zn−1, B(z) = b0 +b1z+· · ·+bn−1zn−1
1 A, B를 2n − 1차 다항식으로 생각해서,
(a0, a1, · · · , a2n−1)의 DFT (p0, p1, · · · , p2n−1)을 구하고
(b0, b1, · · · , b2n−1)의 DFT (q0, q1, · · · , q2n−1)을 구한다.
FFT 두번 하므로 O(n log n) 걸린다.
2 구한 두 DFT 벡터의 pointwise 곱 (p0q0, · · · , p2n−1q2n−1)을
구한다. O(n)
3 이 녀석을 다시 DFT 하면, (2n · c0, 2n · c1, · · · , 2n · c2n−1) 이
얻어진다. O(n log n)
4 c0, c1, · · · , c2n−1
이 C(z) = A(z)B(z)의 계수이다.

View Slide

Remark
DFT 를 구할 때, FFT 계산 시 n을 2의 power라고 가정하고
있음.
실수 연산을 해야하기 때문에 실수 오차가 발생할 수 있다.
(input 다항식의 계수가 정수라고 가정하면) 최종적으로
얻어진 cj
계수들이 정수가 아닐 수도 있는데, 가장 가까운
integer로 반올림하면 된다. (오차 0.5 미만이면 OK)
실수 연산 없이, ﬁnite ﬁeld Zp (n | p − 1) 에서 ω의 역할을
하는 primitive root of unity 를 쓰면 정수 연산만으로도 되는데
(Number Theoretic Transform), 뭔가 엄청 복잡하고 실제로도
나눗셈 연산때문에 더 느림. (CRLS Problem 30-6 참고)
FFT는 팀노트에 넣어가면 좋다.

View Slide

예제 문제.
두 n-length array A = (a0, · · · , an−1), B = (b0, · · · , bn−1) 가
있다. 이 때 A와 B의 convolution C = A ∗ B = (c0, · · · , cn−1)
cj =
j
k=0
akbj−k
를 구하시오. 단, n은 무지무지 큼.

View Slide

예제 문제.
있다. 이 때 A와 B의 convolution C = A ∗ B = (c0, · · · , cn−1)
cj =
j
k=0
akbj−k
앞에서 배운 Fast Polynomial Multiplication 과 거의 똑같이 하면
된다(!). Exercise.

View Slide

예제 문제 (SRM 518 Div1 Hard)
있다. 이 때 A와 B의 xor-convolution C = A ⊕ B = (c0, · · · , cn−1)
ck =
i xor j=k
aibj

View Slide

O(n log n) 에 도전. FFT의 아이디어를 이용해서 풉니다.

View Slide

어떤 DFT F : Zn → Zn 를 찾자.
F(A), F(B) 를 구한다.
F(A), F(B) 를 pointwise 로 곱해서 F(A ⊕ B) 를 구한다.
F(A ⊕ B) 에 F의 역변환을 끼얹어서 A ⊕ B 를 구한다.

View Slide

어떤 DFT F : Zn → Zn 를 찾자.
F(A), F(B) 를 구한다.
F(A), F(B) 를 pointwise 로 곱해서 F(A ⊕ B) 를 구한다.
F(A ⊕ B) 에 F의 역변환을 끼얹어서 A ⊕ B 를 구한다.
그런 F 를 어떻게 찾지?

View Slide

통밥의 정석
n = 2일 때를 가정해 보자.
A = (a, b), B = (c, d) 라고 하면 A ⊕ B = (ac + bd, ad + bc) 이다.

View Slide

통밥의 정석
F(a, b) · F(c, d) = F(ac + bd, ad + bc) 가 되어야 한다.

View Slide

통밥의 정석
F(a, b) = (a − b, a + b) 라고 가정해보자(???).

View Slide

통밥의 정석
F(a, b) = (a − b, a + b) 라고 가정해보자(???).
F(a, b) · F(c, d) = (a − b, a + b) · (c − d, c + d)

View Slide

통밥의 정석
F(a, b) = (a − b, a + b) 라고 가정해보자(???).
F(a, b) · F(c, d) = (a − b, a + b) · (c − d, c + d)
= ((a − b)(c − d), (a + b)(c + d))

View Slide

통밥의 정석
F(a, b) = (a − b, a + b) 라고 가정해보자(???).
F(a, b) · F(c, d) = (a − b, a + b) · (c − d, c + d)
= ((a − b)(c − d), (a + b)(c + d))
= (ac − ad − bc + bd, ac + ad + bc + bd)

View Slide

통밥의 정석
F(a, b) = (a − b, a + b) 라고 가정해보자(???).
F(a, b) · F(c, d) = (a − b, a + b) · (c − d, c + d)
= ((a − b)(c − d), (a + b)(c + d))
= ((ac + bd) − (ad + bc), (ac + bd) + (ad + bc))

View Slide

통밥의 정석
F(a, b) = (a − b, a + b) 라고 가정해보자(???).
F(a, b) · F(c, d) = (a − b, a + b) · (c − d, c + d)
= ((a − b)(c − d), (a + b)(c + d))
= ((ac + bd) − (ad + bc), (ac + bd) + (ad + bc))
= F(ac + bd, ad + bc)
Olleh!

View Slide

통밥의 정석
n = 4일 때에는? 좀 복잡한 듯.

View Slide

통밥의 정석
n = 2일 때 F(a, b) = (a − b, a + b) 였으니까,

View Slide

통밥의 정석
n = 2일 때 F(a, b) = (a − b, a + b) 였으니까, 느낌상
F(X Y ) = (F(X) − F(Y ) F(X) + F(Y ))
일 것 같다. (???) 세상은 원래 이렇게 아름다운 거에요

View Slide

통밥의 정석
n = 2일 때 F(a, b) = (a − b, a + b) 였으니까, 느낌상
F(X Y ) = (F(X) − F(Y ) F(X) + F(Y ))
일 것 같다. (???) 세상은 원래 이렇게 아름다운 거에요
사실 저 추측만 하면, 증명은 간단한 수학적 귀납법으로 어렵지
않다 (자세한 것은 editorial 참고).

View Slide

역변환
F(X Y ) 가 주어질 때, X Y 를 찾을 수 있는가? 즉, X, Y 를
찾을 수 있는가?
F(X Y ) = (F(X) − F(Y ) F(X) + F(Y )) 를 알고
있으므로, F(X) 와 F(Y ) 를 알 수 있다.
길이가 반으로 줄었으므로 재귀적으로 풀어서 X, Y 를 얻는다.
X, Y 를 concatenate 시키면 끝.

View Slide

FFT: 기타 응용
Fast Integer/Polynomial Multiplication
Pattern Matching (Text Algorithms)
Signal Processing
Image Processing
· · ·

View Slide

The End
Thank You Very Much
Any Questions?

View Slide

Fast Fourier Transform Algorithm

Fast Fourier Transform Algorithm

Jongwook Choi

More Decks by Jongwook Choi

Other Decks in Programming

Featured

Transcript