logo

고속 푸리에 변환 📂푸리에해석

고속 푸리에 변환

개요1

이산푸리에변환DFT은 수식적인 정의를 순진하게 따라서 계산할 경우 O(N2)\mathcal{O}(N^{2})의 시간복잡도를 가지지만, 아래에서 설명할 알고리즘에 따라 계산하면 시간복잡도가 O(Nlog2N)\mathcal{O}(N\log_{2}N)로 내려간다. 이산푸리에변환을 빠르게 수행하는 이러한 알고리즘을 고속푸리에변환fast Fourier transform, FFT이라고 한다.

빌드업

두 숫자를 곱하고, 그것을 다시 다른 숫자와 더하는 것을 11번의 연산operation이라고 하자. 그러면 i=0n1anbn\sum\limits_{i=0}^{n-1}a_{n}b_{n}의 값을 얻기 위해선 nn번의 연산이 필요하다.

n=00anbb=a0b0=0+a0×b01 operationn=01anbb=a0b0+a1b1=(0+a0×b01 operation)+a1×b12 operationsn=02anbb=a0b0+a1b1+a2b2=((0+a0×b01 operation)+a1×b12 operations)+a2×b23 operations \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*}

이제 이산 푸리에 변환의 정의를 떠올려보자.

선형변환 FN:CNCN\mathcal{F}_{N} : \mathbb{C}^{N} \to \mathbb{C}^{N}이산 푸리에 변환이라 한다.

FN(a)=a^=[a^0a^1a^N1],a^m=n=0N1ei2πmn/Nan(0m<N)(1) \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}

이때 a=[a0a1aN1]T\mathbf{a} = \begin{bmatrix} a_{0}& a_{1}& \dots& a_{N-1} \end{bmatrix}^{T}.

a^m\hat{a}_{m}을 계산하기 위해서는 NN번의 연산이 필요하고 a^\hat{\mathbf{a}}를 계산하기 위해선 이를 NN번 수행해야하므로, 이산 푸리에 변환을 계산하기 위해서는 총 N2N^{2}번의 연산이 필요하다. 즉 O(N2)\mathcal{O}(N^{2})시간복잡도를 가진다. 이는 컴퓨터 계산의 측면에서 푸리에 변환이 적지 않은 코스트를 가짐을 의미한다.

알고리즘

데이터의 길이 NN합성수 N=N1N2N = N_{1}N_{2}라고 하자. 그리고 인덱스 m,nm, n을 다음과 같이 두자.

m=mN1+m,n=nN2+n m = m^{\prime}N_{1} + m^{\prime \prime},\quad n = n^{\prime}N_{2} + n^{\prime \prime}

그러면 0m,nN210 \le m^{\prime}, n^{\prime \prime} \le N_{2}-1이고, 0m,nN110 \le m^{\prime \prime}, n^{\prime} \le N_{1}-1이다. (1)(1)의 지수부분을 다음과 같이 나타내자.

ei2πmn/N=ei2π(mN1+m)(nN2+n)/(N1N2)=ei2π(mnN1N2+mnN1+mnN2+mn)/(N1N2)=ei2πmnei2π((mn/N2)+(mn/N1)+(mn/N))=ei2π((mn/N2)+(mn/N1)+(mn/N)) \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)(1)에 대입하면,

a^m=n=0N1ei2πmn/Nan=n=0N21n=0N11ei2π((mn/N2)+(mn/N1)+(mn/N))anN2+n=n=0N21[n=0N11ei2πmn/N1anN2+n]ei2π[(mn/N2)+(mn/N)] \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*}

위 식대로라면 각괄호 안의 [n=0N11]\left[ \sum_{n^{\prime}=0}^{N_{1}-1} \right]를 계산하는데 N1N_{1}번의 연산, 각괄호 밖의 n=0N21\sum_{n^{\prime \prime}=0}^{N_{2}-1}를 계산하는데 N2N_{2}번의 연산이 필요하므로, a^m\hat{a}_{m}을 계산하기 위해서는 총 (N1+N2)(N_{1} + N_{2})번의 연산이 필요함을 알수 있다. a^\hat{\mathbf{a}}를 얻기 위해선 이를 NN번 반복해야하므로 총 N(N1+N2)N(N_{1} + N_{2})의 비용이 들어 N2N^{2}보다 줄어들었음을 확인할 수 있다.

각괄호 안을 잘 보면, N1N_{1}이 다시 합성수일 때 같은 논리를 적용할 수 있음을 알 수 있을 것이다. 따라서 일반적으로 데이터의 길이가 N=N1N2NkN = N_{1} N_{2} \cdots N_{k}와 같은 합성수일 때, 시간복잡도는 다음과 같이 줄어든다.

O(N2)O(N(N1+N2++Nk)) \mathcal{O}(N^{2}) \searrow \mathcal{O}\big( N(N_{1} + N_{2} + \cdots + N_{k}) \big)

이제 NN22의 제곱수 N=2kN = 2^{k}라고 해보자. 그러면 log2N=k\log_{2}N = k이고, N2=2kN^{2} = 2^{k}에서 2k(2k)2^{k}(2k)만큼 줄어드는 것이므로, 시간 복잡도는 다음과 같이 줄어든다.

O(N2)O(2Nlog2N) \mathcal{O}(N^{2}) \searrow \mathcal{O}(2N \log_{2}N)

여담

이는 1965년 Cooley와 Tukey2에 의해 제안되었기 때문에 쿨리-튜키 알고리즘Cooley-Tukey algorithm이라고도 불린다. 하지만 이들이 최초로 발명한 것은 아니다. 가우스도 이와 같은 알고리즘을 연구했지만, 제대로 발표하지 않아 이 사실은 후에 알려지게되었다3.


  1. Gerald B. Folland, Fourier Analysis and Its Applications (1992), p252-253 ↩︎

  2. J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series, Mathematics of Computation 19 (1965), 297-301. ↩︎

  3. 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. ↩︎