FFT 이해하기

#blog#computer-science#math#c#ruby

고속 푸리에 변환

문제 풀이 판에서 고속 푸리에 변환 (Fast Fourier Transfrom, FFT) 는 매우 중요하고 한 번은 짚고 넘어가야 하는 주제입니다.

예전부터 이걸 이해를 하고 넘어가자 생각하면서 미루고 미뤄오다 한 번 실제로 코드로 옮기면서, 중요한 개념들을 정리하고자 합니다.

이산 푸리에 변환

FFT는 이름처럼 푸리에 변환을 빠르게 하는 알고리즘입니다. 그렇기에 푸리에 변환이 무엇인지부터 알아야 합니다.

푸리에 변환 (Fourier Transform) 은 어떤 신호를 주파수 성분으로 분해하는 방법입니다. 예를 들어, 음악 신호는 여러 주파수의 소리로 이루어져 있습니다. 푸리에 변환을 통해 이 신호를 각 주파수 성분으로 나눌 수 있습니다.

원래 푸리에 변환은 연속적인 신호에 적용됩니다. 하지만 컴퓨터는 이산적인 데이터를 다루기 때문에, 이산 푸리에 변환 (Discrete Fourier Transform, DFT) 을 사용합니다. DFT는 다음과 같이 정의됩니다:

Xk=∑n=0N−1xne−2πikn/NX_k = \sum_{n=0}^{N-1} x_n e^{-2\pi i kn/N}

왜 이런 식이 나오는지는 3B1B 유튜브 영상 등 좋은 참고 자료가 많으니 참고하시기 바랍니다.

하지만 이 식을 그대로 계산하면 딱 보기에도 O(N2)O(N^2) 의 시간 복잡도를 가집니다. 이걸 더 빠르게, 구체적으로는 O(Nlog⁡N)O(N \log N) 으로 계산하는 방법이 바로 FFT입니다.

FFT의 아이디어

FFT의 핵심 아이디어는 분할 정복 (Divide and Conquer, DnC) 입니다. DFT의 계산을 절반으로 나누어 재귀적으로 계산하는 방식입니다.

구체적으로 이를 실행하는 방법 중 하나는 쿨리-튜키 알고리즘 (Cooley-Tukey algorithm) 입니다. 이 알고리즘은 다음과 같은 단계를 거칩니다:

  1. 입력 신호를 짝수 인덱스와 홀수 인덱스로 나누기
  2. 각각의 부분 신호에 대해 DFT를 재귀적으로 계산하기
  3. 결과를 결합하여 최종 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
  end

C 구현

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(\log N) 의 깊이를 가지며, 각 단계에서 O(N)O(N) 의 작업을 수행하므로 전체 시간 복잡도는 O(Nlog⁡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
end

C 구현

// 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) 을 빠르게 계산하기 위함입니다.
컨볼루션은 두 신호의 겹치는 부분을 곱하고 더하는 연산으로, 신호 AA 와 BB 의 컨볼루션 CC 는 C[n]=∑m=0N−1A[m]B[n−m]C[n] = \sum_{m=0}^{N-1} A[m] B[n-m] 로 정의됩니다.

일반적로으로 이는 성분별 곱셈과 덧셈을 통해 O(N2)O(N^2) 의 시간 복잡도를 가집니다. 행렬을 쓰면 편하게 표현할 수 있지만, 계산량은 동일합니다.
하지만 푸리에 변환의 중요한 성질 중 하나는 컨볼루션 정리 (Convolution Theorem) 입니다.

컨볼루션 정리에 따르면, 두 신호의 컨볼루션은 각각의 푸리에 변환을 취한 후 성분별 곱셈을 하고, 다시 역푸리에 변환을 취하는 것과 같습니다. 즉,

C[n]=A[n]∗B[n]F{C[n]}=F{A[n]}⋅F{B[n]}F−1{F{C[n]}}=F−1{F{A[n]}⋅F{B[n]}}C[n]=F−1{F{A[n]}}⋅F{B[n]}\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(Nlog⁡N)O(N \log N) 의 시간 복잡도로 계산할 수 있습니다.

여기서 사용되는 역푸리에 변환 (Inverse Fourier Transform, IFT) 은 다음과 같이 정의됩니다:

F−1{X[k]}=1N∑n=0N−1X[k]e2πikn/N\mathcal{F}^{-1}\{X[k]\} = \frac{1}{N} \sum_{n=0}^{N-1} X[k] e^{2\pi i k n / N}

이 식은 DFT와 거의 동일하지만, 지수의 부호가 반대이고, 결과에 1/N1/N 을 곱하는 점이 다릅니다.
따라서, 이걸 FFT로 빠르게 계산하는 방법도 DFT와 거의 동일하며, 단지 지수의 부호만 반대로 바꾸고 결과를 NN 로 나누면 됩니다.

이를 코드로 표현하면 다음과 같습니다:

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 }
end

C 구현

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를 구현하라고 대놓고 요구하는 문제입니다.

참고 자료