FFT (Fast Fourier Transform)

개요

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

 

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{(}\mathbf{DFT}(A)\cdot \mathbf{DFT}(B))$ 의 결과값은 두 다항식을 곱한 다항식의 계수가 됩니다.

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

여기서 우리는 DFT를 $\mathrm{\Theta }(nlogn)$에 계산하는 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를 $\mathrm{\Theta }(nlogn)$에 계산할 수 있게 됩니다.

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

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으로 나눠주면 됩니다.

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

$\left[\begin{array}{cccc}{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{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\end{array}\right]$

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

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

$\frac{1}{N}\left[\begin{array}{cccc}{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{array}\right]$

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

 

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

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)+isin(-2\pi /N)$ 이 성립함이 알려져 있기에, 우리는 $w$를 쉽게 계산할 수 있습니다.