FFT (Fast Fourier Transform)

개요

FFT는 이산적인 영역에서 하는 푸리에 변환인 DFT(Discrete Fourier Transform)을 분할정복 기법을 이용하여 $\Theta (n\,log\, n)$ 시간 내에 수행하는 알고리즘 입니다. 이를 보통 PS에서는 두 다항식의 곱셈을 $\Theta (n\,log\, n)$ 에 처리하기 위해 사용합니다. 이 글을 이해하기 위한 사전지식으로는 기초적인 프로그래밍 지식과 복소수 지수에 대한 이해가 필요합니다.

 

DFT (Discrete Fourier Transform)

우선 DFT의 정의부터 살펴봅시다. 수열 $a_{0}, \,a_{1}, \,a_{2} , \,a_{3},\,...\, a_{n-1}$ 에 대하여,

$\omega_{N} = e^{-2\pi i/N} $,  $\hat{a}_{n} = \sum_{k=0}^{N-1}a_{k}w_{N}^{kn}$ 로 정의되는데,

사실 이 정의가 뜻하는 바는, 같은 수열에 대해 다항식, $P(x) = a_{0}+a_{1}x+a_{2}x^2+a_{3}x^3...\, +a_{n-1}x^{n-1}$에 대하여,

$ \hat{a}_{n} = P(w_{N}^{n}) $입니다.

 

그럼 이걸 가지고 뭘 할 수 있을까요?

우선 두 다항식의 곱셈에 대한 문제를 생각해봅시다.

다항식 $A(x)$와 $B(x)$의 계수들에 대해, $\mathbf{DFT}(A),\,\mathbf{DFT}(B)$을 취한 후, 각 결과값들을 곱해주면, $\mathbf{DFT}(AB)$의 결과값인 두 다항식을 곱한 후 DFT를 한 결과와 같음이 자명합니다.

왜냐하면 $A(x), \,B(x)$에 같은 값을 대입하고, 그 값을 곱해주면, $A(x)B(x)$에 그 값을 대입한 결과와 같은 값이 나오기 때문입니다.

그렇다면 만약 역변환인 IDFT(Inverse Discrete Fourier Transform)을 취할 수 있으면(IDFT에 대한 설명은 나중에 나옵니다),

두 다항식을 곱한 결과값을 구할 수 있게 됩니다. 

즉, $\mathbf{IDFT(}\mathbf{DFT}(A)\cdot\mathbf{DFT}(B))$ 의 결과값은 두 다항식을 곱한 다항식의 계수가 됩니다.

하지만 문제가 있는데, DFT를 정의대로 계산하면, $\Theta (n^2)$ 시간이 걸리게 됩니다.

여기서 우리는 DFT를 $\Theta (n\,log\, n)$에 계산하는 FFT를 사용하게 됩니다.

 

FFT(Fast Fourier Transform)

DFT를 빨리 계산하기 위해서 우리는 우함수와 기함수의 성질을 이용합니다.

짝수 차수만 있는 우함수 $E(x)$와 홀수 차수만 있는 기함수 $O(x)$에 대하여,

$E(x) = E(-x),\, O(x) = -O(-x)$가 성립합니다. 

그렇다면 $P(x)$에서 홀수차수는 홀수차수끼리, 짝수차수는 짝수차수끼리 모아서 $E(x),\,O(x)$를 만든 다음,

$P(x) = E(x^2) + xO(x^2) $

$P(-x) = E(x^2) - xO(x^2) $

가 성립하기 때문에, $E(x^2)$의 값과 $O(x^2)$의 값을 아는것으로 계산을 단축할 수 있습니다. 벌써부터 분할정복의 느낌이 나죠?

위에서 설명했던 DFT로 돌아가자면, $w_{N}^{N/2} = -1$임을 이용하여 

$P(w_{n}) = E(w_{n}^2) + w_{n}O(w_{n}^2) $

$P(w_{n+N/2}) = E(w_{n}^2) - w_{n}O(w_{n}^2) $

을 알 수 있습니다. 

단 이 방식이 재귀적으로 적용되기 위해서는 다항식의 길이가 2의 거듭제곱 꼴 이어야합니다. 그럼으로 부족한 만큼 0을 추가해주면 됩니다.

그럼 우리는 FFT를 $\Theta (n\,log\, n)$에 계산할 수 있게 됩니다.

파이썬 코드를 살펴보자면,

def fft(p, w):
    n = len(p)
    if n == 1: #Base Case
        return p

    even, odd = p[0::2], p[1::2] #짝수차수와 홀수차수로 분할
    yeven, yodd = fft(even, w*w), fft(odd, w*w) #재귀적으로 불러옴 w*w인 이유는 w^2가 들어가기 떄문에

    z = complex(1, 0)
    y = [0] * n
    for j in range(n//2): #위에서 구한 공식
        y[j] = yeven[j] + z * yodd[j]
        y[j + n//2] = yeven[j] - z * yodd[j]
        z *= w
    return y

 

IDFT(Inverse Discrete Fourier Transform)

${a}_{n} = \frac{1}{N}\sum_{k=0}^{N-1}\hat a_{k}w_{N}^{-kn}$이 성립합니다. 이 식은 DFT와 굉장히 유사하게 생겼습니다. 사실 w대신 w의 역수를 대입하고, 모든 결과값을 N으로 나눠주면 됩니다.

이것이 왜 성립하는지 잘 이해가 가지 않을 수 있는데, 

$\begin{bmatrix}
 w_{N}^{0}& w_{N}^{0} &w_{N}^{0}  &w_{N}^{0}  \\
 w_{N}^{0}&  w_{N}^{1}&  w_{N}^{2}& w_{N}^{3} \\
 w_{N}^{0}&  w_{N}^{2}&  w_{N}^{4}& w_{N}^{6} \\
 w_{N}^{0}&  w_{N}^{3}& w_{N}^{6} & w_{N}^{9} \\
\end{bmatrix}
\begin{bmatrix}
a_{0} \\
 a_{1}\\
 a_{2}\\
a_{3} \\
\end{bmatrix}$

이 행렬곱은 N = 4일때 DFT의 정의와 정확히 일치하는것을 알 수 있습니다. 

그렇다면 이 행렬의 역행렬은,

$\frac{1}{N}
\begin{bmatrix}
 w_{N}^{0}& w_{N}^{0} &w_{N}^{0}  &w_{N}^{0}  \\
 w_{N}^{0}&  w_{N}^{-1}&  w_{N}^{-2}& w_{N}^{-3} \\
 w_{N}^{0}&  w_{N}^{-2}&  w_{N}^{-4}& w_{N}^{-6} \\
 w_{N}^{0}&  w_{N}^{-3}& w_{N}^{-6} & w_{N}^{-9} \\
\end{bmatrix}$

임이 자명합니다. 계산해보면 항등행렬이 나옵니다.

 

다항식 곱셈 파이썬 코드를 살펴봅시다.

def mult(a, b):
    n = 1
    while n <= len(a) and n <= len(b):
        n *= 2
    n *= 2 #다항식의 크기를 2의 거듭제곱꼴로 맞춰준다
    a = a + [0] * (n - len(a))
    b = b + [0] * (n - len(b))
    w = complex(math.cos(2*pi/n), math.sin(2*pi/n)) #w계산
    a = fft(a, w)
    b = fft(b, w) #a, b다항식을 각각 fft 취해준다.
    c = [0] * n

    for i in range(n):
        c[i] = a[i] * b[i] #fft의 결과값을 곱해준다.
    c = fft(c, 1/w) #IFFT를 취해준다.

    for i in range(n):
        c[i] /= n
        c[i] = round(c[i].real)
    return c

 

$\omega_{N} = e^{-2\pi i/N} = cos(-2\pi /N) + i \,sin(-2\pi / N)$ 이 성립함이 알려져 있기에, 우리는 $w$를 쉽게 계산할 수 있습니다.