12장. 고속 푸리에 변환 (FFT)

길이 N인 DFT를 정의대로 계산하면 곱셈이 N²번 필요합니다. 고속 푸리에 변환(FFT)은 신호를 짝수·홀수 번째로 나누어 되풀이하는 방법으로 이것을 약 N log N번으로 줄입니다. 20세기의 가장 중요한 알고리즘 가운데 하나입니다.

FFT는 새로운 변환이 아닙니다. DFT(4장 정의 4.3)와 똑같은 값을 훨씬 적은 연산으로 계산하는 알고리즘입니다.

12.1 직접 계산의 비용

$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}$을 정의대로 계산하면 $k$마다 곱셈 $N$번, $k$가 $N$개이므로 $N^2$번입니다. $N = 2^{20} \approx 10^6$(오디오 약 24초)이면 약 $1.1 \times 10^{12}$번입니다. 초당 $10^9$번 곱하는 컴퓨터로 18분쯤 걸립니다. FFT는 같은 일을 약 $2 \times 10^7$번, 0.02초에 합니다.

12.2 분할 정복: 기수 2 FFT

$N$이 짝수라 합시다. $w_N = e^{-j2\pi/N}$이라 쓰면 $X[k] = \sum_nx[n]w_N^{kn}$입니다. 합을 짝수 번째와 홀수 번째 표본으로 나눕니다.

$w_N^2 = e^{-j2\pi/(N/2)} = w_{N/2}$이므로 두 합은 각각 길이 $N/2$인 DFT다.

정리 12.1 (FFT 점화식). $E[k]$, $O[k]$를 짝수 표본열 $x[0], x[2], \ldots$과 홀수 표본열 $x[1], x[3], \ldots$의 길이 $N/2$ DFT라 하면, $k = 0, 1, \ldots, N/2 - 1$에 대해

입니다.

증명. 첫 식은 위에서 보였습니다. 둘째 식은 이렇습니다. $E$, $O$는 주기 $N/2$이므로 $E[k + N/2] = E[k]$, $O[k + N/2] = O[k]$이고, $w_N^{k + N/2} = w_N^ke^{-j\pi} = -w_N^k$입니다. $\blacksquare$

길이 $N$인 DFT 하나가 길이 $N/2$인 DFT 둘과 곱셈 $N/2$번, 덧셈 $N$번으로 계산됩니다. $E[k]$와 $w_N^kO[k]$ 한 쌍에서 출력 두 개가 나오는 이 연산을 버터플라이(butterfly)라 합니다. $N = 2^p$이면 같은 일을 반복해 길이 1까지 내려갑니다. 길이 1인 DFT는 자기 자신입니다.

연산량. 길이 $N$의 비용을 $C(N)$이라 하면 $C(N) = 2C(N/2) + cN$입니다. 한 단계의 작업량이 $cN$이고 단계가 $\log_2N$개이므로

입니다. $N = 1024$이면 $N^2 \approx 10^6$, $N\log_2N \approx 10^4$로 100배입니다. $N = 10^6$이면 5만 배입니다. 이 차이 때문에 디지털 신호처리가 실용 기술이 되었습니다.

예제 12.2. $\mathbf{x} = (1, 2, 3, 4)$, $N = 4$, $w_4 = -j$라 합시다. 짝수 표본 $(1, 3)$의 길이 2 DFT는 $E = (1 + 3, 1 - 3) = (4, -2)$, 홀수 표본 $(2, 4)$는 $O = (6, -2)$입니다.

4장 예제 4.5와 같습니다.

그림 12.1. $N = 4$ FFT의 신호 흐름. 화살표 위의 수는 곱하는 값이고(없으면 1), 한 점으로 들어오는 화살표는 더합니다. 입력의 순서가 $x[0], x[2], x[1], x[3]$으로 바뀌어 있습니다(비트 반전 순서).

12.3 구현

점화식을 그대로 옮긴 재귀 구현입니다. 실제로는 라이브러리를 쓰지만 한 번 손으로 따라가 볼 가치가 있습니다.

import numpy as np


def fft(x):                            # len(x) must be a power of 2
    N = len(x)
    if N == 1:
        return np.array(x, dtype=complex)
    E = fft(x[0::2])                   # DFT of even-indexed samples
    O = fft(x[1::2])                   # DFT of odd-indexed samples
    w = np.exp(-2j*np.pi*np.arange(N//2)/N)
    return np.concatenate([E + w*O, E - w*O])


print(fft([1, 2, 3, 4]))               # [10.+0.j -2.+2.j -2.+0.j -2.-2.j]

• 역변환. 역 DFT는 지수의 부호와 $1/N$만 다릅니다. $\mathbf{x} = \frac{1}{N}\operatorname{FFT}(X^*)^*$이므로 같은 코드로 계산합니다.
• 제자리 계산. 실용 구현은 재귀 대신 입력을 비트 반전 순서(첨자를 이진수로 쓰고 뒤집은 순서)로 먼저 재배열한 뒤, 그림 12.1처럼 단계별로 버터플라이를 적용합니다. 추가 메모리가 필요 없습니다.
• 2의 거듭제곱이 아닌 $N$. 3이나 5로 나누는 혼합 기수 알고리즘, 소수 길이를 합성곱으로 바꾸는 방법이 있어 어떤 $N$이든 $O(N\log N)$에 계산됩니다. 다만 작은 소인수만 가진 $N$이 가장 빠릅니다.
• 실수 입력. 켤레 대칭(4장 정리 4.6) 때문에 절반만 계산하면 됩니다(numpy.fft.rfft).

12.4 FFT로 합성곱 계산하기

합성곱 정리(4장 정리 4.13)를 거꾸로 씁니다. 길이 $L$과 $M$인 두 수열의 합성곱을 직접 계산하면 곱셈이 $LM$번이지만, 다음 순서로 하면 $O(N\log N)$입니다.

1. 두 수열 뒤에 0을 채워 길이 $N \ge L + M - 1$로 만듭니다.
2. 각각 FFT합니다.
3. 성분별로 곱합니다.
4. 역 FFT합니다.

1단계의 영 채우기가 필수입니다. DFT의 합성곱은 순환 합성곱이라, 길이가 모자라면 결과의 뒷부분이 앞으로 감겨 들어와 겹칩니다.

예제 12.3 (다항식의 곱). 다항식의 곱은 계수열의 합성곱입니다. $(1 + 2z + 3z^2)(4 + 5z)$의 계수는 $(1, 2, 3) * (4, 5)$로

즉 $4 + 13z + 22z^2 + 15z^3$입니다. FFT로는 $N = 4$로 채워 $(1, 2, 3, 0)$, $(4, 5, 0, 0)$을 변환하고 곱하고 역변환합니다. 큰 정수의 곱셈(자릿수를 계수로 봅니다), 긴 필터의 적용, 상관 계산이 모두 이 방법으로 빨라집니다.

쓰이는 곳. 오디오·영상 처리, 큰 합성곱의 빠른 계산, 스펙트럴 방법, 레이더와 통신 수신기가 모두 FFT를 씁니다. 논문에서 '$O(N\log N)$'이라는 계산량을 보면 대개 FFT를 쓴 것입니다.

연습문제

4
0
1
-2.414
0
0
1
-0.414

19
13
22

204.8

130