[경고] 아래 글을 읽지 않고 "고속 푸리에 변환"을 보면 바보로 느껴질 수 있습니다.
[고속 푸리에 변환의 쉬운 이해]
통신 이론과 신호 처리에 필수적인 고속 푸리에 변환(fast Fourier transform, FFT)은 이산 푸리에 변환(discrete Fourier transform, DFT)의 항들을 우함수(even function)와 기함수(odd function)로 모아서 급수를 빠르게 계산하기 위한 혁신적인 알고리즘이다[1], [2]. FFT는 1805년가우스 28세, 조선 순조 시절 무렵에 가우스Carl Friedrich Gauss(1777–1855)가 먼저 고안해 사용했다는 전설을 가지고 있다.[푸리에 교수가 푸리에 급수(Fourier series)를 공개한 시점이 1807년이다.] 현대적인 FFT는 뼈대 있는 수학자 룽에Carl Runge(1856–1927)의 1903년룽에 47세, 대한제국 시절 기여를 거쳐 미국 수학자 쿨리James Cooley(1926–2016)와 투키John Tukey(1915–2000)가 1965년쿨리 39세, 박정희 정부 시절에 재발견했다[3]. 쿨리와 투키가 FFT를 구현한 방법론은 쿨리–투키 알고리즘(Cooley–Tukey algorithm)으로 부른다. 이 알고리즘의 유도는 기존 DFT에 대한 관찰로부터 시작한다.
(1)여기서 $M$은 표본 개수, $f_m$은 시계열(time sequence), $F_k$는 푸리에 계수(Fourier coefficient)이다. DFT의 모든 항을 항상 이등분할 수 있도록, FFT에서는 보통 $M$을 2의 거듭제곱인 $M$ = $2^N$으로 가정한다. 여기서 $N$은 1이상의 자연수이다. 먼저 식 (1)에서 복소 지수 함수를 $W_M$으로 치환해서 멱급수(power series) 형태로 다시 쓴다.
(2)복소 지수 함수 $W_M$의 거듭제곱은 다음 성질을 가진다.
(3a)
(3b)그 다음에 식 (2)를 이등분해서 DFT 항을 $W_M^k$에 대한 우함수와 기함수인 부분으로 나눈다.
(4)
(5)여기서 $k$ = $0, 1, \cdots, M/2-1$, 원래 DFT 급수를 1번 분해한 $E_1(x)$와 $O_1(x)$는 각각 짝수 및 홀수 밑수를 가진 $f_k$만을 모은 다항식이다. 그러면 식 (3a)의 둘째식을 이용해서 $F_{k+M/2}$를 다시 쓸 수 있다.
(6)여기서 $W_M^{k+M/2}$ = $-W_M^k$이다. 그러면 식 (6)을 다시 계산할 필요없이, 식 (4)에서 구한 $E_1(\cdot)$과 $O_1(\cdot)$의 결과값에서 부호만 ($+$)에서 ($-$)로 바꾸면 된다. 즉, DFT 급수를 우함수와 기함수로 분해함으로써 계산 속도를 2배만큼 증가시킬 수 있다. 이 방식을 염두에 두고 식 (4)를 새롭게 보면, $E_1(x)$와 $O_1(x)$는 항의 개수가 반[= $M/2$]만큼 줄어든 또 다른 DFT 급수이다. 그래서 재귀적인 방법으로 식 (5)를 다시 우함수와 기함수 관련 다항식으로 모아서 계속 연산을 가속화시킨다. 이와 같은 기법이 바로 쿨리–투키 알고리즘의 핵심이. 예를 들어, 재귀 관계에서 $E_1(\cdot)$과 $O_1(\cdot)$은 다음 다항식으로 재분해된다.
(7a)
(7b)
(8a)
(8b)좀더 쉽게 쿨리–투키 알고리즘을 이해하기 위해, 예시적으로 $M$ = $2^1$인 경우부터 출발해 식 (4)를 다시 기술한다.
(9)
[그림 1] 기본 $2 \times 2$ 나비 도표(사진 출처: wikipedia.org)
[그림 2] 변형 $2 \times 2$ 나비 도표
[그림 1]은 나비를 닮아서 나비 도표(butterfly diagram)라고 명한다. 기본 나비 도표의 $f_1$에 $W$를 곱한 변형은 [그림 2]처럼 화살표를 하나 추가해서 그린다. 표본 크기를 더 증가시킨 $M$ = $4$인 경우도 식 (4)를 바탕으로 우함수와 기함수를 분리한다.
(10)
[그림 3] $4 \times 4$ 나비 도표
식 (10)을 나비 도표로 나타낸 결과인 [그림 3]은 약간 복잡하지만 명확한 규칙이 보인다. 시계열 $f_k$를 시작할 때, $k$를 이진수로 나타내어 비트를 역전시킨다. 이 과정은 FFT의 비트 역전(bit reversal)이라 이름 붙인다. 그러면 [그림 3]의 입력은 0, 2, 1, 3의 순서로 들어간다. 그 다음에 [그림 1]처럼 $M$ = $2$에 대한 FFT를 각각 수행하고, 이 결과를 모아서 각 화살표에 대해 $M$ = $4$에 대한 FFT 연산을 다시 적용한다. 식 (4)와 (7)은 비트 역전이 필요한 이유를 잘 보여주고 있다. 우리가 홀짝을 나눈 다음에 표본 개수 $M$과 지수 $k$를 1/2하고, 우함수와 기함수 관련 다항식을 각각 홀짝으로 다시 분류한다. 지수 $k$를 줄이는 단계는 2의 거듭제곱으로 나누면서 계속 반복된다. 최종적으로 항이 2개만 남으면 멈춘 후, 나비 도표인 [그림 1]을 이용해 계산한다. 시계열 $f_k$의 $k$가 줄어서 마지막에 선택된 두 항[그림 1에서 0과 1로 표기]은 $k$를 이진수로 표시할 때 최상위 비트(most significant bit, MSB)가 0이나 1인 경우이다. [그림 1]을 포함한 [그림 3]에서 나비 도표의 둘째 단계는 MSB의 그 다음 비트가 0이나 1을 가진 $k$를 처리한다. 따라서 FFT에서는 [그림 3]처럼 비트 역전을 통해서 입력 시계열 $f_k$를 두 개씩 모아야 한다.
[그림 4] $8 \times 8$ 나비 도표
[그림 4]는 [그림 3]을 2개 나열해서 $8 \times 8$ 나비 도표를 구현한 결과를 보여준다. 시계열 $f_k$의 $k$는 비트 역전을 통해서 결정한다. 입력이 8개라서 우함수와 기함수 부분인 다항식 $E_1(\cdot)$과 $O_1(\cdot)$은 각각 4개의 표본을 가지고 있다. 식 (4)를 사용해서 유도한 $M$ = $8$인 FFT 관계식은 다음과 같다.
(11)여기서 $k$ = $0, 1, 2, 3$, $E_1(\cdot)$과 $O_1(\cdot)$은 [그림 3]으로 계산한다.
표본 개수가 $M$ = $2^N$이 아니더라도 [1]을 참고해서 FFT를 수행할 수 있다. 하지만 $M$이 2의 거듭제곱이 아니면 홀짝 모으기가 번거로워진다. 만약 표본이 작아서 $M$ < $2^N$이면, $M$을 증가시켜 강제로 $M$ = $2^N$을 만들고 모르는 $f_k$는 $0$으로 채워서 FFT를 돌린다. 이런 부가적인 과정은 FFT의 영 채우기(zero padding)로 이름 붙인다.
FFT의 기본 알고리즘은 간단해서 우리가 쉽게 프로그램화할 수 있다. 그런데 많은 자료를 진짜 빠르게 처리하려면 전문적인 FFT 함수 집합 혹은 라이브러리(library)를 써야 한다. 우리가 쉽게 구할 수 있는 상용 FFT 라이브러리 중에 FFTW(FFT in the West)가 있다[4]. FFTW는 MATLAB에도 쓰일 정도로 잘 검증된 연산을 제공한다. 다만 FFTW는 비상업적인 용도로 쓸 때에만 무료로 이용할 수 있다.
[참고문헌]
[1] J. W. Cooley and J. W. Tukey, "An algorithm for the machine calculation of complex Fourier series," Math. Comput., vol. 19, no. 90, pp. 297–301, Apr. 1965.
[2] E. O. Brigham and R. E. Morrow, "The fast Fourier transform," IEEE Spectr., vol. 4, no. 12, pp. 63–70, Dec. 1967.
[3] M. T. Heideman, D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Mag., vol. 1, no. 4, pp. 14–21, Oct. 1984.
[4] M. Frigo and S. G. Johnson, "C Subroutine Library for Computing the DFT," FFTW 3.3.10, 2022. (방문일 2022-10-04)
(1)
(3)
(4b)
(7)
(9)
(10)
(1.1)
(2.1)
(3.1b)
(4.1a)
(4.1b)
