# About
Fast Fourier Transform is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in the frequency domain and vice versa. The DFT is obtained by decomposing a sequence of values into components of different frequencies.
(reference)
# Concept
이산 푸리에 변환(DFT)을 고속으로 처리하는 알고리즘으로 다항식의 연산, 큰 수의 곱연산의 고속화에도 활용될 수 있다.
먼저 고속 푸리에 변환이 무엇을 하는 것인지 직관적으로 이해해보자.
다항식의 표현방식에는 coefficient representation과 point-value representation이 있다.
coefficient representation은 단어 그대로 계수를 벡터로 나열함으로써 다항식을 표현하는 우리에게 가장 직관적인 방법이고 point-value representatoin은 다항식이 지나는 서로 다른 점들로 표현하는 방식이다.
다항식 $a_{0} + a_{1}x^1 + a_{2}x^2 + \ldots + a_{n}x^n$ 에 대해
$FFT(a_{0} + a_{1}x^1 + a_{2}x^2 + \ldots + a_{n}x^n)$ $\Rightarrow$ $((x_{0}, f(x_{0}), (x_{1}, f(x_{1}), \ldots , (x_{n}, f(x_{n}))$
위와 같이 point-value represetation으로 변환하는 역할이라고 생각하면 직관적으로 이해하기 편하다.
# Core Idea
$f(x) = a_{0} + a_{1}x^1 + a_{2}x^2 + \ldots + a_{n-1}x^{n-1}$ 에 대한 point-value represetation을 위해서 서로 다른 점 $n$개가 필요하며 이는 서로 다른 x좌표를 $n$개 잡으면 된다.
다항식을 짝수 차수 항들과 홀수 차수 항들로 나누어 정리해보자.
$f(x) = (a_{0} + a_{2}x^2 + \ldots ) + (a_{1}x + a_{3}x^3 + \ldots)$
$=(a_{0} + a_{2}x^2 + \ldots ) + x(a_{1} + a_{3}x^2 + \ldots)$
$=P_{even}(x^2) + xP_{odd}(x^2)$
묶인 짝수 차수 항들과 홀수 차수 항들에 붙은 $x$의 차수가 동일한 규칙을 가진 다는 것을 눈치 챌 수 있다.
이렇게 동일한 패턴의 쌍들에 대해서 계산량을 줄이는 것이 FFT의 첫번째 중요한 포인트이다.
두번째 포인트는 $P_{even}(x^2) + xP_{odd}(x^2)$에 대해서
$f(-x) = P_{even}((-x)^2) + (-x)P_{odd}((-x)^2)$
$= P_{even}(x^2) - xP_{odd}(x^2)$
가 성립한다는 점이다. $f(x)$와 $f(-x)$를 구하기 위해서는 $P_{even}(x^2)$ 와 $xP_{odd}(x^2)$의 값을 알면 충분하고 유일한 차이는 더하느냐 빼느냐이다.
여기서 이제 위의 두가지 포인트를 묶어버리는 정말 핵심적인 아이디어가 나온다.
$(x_{0}, x_{1}, \ldots , x_{n-1}) \gets$ $((e^{\frac{2\pi i}{n}})^0, (e^{\frac{2\pi i}{n}})^1, \ldots )$
$w =e^{\frac{2\pi i}{n}}$라 하자.
$w^(n) = w^0 = 1$ 이며 $\forall w \in W_{n},$ $-w \in W_{n}$ 이다.
또한 $\mid w \mid = 1$, ${w^0, w^1, \ldots , w^{n-1}}$는 곱셈에 대하여 닫혀있다(자명히 제곱 연산에 대해서도 닫혀있다).
이때 $n = 2^k$꼴을 만족한다면 $w^n = 1$로부터 $w^{\frac{n}{2}} = -1$임을 알 수 있고
$f(w^{i+{n \over 2}}) = f(-w^i)$ $= P_{even}((-w^i)^2) + (-w^i)P_{odd}((-w^i)^2)$ $= P_{even}(w^2i) - w^i P_{odd}(w^2i)$
$f(w^i) = P_{even}(w^2i) + w^iP_{odd}(w^2i)$
으로부터 계산량을 줄일 수 있게 되는 것이다.
# Algorithm
가장 일반적이고 유명한것은 Cooley–Tukey algorithm이다.
쿨리-툴키 알고리즘를 구현할때 radix-2의 형태를 가장 많이 쓴다. 이를 소개해보도록 하겠다.
알고리즘의 원리는 위의 # Core Idea에서 설명한 포인트들과 동일하기에 생략하겠다.
구현시, radix-2의 형태로 맞춰주기 위하여(위의 $n = 2^k$ 조건을 떠올려도 된다)
$n$차 다항식이라면 $n+1$보다 크거나 같은 가장 작은 2의 거듭제곱을 벡터의 길이로 잡는다.
1. 재귀적 구현
$f(w^{i+{n \over 2}}) = f(-w^i)$ $= P_{even}((-w^i)^2) + (-w^i)P_{odd}((-w^i)^2)$ $= P_{even}(w^2i) - w^i P_{odd}(w^2i)$
$f(w^i) = P_{even}(w^2i) + w^i P_{odd}(w^2i)$
을 이용하여 길이가 1이 될때까지 분할정복을 진행하면 된다.
void fft(vector<complex>& X, cpx w) {
int len = X.size();
if(len == 1) return;
vector<complex> even(len >> 1), odd(len >> 1);
for(int i = 0; i < len; i++) {
if(i & 1) odd[i >> 1] = X[i];
else even[i >> 1] = X[i];
}
fft(even, w*w);
fft(odd, w*w);
//ceof = 1 + 0i;
complex coef(1, 0);
for(int i = 0; i < len/2; i++) {
X[i] = even[i] + coef * odd[i];
X[i + len/2] = even[i] - coef * odd[i];
coef *= w;
}
}
2. 비재귀적 구현
같은 논리라도 비재귀적인 구현으로 바꾸면 비교적으로 복잡해지는 경우가 많다. 이도 그중 하나이다.
Radix-2 bit traversal 방식은 우리가 원하는 인덱스의 순서와 일치하기 때문에 이를 이용한다.
void fft(vector<complex> &a) {
int n = a.size();
vector<complex> w(n/2), aux(n);
for(int i=0; i<n/2; i++){
int k = i&-i;
if(i == k){
double ang = 2 * M_PI * i / n;
// w[i] = cos(ang) + i*sin(ans)
w[i] = complex(cos(ang), sin(ang));
}
else w[i] = w[i-k] * w[k];
}
for(int i=n/2; i; i>>=1){
aux = a;
for(int k=0; 2*k<n; k+=i){
for(int j=0; j<i; j++){
cpx u = aux[2*k + j], v = aux[2*k + j + i] * w[k];
a[k + j] = u + v;
a[k + j + n/2] = u - v;
}
}
}
}# Time Complexity
$T(N) = 2T(\frac{N}{2}) + O(N)$
으로부터 $O(Nlog(N))$ 이 성립한다.