Fast Fourier Transform (FFT)
N-point DFT
N-point DFT와 역변환 IDFT는 아래와 같이 정의된다.
$ X[k] \triangleq \sum_{n=0}^{N-1}x[n]e^{-i2\pi nk \over N} \quad (k=0,...,N-1) $
$ x[n] \triangleq {1 \over N} \sum_{k=0}^{N-1}X[k]e^{i2\pi nk \over N} \quad (n=0,...,N-1) $
정의에 따라 모든 frequency bin $k [0,N-1]$에 대해 DFT $X[k]$를 계산하는데 $O(N^2)$의 복잡도를 가진다.
Radix-2 DIT FFT
Radix-2 FFT는 N-point DFT를 두 개의 N/2-Point DFT로 분해할 수 있다는 성질을 활용하여 DFT를 낮은 복잡도로 계산하는 알고리즘이다. N-point DFT의 복잡도는 $O(N^2)$이기 때문에 N-Point DFT를 계산하는 것보다 분해된 두 개의 N/2-Point DFT를 계산하는게 복잡도가 낮다는 idea로 부터 시작한다.
앞의 DFT의 정의에서 $e^{-i2 \pi \over N}$을 표기상의 편의를 위해 $W_N$으로 정의하자.1
$ X[k] = \sum_{n=0}^{N-1}x[n]W_N^{nk} \quad (k=0,...,N-1) $
식을 n에 대한 홀수항($2n+1$)과 짝수항($2n$)으로 나누면 아래와 같이 변경할 수 있다.
$ X[k] = \sum_{n=0}^{{N \over 2}-1} x[2n]W_N^{2nk} + \sum_{n=0}^{{N \over 2}-1} x[2n+1]W_N^{(2n+1)k}$
위 식으로 부터 아래와 같이 N-point DFT를 두 개의 N/2-point DFT로 분해할 수 있다.
$$ \begin{align} X[k] &= \sum_{n=0}^{{N \over 2}-1} x[2n]W_N^{2nk} + \sum_{n=0}^{{N \over 2}-1} x[2n+1]W_N^{(2n+1)k}\\ &= \underbrace{\sum_{n=0}^{{N \over 2}-1} x[2n]W_N^{2nk}}_{X_{even}[k]} + W_N^k\underbrace{\sum_{n=0}^{{N \over 2}-1} x[2n+1]W_N^{2nk}}_{X_{odd}[k]} \quad (k=0,...,{N\over2}-1)\\ X[k+{N\over2}] &= \sum_{n=0}^{{N \over 2}-1} x[2n]W_N^{2n(k+{N\over2})} + W_N^{k+{N\over2}}\sum_{n=0}^{{N \over 2}-1} x[2n+1]W_N^{2n(k+{N\over2})}\quad (W_N^{N\over2} = e^{-i\pi}= -1, W_N^{N} = e^{-i2\pi} = 1)\\ &= \underbrace{\sum_{n=0}^{{N\over2}-1}x[2n]W_N^{2nk}}_{X_{even}[k]}-W_N^k\underbrace{\sum_{n=0}^{{N\over2}-1}x[2n+1]W_N^{2nk}}_{X_{odd}[k]} \quad (k=0,...,{N\over2}-1) \end{align} $$
이러한 방법을 다시 재귀적으로 적용할 수 있다. N/2-Point DFT는 N/4-Point DFT로부터 구할 수 있으며, ..., 4-Point DFT는 2-Point DFT로부터 구할 수 있다.
두 개의 N/2-Point DFT을 합쳐서 N-Point DFT를 만드는데 O(N)이 걸리므로 다음과 같은 점화식을 세워 시간복잡도를 계산할 수 있다.
$$\begin{align}
T(N) &= 2T({N\over2})+O(N)\\
&= O(N\log_2N)
\end{align}$$
앞에서의 분해 과정을 이용하여 DFT를 계산하는 알고리즘을 Radix-2 DIT(decimation in time) FFT라고 부른다.2
여기서 decimation이 어떤 의미일까? decimation은 로마군의 형벌로, 제비뽑기로 병사 10명 중에 1명을 처형하는 형벌을 의미한다고 한다. 신호처리에서는 일부 sample을 버려서 downsampling을 한다는 의미로 쓰인다. 위 그림에서는 트리의 하위로 갈 때 마다 매 2개의 sample마다 한개의 sample을 버리는 decimation을 한다.
이러한 의미로 보면 Radix-2 FFT는 서로 다른 shift가 적용된 downsampling된 DFT들로부터 원래의 DFT를 재귀적으로 합성해나가는 알고리즘으로 생각할 수 있다.
구현 1 (recursive)
#include <complex>
#include <vector>
using namespace std;
void FFT_recursive(vector<complex<double>> &x) {
const double pi = 3.14159265358979323846;
int N = x.size(); // assume n is a power of 2
if (N == 1) return;
vector<complex<double>> even(N/2), odd(N/2);
for (int i = 0; i < N; i++) {
if (i%2 == 1) odd[i/2] = x[i];
else even[i/2] = x[i];
}
FFT_recursive(even);
FFT_recursive(odd);
for (int k = 0; k < N/2; k++) {
complex<double> Wk = polar(1.0, -2.0 * pi * k / N);
x[k] = even[k] + Wk*odd[k];
x[k+N/2] = even[k] - Wk*odd[k];
}
}앞의 FFT 수식을 그대로 코드로 옮기면 된다.
이 구현은 직관적인 장점이 있지만 복사 과정으로 인해 실제 성능은 좋지 못하다.
구현 2 (in-place)
앞의 구현에서 문제였던 복사과정을 없앨 수 있는 방법이 있다. 처음부터 입력 데이터를 순서를 바꾸어 결과값이 필요한 위치에 미리 배치해 놓는다면 복사 과정을 없앨 수 있다.
위의 그림을 예로 들면 입력 데이터의 순서를 리프 노드의 인덱스(0,8,4,12...,7,15)대로 재배열을 해주면 된다.
리프노드의 인덱스는 비트 순서를 뒤집었을 때의 값을 기준으로 오름차순으로 정렬한 것인데, 비트를 뒤집은 값에 해당하는 인덱스에 배치된다는 규칙을 찾을 수 있다. 예를 들어 3번 인덱스의 값은 0011을 뒤집은 1100, 12번 인덱스에 배치하면 된다.
이러한 인덱스의 순열을 bit-reversal permutation이라고 한다. 구현시 bit reversal은 대칭적(0011은 1100으로, 1100은 0011로)이므로 절반에 대해서만 swap을 수행하면 된다.
#include <vector>
#include <complex>
using namespace std;
void FFT_inplace(vector<complex<double>> &x, bool inverse) {
const double pi = 3.14159265358979323846;
int N = x.size(); // assume n is a power of 2
// bit reversal permutation
for (int i = 1, rev = 0; i < N; i++) {
int bit = (N >> 1);
while (rev >= bit) rev -= bit, bit >>= 1;
rev += bit;
// bit reversal is symmetric
if (i < rev) swap(x[i], x[rev]);
}
// iterative implementation
for (int i = 2; i <= N; i <<= 1) {
for (int j = 0; j < N; j += i)
for (int k = 0; k < i / 2; k++) {
complex<double> Wk = polar(1.0, (inverse ? -1 : 1) * -2.0 * pi * k / i);
complex<double> even = x[j + k], odd = x[j + k + i / 2];
x[j + k] = even + Wk * odd;
x[j + k + i / 2] = even - Wk * odd;
}
}
// IFFT scaling
if (inverse) for (int i = 0; i < N; i++) x[i] /= N;
}
참고할만한 내용
DIF(decimation in frequency)
반대로 frequency domain에서 위의 분해 과정을 진행할 수 있는데 이를 DIF(decimation in frequency)라고 부른다. 구현상으로는 DIT는 입력을 홀수항과 짝수항으로 절반을 나누고, DIF는 입력을 앞 절반의 항과 뒤 절반의 항으로 나누는 차이가 있다. DIT는 입력이 bit reversal이 되며, DIF는 출력이 bit reversal이 된다.