고속 푸리에 변환

1. 개요

고속 푸리에 변환 Fast Fourier Transform
약칭	FFT
분야	수치해석 × 신호처리 × 계산물리
계산 복잡도	O(N log N)
기수-2 복소 곱셈	(N/2)log2N − N + 1
대표 알고리즘	쿨리-투키(Cooley–Tukey) 분할정복
오차 규모	O(u log N) — 직접 DFT 의 O(uN) 보다 작다
주요 응용	스펙트럴 방법, 에발트 합산, 신호처리

N²을 N log N으로 바꾸는 순간, 세상의 절반이 계산 가능해졌다.

고속 푸리에 변환(Fast Fourier Transform, FFT)은 이산 푸리에 변환(DFT)을 $O(N^2)$이 아닌 $O(N \log N)$에 계산하는 알고리즘의 총칭이다. 변환의 결과값 자체는 DFT와 완전히 동일하다 — FFT는 새로운 변환이 아니라, DFT를 빠르게 계산하는 방법일 뿐이다. 이름에 “고속”이 붙은 이유가 그것.

DFT는 시간(또는 공간) 영역의 신호를 주파수 영역으로 갈아 넣는 도구로, 연속 신호를 다루는 푸리에 변환의 유한·이산판이다. 그런데 정의대로 곧이곧대로 계산하면 결과값 하나마다 $N$번의 곱셈이 필요해서 총 $N^2$번. 데이터가 백만 개면 $10^{12}$번이라 슈퍼컴퓨터도 땀을 뻘뻘 흘린다. FFT는 이걸 약 $2 \times 10^7$번으로 줄인다. 대략 5만 배 빨라지는 셈이니, “일단 FFT부터 돌려”가 국룰이 된 데는 이유가 있다.¹

2. 이산 푸리에 변환의 정의

길이 $N$의 복소수 수열 $x_0, x_1, \dots, x_{N-1}$에 대한 DFT는 다음과 같이 정의된다.

$$X_k = \sum_{n=0}^{N-1} x_n \, e^{-2\pi i \, kn / N}, \qquad k = 0, 1, \dots, N-1$$

여기서 $W_N = e^{-2\pi i / N}$을 회전인자(twiddle factor)라 부른다. 이 식을 행렬-벡터 곱으로 보면 $N \times N$ 행렬을 벡터에 곱하는 꼴이라 곧이곧대로는 $O(N^2)$이다. 역변환은 지수의 부호를 뒤집고 $1/N$을 곱하면 된다.

$$x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k \, e^{+2\pi i \, kn / N}$$

정규화 상수를 앞에 어떻게 나눠 붙이느냐는 관례가 라이브러리마다 달라서, FFTW는 아예 아무 정규화도 하지 않는다(정변환·역변환을 연달아 걸면 $N$배가 된다). 라이브러리를 바꿨는데 진폭이 $N$배 틀어졌다면 대개 이 관례 차이다.

3. 쿨리-투키 분할정복

FFT의 핵심 아이디어는 회전인자의 대칭성과 주기성을 이용한 분할정복이다. 가장 널리 쓰이는 기수-2(radix-2) 쿨리-투키 알고리즘은 $N$이 2의 거듭제곱일 때, 수열을 짝수 인덱스와 홀수 인덱스로 반으로 쪼갠다.

$$X_k = \underbrace{\sum_{m=0}^{N/2-1} x_{2m} W_{N/2}^{km}}_{E_k} + W_N^k \underbrace{\sum_{m=0}^{N/2-1} x_{2m+1} W_{N/2}^{km}}_{O_k}$$

즉 크기 $N$짜리 DFT 하나가 크기 $N/2$짜리 DFT 두 개($E_k$, $O_k$)로 쪼개진다. 게다가 $W_N^{k+N/2} = -W_N^k$라는 대칭성 덕분에 $X_k$와 $X_{k+N/2}$를 한 번에 얻는다.

$$X_k = E_k + W_N^k O_k, \qquad X_{k+N/2} = E_k - W_N^k O_k$$

이 쪼개기를 재귀적으로 $\log_2 N$번 반복하면, 각 단계마다 $O(N)$의 연산이 들어 전체 $O(N \log N)$이 된다. 분할정복의 교과서적 사례로, 컴퓨터과학 시간에도 단골로 등장한다.

4. 버터플라이 연산

위의 $X_k = E_k + W_N^k O_k$, $X_{k+N/2} = E_k - W_N^k O_k$ 한 쌍을 데이터 흐름도로 그리면 나비 날개처럼 생겨서 버터플라이(butterfly)라 부른다. FFT의 실제 구현은 이 버터플라이를 격자처럼 층층이 쌓아 올린 것에 불과하다. 입력을 비트 반전(bit-reversal) 순서로 재배치하면 추가 메모리 없이 제자리(in-place)에서 계산할 수 있어, 하드웨어 구현에도 유리하다.²

버터플라이 하나는 복소수 곱셈 한 번과 덧셈·뺄셈 두 번으로 끝난다. 단계마다 $N/2$개의 버터플라이가 있고 단계가 $\log_2 N$개이니, 곱셈 횟수는 대략 $\frac{N}{2}\log_2 N$이다.

4.1. 검산: 곱셈 횟수의 정확한 값

“대략 $\frac{N}{2}\log_2 N$“은 상한이고, 각 단계의 $k = 0$ 버터플라이는 $W_N^0 = 1$이라 곱셈이 아예 필요 없다. 이 자명한 경우만 빼고 실제로 세면 딱 떨어지는 닫힌 형태가 나온다.

$$M(N) = \frac{N}{2}\log_2 N - N + 1$$

기수-2 재귀를 계수기로 감싸 $N = 2 \dots 16384$($p = 1 \dots 14$) 전부 확인하면 14개 값이 모두 정확히 일치한다($N=1024$에서 4097, $N=16384$에서 98305). $W^{N/4} = -i$처럼 곱셈 없이 처리되는 회전인자를 더 걷어내면 실제 고성능 구현의 횟수는 이보다 더 줄어든다 — $N^2$과의 격차는 여기서 벌어진다.

4.2. 검산: 결과가 정말 DFT와 같은가

정의식 그대로의 $O(N^2)$ DFT와 기수-2 FFT를 같은 무작위 복소 입력에 돌려 비교하면 최대 상대차가 $N=16$에서 $3.3\times10^{-15}$, $N=1024$에서 $2.8\times10^{-13}$이다. 둘은 같은 변환이고, 차이는 반올림뿐이다. 실수 입력에 대한 켤레 대칭성 $X_{N-k} = \overline{X_k}$도 $N=256$에서 $5.4\times10^{-14}$ 이내로 성립하며, 순환 합성곱을 직접 계산한 것과 역FFT(FFT(a)·FFT(b))로 계산한 것의 최대 상대차는 $N=64$에서 $3.9\times10^{-16}$이다(합성곱 정리가 수치적으로도 성립한다).

5. 반올림 오차: 빠른 쪽이 정확하기까지 하다

FFT가 DFT보다 연산을 덜 하므로 오차도 덜 쌓인다는 건 직관에 반하지만 사실이다. 정규화 오차 상한이 직접 합산은 $O(uN)$, 기수-2 FFT는 $O(u\log N)$ 규모다. 회전인자를 양쪽 모두 $kn \bmod N$으로 축약해 합산 오차만 분리하고, 카한 합산으로 만든 준기준해에 대한 상대 RMS 오차를 재면 이렇게 나온다.

$N$	직접 DFT	FFT	비
32	1.94e-16	2.37e-16	0.82
128	3.93e-16	3.18e-16	1.24
512	7.93e-16	3.63e-16	2.18
1024	1.07e-15	3.98e-16	2.69

직접 DFT의 오차는 $N$이 32배 커지는 동안 5.5배 자랐다 — $\sqrt{32} = 5.66$이니 $\sqrt{N}$ 성장과 맞는다. FFT 쪽은 $2.4\times10^{-16}$에서 $4.0\times10^{-16}$으로 거의 평평하다. 다만 격차 자체는 2~3배 수준으로 온건하다는 점도 같이 기억해야 한다. FFT를 쓰는 이유는 여전히 속도이고, 정확도는 보너스다.

정작 자릿수를 날려 먹는 건 알고리즘이 아니라 회전인자를 만드는 방식이다. $e^{-2\pi i kn/N}$을 $kn$을 줄이지 않고 그대로 계산하면 각도의 절대오차가 $kn \sim N^2$에 비례해 커진다. 같은 직접 DFT의 상대오차가 $N=1024$에서 $1.7\times10^{-13}$까지 뛰어 축약판보다 440배 나빠졌다. 회전인자는 $kn \bmod N$으로 줄이거나 표로 미리 만들어 두는 것이 정석이다(부동소수점 연산).³

6. 실수 신호와 변형들

현실의 신호는 대개 실수다. 실수 입력은 위에서 확인한 켤레 대칭성을 가져 계산의 절반이 낭비되는데, 이를 활용한 실수 FFT(rFFT)는 복소 FFT보다 약 두 배 빠르고 메모리도 절반이다. $N$이 2의 거듭제곱이 아니어도 겁먹을 필요는 없다 — 혼합 기수(mixed-radix)나 소인수 분해(prime-factor) 알고리즘, 그리고 임의의 $N$(특히 큰 소수)을 처리하는 블루스타인 알고리즘이 있어서, FFTW 같은 라이브러리는 아무 길이나 던져도 알아서 최적 경로를 잡아준다.⁴ 그래도 $N$을 2·3·5의 곱으로 맞춰 주면 눈에 띄게 빠르다.

관련 변환으로는 이산 코사인 변환(DCT, JPEG의 심장), 이산 사인 변환(DST) 등이 있으며, 모두 FFT로 빠르게 계산된다. 유한 구간 신호를 그냥 변환하면 양 끝의 불연속이 스펙트럼 전체로 새는 누설(leakage)이 생기므로, 실측 스펙트럼 분석에서는 윈도 함수를 씌우는 것이 기본이다.

7. 응용

FFT가 없었다면 현대 계산과학의 상당 부분이 멈춘다.

• 스펙트럴 방법 — 편미분방정식을 푸리에 기저로 전개해 미분을 곱셈으로 바꾼다. 직접수치모사의 주력 도구.
• 푸아송 솔버 — 주기 경계의 포아송 방정식은 푸리에 공간에서 대각화되어, FFT 두 번이면 즉시 풀린다. 주기경계조건 아래 유동 해석의 단골.
• 에발트 합산 — 분자동역학의 장거리 쿨롱 상호작용을 PME(Particle-Mesh Ewald)로 계산할 때 격자 위에서 FFT를 돌린다.
• 신호처리 — 스펙트럼 분석, 필터링, 그리고 두 신호의 합성곱을 곱셈으로 바꾸는 합성곱 정리의 핵심.
• 큰 수 곱셈 — 자릿수를 계수로 보면 곱셈이 합성곱이 되어, 초거대 정수 곱셈도 FFT로 처리한다.

푸리에 급수·테일러 급수처럼 함수를 기저로 전개하는 발상이지만, FFT는 그 전개를 실제로 계산 가능하게 만든 실용의 승리다.

8. 관련 문서

• 푸리에 변환 · 이산 코사인 변환 · 윈도 함수
• 합성곱 정리 · 블루스타인 알고리즘
• 스펙트럴 방법 · 직접수치모사
• 에발트 합산 · 분자동역학
• 포아송 방정식 · 주기경계조건
• 테일러 급수 · 직교다항식
• 부동소수점 연산

9. Footnotes

1. 동가라와 설리번이 Computing in Science & Engineering 2000년 1·2월호에 꼽은 “20세기의 알고리즘 톱 10”에 FFT가 당당히 한 자리를 차지했다. 나머지 아홉 개를 다 합쳐도 인용 횟수로는 FFT를 못 이긴다는 우스갯소리가 있다. ↩

2. 쿨리(James Cooley)와 투키(John Tukey)가 1965년 Mathematics of Computation 에 논문을 낸 게 공식적인 시작이지만, 사실 그 뿌리는 1805년 가우스까지 거슬러 올라간다. 가우스는 소행성 궤도의 삼각 보간을 하다 이 알고리즘을 먼저 적어 두었는데, 라틴어 노트에 묻어 두고 출판을 안 해서 150년간 잊혀 있었다(하이데만·존슨·버러스가 1984년에 이 계보를 정리했다). 논문을 안 내면 이렇게 된다. ↩

3. 그래서 “FFT가 DFT보다 정확하다”를 자랑할 때 실제로 자랑하는 것의 절반은 회전인자를 재귀 구조상 작은 각도로만 쓰게 된다는 점이다. 알고리즘을 바꾸는 것보다 각도를 축약하는 게 먼저다. ↩

4. FFTW의 이름은 “Fastest Fourier Transform in the West”의 약자다. 개발자들의 겸손함이 돋보인다. 참고로 정말 빠르긴 하다 — 실행 전에 여러 분해 경로를 실측해 보고 계획(plan)을 고르는 방식으로, 기계마다 다른 캐시 구조까지 알아서 흡수한다. ↩