FFT 이해하기
고속 푸리에 변환
문제 풀이 판에서 고속 푸리에 변환 (Fast Fourier Transfrom, FFT) 는 매우 중요하고 한 번은 짚고 넘어가야 하는 주제입니다.
예전부터 이걸 이해를 하고 넘어가자 생각하면서 미루고 미뤄오다 한 번 실제로 코드로 옮기면서, 중요한 개념들을 정리하고자 합니다.
이산 푸리에 변환
FFT는 이름처럼 푸리에 변환을 빠르게 하는 알고리즘입니다. 그렇기에 푸리에 변환이 무엇인지부터 알아야 합니다.
푸리에 변환 (Fourier Transform) 은 어떤 신호를 주파수 성분으로 분해하는 방법입니다. 예를 들어, 음악 신호는 여러 주파수의 소리로 이루어져 있습니다. 푸리에 변환을 통해 이 신호를 각 주파수 성분으로 나눌 수 있습니다.
원래 푸리에 변환은 연속적인 신호에 적용됩니다. 하지만 컴퓨터는 이산적인 데이터를 다루기 때문에, 이산 푸리에 변환 (Discrete Fourier Transform, DFT) 을 사용합니다. DFT는 다음과 같이 정의됩니다:
$$X_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i kn/N}$$
왜 이런 식이 나오는지는 3B1B 유튜브 영상 등 좋은 참고 자료가 많으니 참고하시기 바랍니다.
하지만 이 식을 그대로 계산하면 딱 보기에도 $O(N^2)$ 의 시간 복잡도를 가집니다. 이걸 더 빠르게, 구체적으로는 $O(N \log N)$ 으로 계산하는 방법이 바로 FFT입니다.
FFT의 아이디어
FFT의 핵심 아이디어는 분할 정복 (Divide and Conquer, DnC) 입니다. DFT의 계산을 절반으로 나누어 재귀적으로 계산하는 방식입니다.
구체적으로 이를 실행하는 방법 중 하나는 쿨리-튜키 알고리즘 (Cooley-Tukey algorithm) 입니다. 이 알고리즘은 다음과 같은 단계를 거칩니다:
- 입력 신호를 짝수 인덱스와 홀수 인덱스로 나누기
- 각각의 부분 신호에 대해 DFT를 재귀적으로 계산하기
- 결과를 결합하여 최종 DFT를 계산하기
재귀 구현
이런 재귀적인 아이디어의 과정은 코드로 다음과 같이 표현할 수 있습니다:
Ruby 구현
# Cooley-Tukey FFT algorithm
def recursive_fft(data)
n = data.length
return data if n <= 1
# Divide and conquer
even = recursive_fft(data.each_slice(2).map(&:first))
odd = recursive_fft(data.each_slice(2).map(&:last))
combined = Array.new(n)
(0...n / 2).each do |k|
t = Math::E**(-2.0i * Math::PI * k / n) * odd[k] # Twiddle factor
combined[k] = even[k] + t
combined[k + n / 2] = even[k] - t
end
combined
endC 구현
void fft_recursive(_complex *x, int N) { // _complex is a struct for complex numbers, and c_add, c_sub, c_mul, c_exp are functions for complex arithmetic
if (N <= 1)
return; // Base case
// Divide
_complex even[N / 2];
_complex odd[N / 2];
for (int i = 0; i < N / 2; i++) {
even[i] = x[i * 2];
odd[i] = x[i * 2 + 1];
}
// Conquer
fft_recursive(even, N / 2);
fft_recursive(odd, N / 2);
// Combine
for (int k = 0; k < N / 2; k++) {
_complex t = c_exp(-2 * M_PI * k / N); // e^(-2πik/N)
t = c_mul(t, odd[k]);
x[k] = c_add(even[k], t);
x[k + N / 2] = c_sub(even[k], t);
}
}여기서 보이듯, 각각의 구현은 재귀적으로 입력 데이터를 절반으로 나누는 과정에서 $O(\log N)$ 의 깊이를 가지며, 각 단계에서 $O(N)$ 의 작업을 수행하므로 전체 시간 복잡도는 $O(N \log N)$ 이 됩니다.
비재귀 구현
재귀 구현이 일반적으로 직관적이지만, 많은 문제들의 경우 비재귀 구현을 기준으로 세팅된 제한을 가집니다.
비재귀 구현의 경우에는 비트 반전 순서 (bit-reversed order) 로 데이터를 재배열하는 과정이 필요합니다.
비트 반전 순서는 인덱스의 이진 표현을 뒤집는 것을 의미합니다. 예를 들어, 길이가 8인 배열에서 인덱스 3 (011)은 비트 반전 후 6 (110)이 됩니다.
이 과정을 통해 데이터를 재배열한 후, 반복문을 통해 FFT를 수행할 수 있습니다.
Ruby 구현
# Helper function for bit-reversed order permutation
def bit_reversed_indices(n)
indices = (0...n).to_a
bits = Math.log2(n).to_i
indices.map do |i|
reversed = 0
bits.times do |j|
reversed = (reversed << 1) | ((i >> j) & 1)
end
reversed
end
end
# Iterative FFT algorithm
def iterative_fft(data)
n = data.length
# Bit-reversed order permutation
indices = bit_reversed_indices(n)
permuted = indices.map { |i| data[i] }
data = permuted
m = Math.log2(n).to_i
(1..m).each do |s|
m_val = 1 << s
wm = Math::E**(-2.0i * Math::PI / m_val)
(0...n).step(m_val) do |k|
w = 1
(0...(m_val / 2)).each do |j|
t = w * data[k + j + m_val / 2]
u = data[k + j]
data[k + j] = u + t
data[k + j + m_val / 2] = u - t
w *= wm
end
end
end
data
endC 구현
// Helper function for bit-reversed order permutation
void bit_reversed_permutation(_complex *x, int N) {
int bits = (int)log2(N - 1) + 1; // Number of bits needed, considering the floating-point precision
for (int i = 0; i < N; i++) {
int j = 0;
for (int k = 0; k < bits; k++) {
j = (j << 1) | ((i >> k) & 1);
}
if (j > i) {
_complex temp = x[i];
x[i] = x[j];
x[j] = temp;
}
}
}
// Iterative FFT algorithm
void fft_iterative(_complex *x, int N) {
bit_reversed_permutation(x, N);
int m = (int)log2(N - 1) + 1;
for (int s = 1; s <= m; s++) {
int m_val = 1 << s;
_complex wm = c_exp(-2 * M_PI / m_val);
for (int k = 0; k < N; k += m_val) {
_complex w = 1;
for (int j = 0; j < m_val / 2; j++) {
_complex t = c_mul(w, x[k + j + m_val / 2]);
_complex u = x[k + j];
x[k + j] = c_add(u, t);
x[k + j + m_val / 2] = c_sub(u, t);
w = c_mul(w, wm);
}
}
}
}역변환
FFT로 변환을 빠르게 하는 근본적인 목적은 일반적으로 컨볼루션 (Convolution) 을 빠르게 계산하기 위함입니다.
컨볼루션은 두 신호의 겹치는 부분을 곱하고 더하는 연산으로, 신호 $A$ 와 $B$ 의 컨볼루션 $C$ 는 $C[n] = \sum_{m=0}^{N-1} A[m] B[n-m]$ 로 정의됩니다.
일반적로으로 이는 성분별 곱셈과 덧셈을 통해 $O(N^2)$ 의 시간 복잡도를 가집니다. 행렬을 쓰면 편하게 표현할 수 있지만, 계산량은 동일합니다.
하지만 푸리에 변환의 중요한 성질 중 하나는 컨볼루션 정리 (Convolution Theorem) 입니다.
컨볼루션 정리에 따르면, 두 신호의 컨볼루션은 각각의 푸리에 변환을 취한 후 성분별 곱셈을 하고, 다시 역푸리에 변환을 취하는 것과 같습니다. 즉,
$$\begin{align*} C[n] &= A[n] * B[n] \\ \mathcal{F}\{C[n]\} &= \mathcal{F}\{A[n]\} \cdot \mathcal{F}\{B[n]\} \\ \mathcal{F}^{-1}\{\mathcal{F}\{C[n]\}\} &= \mathcal{F}^{-1}\{\mathcal{F}\{A[n]\} \cdot \mathcal{F}\{B[n]\}\} \\ C[n] &= \mathcal{F}^{-1}\{\mathcal{F}\{A[n]\}\} \cdot \mathcal{F}\{B[n]\} \end{align*}$$
이를 통해 컨볼루션을 $O(N \log N)$ 의 시간 복잡도로 계산할 수 있습니다.
여기서 사용되는 역푸리에 변환 (Inverse Fourier Transform, IFT) 은 다음과 같이 정의됩니다:
$$\mathcal{F}^{-1}{X[k]} = \frac{1}{N} \sum_{n=0}^{N-1} X[k] e^{2\pi i k n / N}$$
이 식은 DFT와 거의 동일하지만, 지수의 부호가 반대이고, 결과에 $1/N$ 을 곱하는 점이 다릅니다.
따라서, 이걸 FFT로 빠르게 계산하는 방법도 DFT와 거의 동일하며, 단지 지수의 부호만 반대로 바꾸고 결과를 $N$ 로 나누면 됩니다.
이를 코드로 표현하면 다음과 같습니다:
Ruby 구현
def iterative_ifft(data)
n = data.length
# Conjugate the complex numbers
conjugated = data.map { |x| x.conj }
# Forward FFT
transformed = iterative_fft(conjugated)
# Conjugate again and scale
transformed.map { |x| x.conj / n }
endC 구현
void ifft_iterative(_complex *x, int N) {
// Conjugate the complex numbers
for (int i = 0; i < N; i++) {
x[i].imag = -x[i].imag;
}
// Forward FFT
fft_iterative(x, N);
// Conjugate again and scale
for (int i = 0; i < N; i++) {
x[i].imag = -x[i].imag;
x[i].real /= N;
x[i].imag /= N;
}
}활용도
실생활의 예로써는, 음악 신호 처리, 이미지 처리, 데이터 압축 등의 분야에서 FFT가 널리 사용됩니다.
여기서 본질적으로 들어가면, 아까 설명한 내용으로부터 컨볼루션을 빠르게 계산하는 것을 통해 여러 문제를 해결할 수 있습니다.
즉, 컨볼루션에 해당하는 다항식의 곱셈 연산, 행렬의 곱셈, 큰 수의 곱셈 등에서 FFT를 활용할 수 있습니다.
대표적으로 큰 수 곱셈 문제는 아예 C/C++ 로 FFT를 구현하라고 대놓고 요구하는 문제입니다.