이산 푸리에 변환과 FFT

주기 데이터의 삼각 보간

주기를 $2\pi$로 맞추고 한 주기를 $N$등분한 점 $x_j = 2\pi j/N$($j = 0, \dots, N - 1$)에서 값 $f_j$를 안다고 합시다. 복소 지수 $e^{ikx} = \cos kx + i\sin kx$로 쓰면 계산이 훨씬 간단해집니다. 핵심은 등간격 점에서 지수 함수들이 이산적으로 직교한다는 사실입니다.

7장의 언어로 말하면 등간격 점 위의 내적 $\langle u, v\rangle = \sum_j u(x_j)\overline{v(x_j)}$에서 $e^{ikx}$들이 직교 기저입니다. 그래서 보간 계수가 정규 방정식을 풀 필요 없이 내적 하나로 나옵니다. 단, $k$와 $k + N$은 표본점에서 구별되지 않습니다($e^{i(k+N)x_j} = e^{ikx_j}$). 이것이 3절의 앨리어싱입니다.

DFT의 정의, 행렬, 단위근

표본의 위치 $x_j$를 잊고 숫자의 배열만 남기면 정의가 깔끔해집니다.

$\omega_N$은 1의 원시 $N$제곱근, 곧 복소 평면의 단위원을 $N$등분한 점입니다(MathBook: 복소수는 회전이다). $F_N$의 $k$번째 행은 단위원을 시계 방향으로 한 걸음에 $k/N$바퀴씩 도는 수열이고, $X_k$는 신호를 그 회전과 짝지어 더한 값입니다. 신호에 그 빠르기로 도는 성분이 있으면 화살표들이 한 방향으로 모여 $|X_k|$가 커지고, 없으면 서로 상쇄됩니다. 역변환이 켤레를 쓰는 이유는 정리 8.1 그대로입니다. $\overline{F_N}F_N = NI$이므로 $F_N/\sqrt N$은 유니터리 행렬입니다. 여기서 두 가지 성질이 바로 나옵니다.

표본 간격이 $\Delta t$초이면 $X_k$는 주파수 $k/(N\Delta t)$ Hz 성분입니다. 곧 주파수 눈금의 간격(주파수 분해능)은 $\Delta f = 1/(N\Delta t) = 1/T$, 관측 시간의 역수입니다. 실수 신호에서는 $k$와 $N - k$가 켤레라서 $k = 0, \dots, N/2$만 보면 충분하고, 진폭 $A$인 사인파 $A\cos(2\pi k_0j/N)$는 $|X_{k_0}| = |X_{N-k_0}| = AN/2$로 나타납니다. 그래서 진폭을 읽을 때는 $2|X_k|/N$을 그립니다.

표본화와 앨리어싱, 나이퀴스트

표본화 주파수 $f_s = 1/\Delta t$로 $\cos(2\pi f t)$를 뜨면 표본은 $\cos(2\pi fj/f_s)$입니다. $f$에 $f_s$의 정수배를 더하거나 부호를 바꿔도 표본은 같습니다($\cos(2\pi(f + mf_s)j/f_s) = \cos(2\pi fj/f_s + 2\pi mj)$). 그래서 표본이 보여 줄 수 있는 주파수는 $[0, f_s/2]$뿐이고, 그 밖의 주파수는 이 구간 안으로 접혀 들어옵니다. 이것이 앨리어싱(aliasing)입니다. 7 Hz를 $f_s = 10$으로 뜨면 $|7 - 10| = 3$ Hz로 보입니다.

겉보기 주파수는 $f_{\text{alias}} = |f - f_s\cdot\mathrm{round}(f/f_s)|$입니다. 실무의 처방은 두 가지입니다. 표본화 전에 아날로그 저역 통과 필터(앤티앨리어싱 필터)로 $f_s/2$ 위의 성분을 지우거나, 관심 대역보다 충분히 빠르게 뜹니다. CD 오디오가 사람의 가청 한계 약 20 kHz의 두 배를 조금 넘는 44.1 kHz를 쓰는 이유입니다. 일단 표본에 섞여 들어온 앨리어스는 어떤 계산으로도 되돌릴 수 없습니다(MathBook: 연속을 격자로 자르기).

쿨리–튜키 FFT: 짝·홀 분할과 나비 연산

$N$이 짝수일 때 합을 짝수 번째와 홀수 번째로 나눕니다. $\omega_N^2 = \omega_{N/2}$이므로

이를 $N = 2^p$까지 되풀이하는 것이 기수 2 쿨리–튜키 FFT(Cooley–Tukey FFT)입니다. 비용을 $T(N)$이라 하면 $T(N) = 2T(N/2) + cN$, $T(1) = 0$이므로 $T(N) = cN\log_2N$입니다. 정확히는 $\log_2N$단계마다 나비가 $N/2$개, 복소 곱셈 $\frac N2\log_2N$번입니다($\omega^0 = 1$ 같은 자명한 곱셈을 빼면 더 적다).

재귀를 펼쳐 놓으면 입력이 $x_0, x_4, x_2, x_6, x_1, x_5, x_3, x_7$($N = 8$) 순서로 맨 아래에 놓입니다. 짝·홀 분할을 세 번 했으니 인덱스를 2진수로 쓰고 비트를 뒤집은 순서입니다($1 = 001_2 \to 100_2 = 4$). 그래서 실제 구현은 먼저 배열을 비트 반전 순서로 바꾼 뒤, 길이 2, 4, 8, …의 나비를 제자리에서 계산합니다. 추가 메모리가 필요 없습니다.

함수 FFT(x)                             // 제자리, N = 2^p
  for i = 0 … N−1                        // 1단계: 비트 반전 순서로 재배열
    j ← bitreverse(i, p)
    if i < j then swap(x[i], x[j])
  len ← 2
  while len ≤ N                         // log₂N 단계
    w ← exp(−2πi / len)                   // 이 단계의 회전 인자 밑
    for s = 0, len, 2·len, … < N          // 블록마다
      t ← 1
      for k = 0 … len/2 − 1              // 나비 len/2개
        a ← x[s+k];  b ← t · x[s+k+len/2]
        x[s+k] ← a + b;  x[s+k+len/2] ← a − b
        t ← t · w
    len ← 2·len
  반환 x                                // 역변환: w의 부호를 바꾸고 N으로 나눈다

회전 인자 $t \leftarrow t\cdot w$를 반복 곱으로 만들면 반올림 오차가 쌓이므로, 정밀한 구현은 미리 $\cos$, $\sin$ 표를 만들어 둡니다. 그렇게 하면 FFT의 상대 오차는 $O(\varepsilon_{\text{mach}}\log N)$으로, 정의대로 계산한 DFT보다 오히려 정확합니다. $N$이 2의 거듭제곱이 아니어도 $N = N_1N_2$로 인수분해되면 같은 분할이 되고(혼합 기수), 소수 $N$도 블루스타인 방법 등으로 $O(N\log N)$에 계산됩니다. FFTW 같은 라이브러리는 기계에 맞는 분할 순서를 실행 중에 골라 씁니다.

빠른 합성곱과 다항식·큰 수 곱셈

다항식 $a(x) = \sum a_ix^i$와 $b(x) = \sum b_jx^j$의 곱의 계수는 $c_k = \sum_{i + j = k}a_ib_j$, 곧 계수 배열의 (선형) 합성곱(convolution) $\mathbf c = \mathbf a * \mathbf b$입니다. 정수 $123 \times 45$는 $x = 10$에서의 다항식 곱 $(3 + 2x + x^2)(5 + 4x)$입니다. 직접 계산은 $O(nm)$인데, DFT는 합성곱을 곱셈으로 바꿉니다.

원형 합성곱은 끝을 넘어간 항이 앞으로 감겨 들어옵니다. 다항식 곱셈에 필요한 것은 선형 합성곱이므로, 두 배열 뒤에 0을 붙여 길이를 $N \ge n + m - 1$로 늘리면 감겨 들 항이 모두 0이 되어 두 합성곱이 같아집니다.

함수 FFT_곱셈(a, b)                       // 길이 n, m의 계수 배열
  N ← n + m − 1 이상인 가장 작은 2의 거듭제곱
  A ← FFT(a 뒤에 0을 붙여 길이 N)
  B ← FFT(b 뒤에 0을 붙여 길이 N)
  C ← A ⊙ B                                 // N번의 곱셈
  c ← IFFT(C)의 앞 n + m − 1개
  반환 c                                    // 정수 계수면 반올림. 큰 수면 받아올림 처리

비용은 FFT 세 번과 곱셈 $N$번, 합쳐서 $O(N\log N)$입니다. 큰 수 곱셈에서는 합성곱 값이 $n\cdot 9^2$까지 커지므로 부동소수점 반올림이 정수 반올림을 망치지 않도록 한 칸에 담는 자릿수를 조절하거나, 반올림 오차가 없는 수론적 변환(유한체에서의 DFT)을 씁니다. 같은 원리로 신호의 FIR 필터, 이미지 흐림(MathBook: 모든 것은 주파수의 합), 확률 분포의 합의 분포를 FFT로 계산합니다. 필터 길이가 수십 이하이면 직접 합성곱이 더 빠르고, 수백 이상이면 FFT가 이깁니다.

스펙트럼 누설, 창 함수, 제로 패딩

DFT는 관측한 $N$개 표본이 한 주기인 주기 신호라고 가정합니다. 주기가 정수 개 들어맞지 않으면 끝과 처음이 이어지는 곳에 불연속이 생기고, 불연속의 계수는 $1/k$로 느리게 줄었습니다(7장 정리 7.10). 수학적으로는 무한한 사인파에 길이 $N$의 사각 창을 곱한 것이고, 곱셈은 주파수 영역에서 창의 스펙트럼 $|\sin(\pi k)/\sin(\pi k/N)|$과의 합성곱입니다. 사인파 하나의 에너지가 모든 칸으로 새어 나가는 이 현상을 스펙트럼 누설(spectral leakage)이라 합니다.

처방은 신호의 양 끝을 부드럽게 0으로 줄이는 창 함수(window function) $w_j$를 곱한 뒤 DFT하는 것입니다. 끝의 불연속이 사라져 누설이 빠르게 줄어드는 대신, 봉우리의 폭(주엽)이 넓어집니다. 주파수 분해능과 누설 억제를 맞바꾸는 셈입니다.

표의 값은 창의 스펙트럼을 아주 촘촘히 계산해 잰 것입니다. "칸 사이 손실"은 주파수가 두 칸의 정중앙(예: 10.5)에 있을 때 봉우리가 얼마나 낮아지는지입니다. 사각 창으로 $N = 64$, 10.5 Hz를 재면 진폭 1이 0.646(−3.79 dB)으로 읽히고 약 20칸 떨어진 곳에서도 −26 dB가 남지만, 한 창은 0.849(−1.42 dB)로 읽히고 같은 곳에서 −93 dB까지 떨어집니다.

제로 패딩(zero padding)은 표본 뒤에 0을 붙여 더 긴 DFT를 하는 것입니다. 같은 연속 스펙트럼(창을 곱한 신호의 이산 시간 푸리에 변환)을 더 촘촘한 간격으로 표본하는 효과라, 봉우리의 위치와 높이를 더 정확히 읽을 수 있습니다. 하지만 분해능, 곧 가까운 두 주파수를 구별하는 능력은 관측 시간 $T$로 정해진 $1/T$에서 나아지지 않습니다. 정보는 0을 붙인다고 늘지 않습니다. 두 성분이 1칸보다 가까우면 더 오래 측정하는 수밖에 없습니다.

응용: 오디오, 이미지 압축, 스펙트럴 미분

오디오 스펙트럼. 음악 분석기와 음성 인식은 신호를 짧은 구간(예: 44.1 kHz에서 2048표본, 약 46 ms)으로 잘라 한 창을 곱하고 FFT를 반복합니다. 이것이 단시간 푸리에 변환이고, 그 크기를 시간–주파수 평면에 그린 것이 스펙트로그램입니다. 구간 길이가 주파수 분해능($44100/2048 \approx 21.5$ Hz)과 시간 분해능을 맞바꿉니다.

이미지 압축 맛보기. JPEG은 영상을 $8\times8$ 블록으로 나눠 각 블록에 2차원 이산 코사인 변환(DCT)을 합니다. DCT는 신호를 좌우 대칭으로 늘여 붙인 뒤의 DFT라, 블록 경계의 불연속이 없어 누설이 적고 에너지가 낮은 주파수 몇 개에 모입니다. 높은 주파수 계수를 거칠게 양자화하거나 버려서 용량을 줄입니다(MathBook: 모든 것은 주파수의 합). 7장의 체비쇼프 계수가 바로 DCT로 계산되었다는 것도 같은 이야기입니다.

스펙트럴 미분. 주기 함수를 미분하면 $e^{ikx}$의 계수가 $ik$배가 됩니다. 그래서 FFT → $ik$ 곱하기 → 역 FFT로 표본점에서의 도함수를 얻습니다(짝수 $N$에서 나이퀴스트 계수는 0으로 둔다). 매끄러운 주기 함수에서는 7장의 해석적 함수처럼 계수가 기하급수로 줄기 때문에 오차가 $N$에 대해 지수적으로 줄어듭니다. 9장의 차분 공식과 비교해 봅시다.

그 밖에. 7장에서 미룬 등간격 주기 데이터의 삼각 최소제곱은 DFT 계수를 $|k| \le n$까지만 남기는 것입니다(이산 직교성 덕분에 정규 방정식이 대각). 2차원 FFT는 푸아송 방정식을 주기 경계에서 $O(N^2\log N)$에 풀고(16장), MRI는 측정한 k-공간 데이터를 2차원 역 FFT로 영상으로 바꿉니다.

이 도구가 쓰이는 곳