고속 푸리에 변환
개요1
이산푸리에변환DFT은 수식적인 정의를 순진하게 따라서 계산할 경우 $\mathcal{O}(N^{2})$의 시간복잡도를 가지지만, 아래에서 설명할 알고리즘에 따라 계산하면 시간복잡도가 $\mathcal{O}(N\log_{2}N)$로 내려간다. 이산푸리에변환을 빠르게 수행하는 이러한 알고리즘을 고속푸리에변환fast Fourier transform, FFT이라고 한다.
빌드업
두 숫자를 곱하고, 그것을 다시 다른 숫자와 더하는 것을 $1$번의 연산operation이라고 하자. 그러면 $\sum\limits_{i=0}^{n-1}a_{n}b_{n}$의 값을 얻기 위해선 $n$번의 연산이 필요하다.
$$ \begin{align*} \sum\limits_{n=0}^{0} a_{n}b_{b} &= a_{0}b_{0} = \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \\ \sum\limits_{n=0}^{1} a_{n}b_{b} &= a_{0}b_{0} + a_{1}b_{1} = \overbrace{\big( \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \big) {\color{#5882FA} + } a_{1} {\color{#5882FA} \times} b_{1}}^{\color{#5882FA}2 \text{ operations}} \\ \sum\limits_{n=0}^{2} a_{n}b_{b} &= a_{0}b_{0} + a_{1}b_{1} + a_{2}b_{2} = \overbrace{\bigg( \overbrace{\big( \overbrace{0 {\color{red} +} a_{0} {\color{red} \times} b_{0}}^{\color{red} 1 \text{ operation}} \big) {\color{#5882FA} + } a_{1} {\color{#5882FA} \times} b_{1}}^{\color{#5882FA}2 \text{ operations}} \bigg) {\color{#FE9A2E} + } a_{2} {\color{#FE9A2E} \times} b_{2}}^{\color{#FE9A2E}3 \text{ operations}} \\ \end{align*} $$
이제 이산 푸리에 변환의 정의를 떠올려보자.
선형변환 $\mathcal{F}_{N} : \mathbb{C}^{N} \to \mathbb{C}^{N}$을 이산 푸리에 변환이라 한다.
$$ \mathcal{F}_{N}(\mathbf{a}) = \hat{\mathbf{a}} = \begin{bmatrix} \hat{a}_{0} \\ \hat{a}_{1} \\ \dots \\ \hat{a}_{N-1} \end{bmatrix} ,\quad \hat{a}_{m} = \sum_{n=0}^{N-1}e^{-i2\pi mn /N} a_{n}\quad (0\le m < N) \tag{1} $$
이때 $\mathbf{a} = \begin{bmatrix} a_{0}& a_{1}& \dots& a_{N-1} \end{bmatrix}^{T}$.
$\hat{a}_{m}$을 계산하기 위해서는 $N$번의 연산이 필요하고 $\hat{\mathbf{a}}$를 계산하기 위해선 이를 $N$번 수행해야하므로, 이산 푸리에 변환을 계산하기 위해서는 총 $N^{2}$번의 연산이 필요하다. 즉 $\mathcal{O}(N^{2})$의 시간복잡도를 가진다. 이는 컴퓨터 계산의 측면에서 푸리에 변환이 적지 않은 코스트를 가짐을 의미한다.
알고리즘
데이터의 길이 $N$을 합성수 $N = N_{1}N_{2}$라고 하자. 그리고 인덱스 $m, n$을 다음과 같이 두자.
$$ m = m^{\prime}N_{1} + m^{\prime \prime},\quad n = n^{\prime}N_{2} + n^{\prime \prime} $$
그러면 $0 \le m^{\prime}, n^{\prime \prime} \le N_{2}-1$이고, $0 \le m^{\prime \prime}, n^{\prime} \le N_{1}-1$이다. $(1)$의 지수부분을 다음과 같이 나타내자.
$$ \begin{align*} e^{-i2\pi mn /N} &= e^{-i2\pi (m^{\prime}N_{1} + m^{\prime \prime})(n^{\prime}N_{2} + n^{\prime \prime})/(N_{1}N_{2})} \\ &= e^{-i2\pi (m^{\prime}n^{\prime}N_{1}N_{2} + m^{\prime}n^{\prime \prime}N_{1} + m^{\prime \prime}n^{\prime}N_{2} + m^{\prime \prime}n^{\prime \prime})/(N_{1}N_{2})} \\ &= e^{-i2\pi m^{\prime}n^{\prime}} e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)} \\ &= e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)} \\ \end{align*} $$
이제 이를 $(1)$에 대입하면,
$$ \begin{align*} \hat{a}_{m} &= \sum_{n=0}^{N-1}e^{-i2\pi mn /N} a_{n} \\ &= \sum_{n^{\prime \prime}=0}^{N_{2}-1} \sum_{n^{\prime}=0}^{N_{1}-1} e^{-i2\pi \big( (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime}/N_{1}) + (m^{\prime \prime}n^{\prime \prime}/N) \big)}a_{n^{\prime}N_{2}+n^{\prime \prime}} \\ &= \sum_{n^{\prime \prime}=0}^{N_{2}-1} \left[ \sum_{n^{\prime}=0}^{N_{1}-1} e^{-i2\pi m^{\prime \prime}n^{\prime}/N_{1}}a_{n^{\prime}N_{2}+n^{\prime \prime}} \right] e^{-i2\pi [ (m^{\prime}n^{\prime \prime}/N_{2}) + (m^{\prime \prime}n^{\prime \prime}/N) ] } \end{align*} $$
위 식대로라면 각괄호 안의 $\left[ \sum_{n^{\prime}=0}^{N_{1}-1} \right]$를 계산하는데 $N_{1}$번의 연산, 각괄호 밖의 $\sum_{n^{\prime \prime}=0}^{N_{2}-1}$를 계산하는데 $N_{2}$번의 연산이 필요하므로, $\hat{a}_{m}$을 계산하기 위해서는 총 $(N_{1} + N_{2})$번의 연산이 필요함을 알수 있다. $\hat{\mathbf{a}}$를 얻기 위해선 이를 $N$번 반복해야하므로 총 $N(N_{1} + N_{2})$의 비용이 들어 $N^{2}$보다 줄어들었음을 확인할 수 있다.
각괄호 안을 잘 보면, $N_{1}$이 다시 합성수일 때 같은 논리를 적용할 수 있음을 알 수 있을 것이다. 따라서 일반적으로 데이터의 길이가 $N = N_{1} N_{2} \cdots N_{k}$와 같은 합성수일 때, 시간복잡도는 다음과 같이 줄어든다.
$$ \mathcal{O}(N^{2}) \searrow \mathcal{O}\big( N(N_{1} + N_{2} + \cdots + N_{k}) \big) $$
이제 $N$을 $2$의 제곱수 $N = 2^{k}$라고 해보자. 그러면 $\log_{2}N = k$이고, $N^{2} = 2^{k}$에서 $2^{k}(2k)$만큼 줄어드는 것이므로, 시간 복잡도는 다음과 같이 줄어든다.
$$ \mathcal{O}(N^{2}) \searrow \mathcal{O}(2N \log_{2}N) $$
여담
이는 1965년 Cooley와 Tukey2에 의해 제안되었기 때문에 쿨리-튜키 알고리즘Cooley-Tukey algorithm이라고도 불린다. 하지만 이들이 최초로 발명한 것은 아니다. 가우스도 이와 같은 알고리즘을 연구했지만, 제대로 발표하지 않아 이 사실은 후에 알려지게되었다3.
Gerald B. Folland, Fourier Analysis and Its Applications (1992), p252-253 ↩︎
J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series, Mathematics of Computation 19 (1965), 297-301. ↩︎
M. T. Heideman, D. H. Johnson, and C. S. Burms, Gauss and the history of the fast Fourier transform, Archive for the History of the Exact Sciences 34 (1985), 264-277. ↩︎