고속 푸리에 변환(FFT: Fast Fourier Transform)

[경고] 아래 글을 읽지 않고 "고속 푸리에 변환"을 보면 바보로 느껴질 수 있습니다.

1. 이산 푸리에 변환

[고속 푸리에 변환의 쉬운 이해]

통신 이론과 신호 처리에 필수적인 고속 푸리에 변환(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)의 둘째식과 식 (3b)의 첫째식을 이용해서 $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)

여기서 $k$ = $0,1,\cdots ,M/4-1$이다. 식 (8a)에서 $E_2|E_1(x)$와 $O_2|E_1(x)$는 각각 $E_1(x)$를 구성하는 $x$의 짝수차와 홀수차 항을 모아서 만든 다항식이다. 비슷하게 $E_2|O_1(x)$와 $O_2|O_1(x)$는 각각 $O_1(x)$의 짝수차와 홀수차 항으로 구성한다.

좀더 쉽게 쿨리–투키 알고리즘을 이해하기 위해, 예시적으로 $M$ = $2^1$인 경우부터 출발해 식 (4)를 다시 기술한다.

                  (9)

표본 크기 $M$이 2인 FFT 결과를 가시화하려 나비 모양처럼 생긴 도표를 [그림 1]처럼 도입한다.

[그림 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)로 이름 붙인다. 영 채우기의 조건은 붙여넣은 $f_k$의 크기가 매우 작아서 0으로 근사할 수 있다는 가정이다. 예를 들어, $M_0$개로 영 채우기를 하면 표본수는 $M+M_0$으로 늘고 전체 표본화 시간은 $T$ = $(M+M_0)\mathrm{\Delta }t$로 증가하다. 여기서 $\mathrm{\Delta }t$는 표본화 주기(sampling period)이다. 이로 인해 주파수 영역의 각주파수는 $p_k$ = $k\mathrm{\Delta }p$로 변한다. 여기서 $k$ = $0,1$, $\cdots$, $M+M_0-2,M+M_0-1$; 각주파수 간격은 $\mathrm{\Delta }p$ = $2\pi /[(M+M_0)\mathrm{\Delta }t]$가 되어 이전보다 줄어든다. 그래서 시간 영역의 표본수를 늘리면, 그만큼 주파수 영역의 표본 간격이 축소된다. 하지만 $p_k$의 범위는 $(M+M_0)\mathrm{\Delta }p$ = $2\pi /\mathrm{\Delta }t$이므로, 영 채우기로는 $p_k$의 최대 한계를 늘릴 수 없다. 여기서 $f_s$ = $1/\mathrm{\Delta }t$와 ${\omega }_s$ = $2\pi /\mathrm{\Delta }t$는 각각 표본화 주파수(sampling frequency)와 표본화 각주파수(sampling angular frequency)로 부른다.

주파수 영역에도 영 채우기가 가능하다. 기존 $p_k$의 마지막 값인 $p_{M-1}$이 0에 가깝다고 생각해서, $p_M$ 이상을 0으로 간주하고 원하는 만큼 0을 추가한다. 이 경우 표본화 주기 $\mathrm{\Delta }t$의 고정 유무에 따라 시간 영역에 나타나는 효과가 달라진다. 먼저 $\mathrm{\Delta }t$를 상수라 가정하고 시간 영역의 영 채우기와 동일한 논리를 적용한다. 그러면 영을 채운 회수 $M_0$을 아무리 늘려도 $p_k$의 영역[= ${\omega }_s$ = $2\pi /\mathrm{\Delta }t$]은 증가하지 않고 각주파수 간격 $\mathrm{\Delta }p$[= $2\pi /[(M+M_0)\mathrm{\Delta }t]$]만 좁아진다. 이때 시간 영역에서는 전체 표본화 시간[= $(M+M_0)\mathrm{\Delta }t$]이 주파수 영역에서 영 채우기한 만큼 늘어난다. 반대로 $\mathrm{\Delta }t$가 바뀔 수 있지만 $\mathrm{\Delta }p$는 고정이면, 데칼코마니(décalcomanie)처럼 시간 영역의 영 채우기와 대칭적으로 생각한다. 영 채우기로 인해 표본화 각주파수는 $(M+M_0)\mathrm{\Delta }p$로 늘어나지만 $\mathrm{\Delta }p$ 자체는 고정된다. 이 결과를 시간으로 가져가서 $t_m$ = $m\mathrm{\Delta }t$ = $2\pi m/[(M+M_0)\mathrm{\Delta }p]$로 둔다. 따라서 표본화 주기는 $\mathrm{\Delta }t$ = $2\pi /[(M+M_0)\mathrm{\Delta }p]$로 축소되고, 전체 표본화 시간은 $(M+M_0)\mathrm{\Delta }t$ = $2\pi /\mathrm{\Delta }p$로 변화 없다.

FFT의 기본 알고리즘은 간단해서 우리가 쉽게 프로그램할 수 있다. 그런데 많은 자료를 진짜 빠르게 처리하려면 전문적인 FFT 함수 집합 혹은 라이브러리(library)를 써야 한다. 우리가 쉽게 구할 수 있는 상용 FFT 라이브러리 중에 FFTW(FFT in the West)가 있다[4]. FFTW는 MATLAB에도 쓰일 정도로 잘 검증된 연산을 제공한다. 다만 FFTW는 비상업적인 용도로 쓸 때에만 무료로 이용할 수 있다. NVIDIA GPU를 사용하는 경우는 cuFFT[CUDA(Compute Unified Device Architecture) 기반 FFT]를 사용할 수 있다[5]. 여러 문헌에 나오는 결과를 종합하면, 표본수가 커질수록 cuFFT는 CPU를 쓰는 FFTW보다 훨씬 빠른 성능을 보인다.

FFT는 중앙화 고속 푸리에 변환(centered FFT, CFFT: CDFT의 고속화), 대칭 고속 푸리에 변환(symmetric FFT, SFFT: SDFT의 고속화), 일반화 고속 푸리에 변환(generalized FFT, GFFT: GDFT의 고속화) 등과 같은 변종이 많이 있다. 하지만 CFFT, SFFT를 위한 고속 알고리즘을 다시 제안하지 않고, 이미 구현된 FFT를 바꾸어서 CFFT, SFFT를 구현한다. CFFT는 순환 편이(circular shift) 연산자 $\mathrm{shift}(\cdot )$를 도입해서 FFT로 재구성한다.

                          (12)

여기서 $M_e$ = $\lfloor M/2\rfloor$, $M_o$ = $\lfloor (M+1)/2\rfloor$, $M$ = $M_e+M_o$; $\lfloor x\rfloor$ 혹은 $[x]$는 $x$를 넘지 않는 바닥 함수(floor function), $\mathrm{shift}(\cdot )$는 표본 수열 $f_m$을 $M_e$만큼 편이시킨다. 이 개념을 CDFT(centered DFT) 공식에 넣어서 FFT를 포함하도록 변형한다.

                          (13a)

CFFT의 역변환인 ICFFT(inverse centered FFT)도 IFFT(inverse FFT)와 $\mathrm{shift}(\cdot )$로 계산한다.

                          (13b)

표본 개수 $M$이 홀수인 경우는 SFFT와 CFFT의 결과는 같다. 하지만 $M$이 짝수가 되면 $m$ = $0$을 기준으로 음과 양인 $m$이 대칭적이지 않아서 $+0.5$만큼 $m$을 움직여야 한다. 즉, 홀수인 $M$에서는 $t$ = $m\mathrm{\Delta }t$ ($m$ = $0,\pm 1,\cdots ,\pm M_e$)를 써도 되지만 짝수에서는 $t$ = $(m+0.5)\mathrm{\Delta }t$ ($-M_e\le m\le M_e-1$)로 만든다. 따라서 SFFT는 식 (13a)처럼 FFT와 $\mathrm{shift}(\cdot )$를 활용하지만 FFT 전에 $+0.5$가 만드는 위상 기여를 곱한다.

                          (14a)

여기서 ${\mathrm{\Delta }}_e(M)$은 $M$의 짝수 판별식(discriminant of even number)이다. 역변환인 ISFFT(inverse symmetric FFT)도 FFT로 표현한다.

                          (14b)

이와 같이 FFT는 약방의 감초처럼 우리가 상상하는 다양한 푸리에 변환에 고속화라는 감칠맛을 더한다.  

[참고문헌]

[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)

[5] cuFFT Library User's Guide, cuFFT API Reference, NVIDIA Corp., California, USA, 2022. (방문일 2022-10-04)