고속 푸리에 변환(FFT) for ps

 

 

고속 푸리에 변환은 합성곱(convolution)을 빠르게 처리하기 위해 사용되는 알고리즘이다.

 

$$ F \ast G = \sum_k F[k] G[t-k] $$

 

합성곱은 이산 공간에서 위와 같이 정의된다. 변수의 합을 $t$로 고정시켜 놓고, $k$를 움직여 가며 누적한 결과라고 이해하면 되겠다. 가장 대표적인 상황은 두 다항식을 곱했을 때, 차수가 $t$인 항의 계수를 구하는 경우일 것이다. 합성곱을 정의 그대로 구현하면 당연히 $O(N^2)$이 소요되겠지만, fft를 사용하면 $O(N \log N)$만에 합성곱을 구할 수 있다.

 

 

이산 푸리에 변환

FFT를 이용해 합성곱을 구하는 과정은 다음과 같이 진행된다.

 

• 길이 $N$인 배열 $F, G$ 각각을 먼저 푸리에 변환($F_f , G_f$)한다. 변환 시 배열의 길이는 유지된다.
• 푸리에 변환에 대해 $F_f \ast G_f [j] = F_f[j]G_f[j]$가 성립하므로, 같은 인덱스의 값끼리 단순히 곱해 준다.
• 얻은 배열을 푸리에 역변환하여 합성곱을 구한다.

 

따라서 먼저 기본적인 푸리에 변환의 방법인 이산 푸리에 변환(DFT)에 대해 알아볼 필요가 있다. 이산 푸리에 변환은 (배열의 값들을 계수로 하는) 다항식에 특정 값을 집어 넣어 일종의 함숫값을 구하는 과정이다. 여기서 투입되는 특정 값은 $N$번 곱했을 때 1이 되는 값으로, 복소평면에서의 회전을 생각하면 $\exp(-i\frac{2\pi}{N}k)$이다.

 

$$ f(x) = F[0] + F[1]x + \cdots + F[N-1]x^{N-1} $$

$$ F_f[k] = f(\exp(-i\frac{2\pi}{N}k)) = \sum_{n=0}^{N-1} F[n] \exp(-i\frac{2\pi}{N}kn) $$

 

푸리에 변환은 다음과 같이 역변환이 가능하다. 식이 매우 유사한데, 기존 투입 값의 켤레복소수를 대입한 후에 $N$으로 나눠 주면 원래의 배열로 돌아온다.

 

$$ F[k] = \frac{1}{N} \sum_{n=0}^{N-1} F_f[n] \exp(i\frac{2\pi}{N}kn) $$

 

 

쿨리-튜키 알고리즘

분할 정복을 통해 $F_f[k]$를 빠르게 구할 수 있는데, 이 방법을 쿨리-튜키 알고리즘(Cooley-Tukey FFT)이라 한다. 편의상 이제부터 $\exp(-i\frac{2\pi}{N}k) = w^k_N$으로 표기한다. 다항식의 각 항들을 짝수 차수의 항과, 홀수 차수의 항끼리 묶어서 다음과 같이 표현할 수 있다. 길이 N은 단순화를 위해 2의 거듭제곱으로 가정하므로 짝수이다.

 

$$ f(w^k_N) = (F[0] + F[2] w^{2k}_N + \cdots + F[N-2] w^{(N-2)k}_N ) + (F[1] w^{k}_N + \cdots + F[N-1] w^{(N-1)k}_N ) $$

$$ F_f[k] = f(w^k_N) = f_{even} ({(w^k_{N})}^2) + w^k_N f_{odd} ( {(w^k_{N})}^2 )$$

 

추가로, 배열에서 다음 $N/2$번째 값은 복소평면 상에서 180도 회전한 상태이므로 $-w^k_N$을 대입한 값과 같다. 따라서 다음과 같이 나타낼 수 있으므로, $N/2$ 이전의 $k$에 대해서만 계산해도 충분하다.

 

$$ F_f[k+N/2] = f(-w^k_N) = f_{even} ({(w^k_{N})}^2) - w^k_N f_{odd} ({(w^k_{N})}^2 )$$

 

쿨리-튜키 알고리즘에 따라 배열의 홀수 항들과 짝수 항들의 변환 결과를 각각 먼저 구해 주고, 구한 값들을 바탕으로 푸리에 변환한 결과를 돌려 주면 되겠다. 코드는 아래와 같다.

def fft(a, w_n):
    n = len(a)
    if n==1: return a

    a_odd = fft(a[1::2], w_n**2)
    a_even = fft(a[0::2], w_n**2)

    res = [0]*n
    w = 1
    for k in range(n//2):
        res[k] = a_even[k] + w*a_odd[k]
        res[k+n//2] = a_even[k] - w*a_odd[k]
        w *= w_n
    return res

w_n = cmath.exp(complex(0, -2*cmath.pi/N))
inv_w_n = cmath.exp(complex(0, 2*cmath.pi/N))

ff = fft(f, w_n)
gf = fft(g, w_n)
r = [ff[i]*gf[i] for i in range(N)]
r = fft(r, inv_w_n)
for i in range(N): r[i] = round(r[i].real / N)