[나만 보는 정리노트] FFT 기초(내가 배워가는 글, 교육용 X)

!!주의!! 푸리에 변환을 거의 이해하지 못한 사람의 발버둥입니다. 전혀 정확하지 않습니다. 오류가 매우 많습니다. 남들한테 보여주려고 업로드한 글이 아닙니다.

1. 푸리에 급수
   $f(x)=\displaystyle\sum_{k=0}^{\infin}{\frac{<f(x), \cos kx>}{||\cos kx||^2}\cos kx+\frac{<f(x), \sin kx>}{||\sin kx||^2}\sin kx}$
   $<⋅, ⋅>$는 함수의 내적으로, $\int_{a}^{b}{f(x)\times g(x)dx}$로 나타낸다. 내적은 곧 유사도를 말하기도 한다. 위끝과 아래끝은 나타내지 않는데, 아마도 푸리에 급수에서는 한 주기, 푸리에 변환에서는 $+\infin$부터 $-\infin$인 것 같다.
   ∥⋅∥는 노름(norm)이다. $||v||:=|v\cdot v|$로 정의된다.

진동수 $\frac{k}{2\pi}$인 사인파와 코사인파의 성분의 합으로 $f(x)$를 나타낸다.
즉 기저가 되는 사인파와 코사인파의 진동수를 $k$와 함께 증가시켜가며 내적을 해서 각 진동수에 대한 성분을 원래 파형에서 추출한 뒤, 이들의 합으로 원함수를 풀어서 나타낸 것이 푸리에 급수이다.

달리 표현하면
$f(x)=\displaystyle\sum_{k=-\infin}^{\infin}{c_ke^{ikx}}$이다. 여기서 $c_k$는 진동수 $k$ 성분이 얼마나 들어있는지, 즉 진폭이 얼마인지를 나타낸다. $e^{ikx}$는 직교 기저이기 때문에 $e^{ikx}$들 간 $\cos \theta$ 값이 0이므로 내적이 0이고, 따라서 진동수들 간 간섭이 없다.

이를 달리 표현했을 때 $f(x)=\frac{1}{2\pi}\displaystyle\sum_{k=-\infin}^{\infin}{<f(x), \psi _k> \psi _k}$이다.
$2\pi$는 정규화를 위한 숫자이다. 정규화란 벡터를 단위벡터로 바꾸는 것이다. 수식으로 나타내면 $\frac{\overrightarrow v}{|v|}$이다.
$\psi _k$는 $e^{ikx}$를 나타낸 것이다.
$<f(x), \psi _k>$는 $c_k$에 대응한다.

2. 푸리에 변환
   이렇게 푸리에 급수를 이용하여 어떤 성분이 포함되어 있는지를 알 수 있다. 그러나 각 성분의 진동수에 따라서 원함수의 한 주기에서도 나타나는 횟수가 다르므로 정확한 진폭의 측정에 영향을 줄 수가 있다. 그렇기에 주기 $T$를 $\infin$로 보내서 모든 성분이 한 번만 나타나게 한다. 이렇게 하면 진동수가 진폭 측정에 영향을 미치지 않으므로 각 진동수에 대한 크기를 정확하게 알아낼 수 있다.

지금까지의 정보를 바탕으로 푸리에 변환식을 살펴보자.
$F(x)=\int_{-\infin}^{+\infin}{f(t)e^{-i2\pi xt}}dx$
인테그랄을 계산하는 것은 내적이다. 푸리에 급수에서 한 주기씩 계산하는 인테그랄의 범위를 무한대로 보내는 것이 주기를 무한대로 보낸 것을 의미한다.
$e^{-2\pi xt}$는 오일러 공식에 따라 사인함수와 코사인함수를 묶어서 표현한 것이다.
사인/코사인의 주기가 $\frac{1}{x}$이므로 진동수가 $x$인 경우를 의미한다. 즉 진동수가 $x$일 때 진동수가 똑같이 $x$인 $\sin, \cos$에 대하여 내적을 해 유사도를 구함으로써 $f(x)$에 진동수 $x$인 성분이 얼마나 포함되어 있는지를 구하는 것이다.

3. 이산 푸리에 변환
   하지만 푸리에 변환은 적분이 사용되기 때문에 연속적인 데이터에만 적용할 수 있다. 컴퓨터가 다루는 데이터는 모두 이산적이므로 이산 푸리에 변환을 사용해야 한다. 이산 푸리에 변환도 그냥 푸리에 변환과 본질은 똑같다. 다만 인테그랄이 시그마로 바뀌었을 뿐이다.

$F_k=\displaystyle\sum_{j=0}^{N-1}{f_je^{\frac{i2\pi jk}{N}}}$

데이터의 점들 간의 내적을 대체제로 사용하는 것이다.

.
..
...

그런데, DFT가 시간 대역을 주파수로 바꿔주기만 하는 것이라면 무슨 의미가 있을까?

행렬곱은 곱해줄 행과 열 간의 내적과 같다. DFT와 근본이 같은 것이다. 그렇기에 DFT를 행렬곱으로 나타낼 수 있는데, 이를 '푸리에 행렬'이라고 부르기도 한다.

DFT의 식을 $w=e^{-i\frac{2\pi}{N}}$를 이용해 간단히 풀어보면
$X[i]=x[0]w^0+x[1]w^i+x[2]w^{2i}+...+x[N-1]w^{i(n-1)}$이다.
이를 행렬로 나타낸다면

$\begin{bmatrix}X_i\\ \end{bmatrix}=\begin{bmatrix}1&w^i&w^{2i}&...&w^{i(n-1)}\\ \end{bmatrix}\begin{bmatrix}x_0 \\ x_1\\ x_2 \\ ...\\x_{N-1} \end{bmatrix}$
이다.

이 행렬은 $N\times N$의 크기일 것이므로, 행렬곱을 위해 곱셈을 $N^2$번 해야 한다. DFT로 시간 대역 데이터를 주파수 대역으로 변환하는 것은 $O(N^2)$의 시간복잡도를 가짐을 알 수 있다.

4. 고속 푸리에 변환
   $O(N^2)$... 불만족스러운 시간복잡도이다. 이를 빠르게 하기 위해서 FFT가 개발되었다. 가장 많이 알려진 쿨리-튜키 알고리즘을 바탕으로 서술하겠다.

다항식을 표현하기 위한 대표적인 방법으로는 계수 표기법(Coefficient Representation)과 값 표기법(Value Representation)이 있다.

계수 표기법은 아주 간단하다. n차다항식의 계수 $n+1$개를 수열로 나타내주는 것이다. $3x^2+2x+1$은 $[1,2,3]$ 따위로 나타내면 된다.

값 표기법은 이보다는 조금 복잡하지만, 어렵지 않다. n차다항식에는 상수항을 포함해서 $n+1$개의 항이 있을 것이다. 따라서 $n+1$개의 지나는 점의 좌표를 알아낸다면 함수를 특정할 수 있다. 이를 라그랑주 보간법이라고 한다. 이 라그랑주 보간법에서 착안하여 n차다항식을 $n+1$개의 순서쌍들로 표현하는 방식이 값 표기법이다. $(x-1)^2$의 경우 $[(1, 0),(0,1),(2,1)]$로 나타낼 수 있을 것이다.

만약 우리가 $C(x)=A(x)\times B(x)$를 계산하고 싶다고 하자. 계수 표기법으로 할 경우 당연히 $O(N^2)$의 시간복잡도가 될 것이다. 만약 값 표기법으로 한다면 어떨까?

마찬가지로 $x_0,x_1,x_2,x_3,...,x_{n}$에 대해 $n+1$번의 곱셈을 해줘야 하므로 $O(N^2)$이다. 하지만 이를 $O(N\log N)$으로 줄일 수 있는 방법이 있다.

알다시피, 홀수차함수는 기함수이고 짝수차함수는 우함수이다. 기함수와 우함수는 그 특성상 $x_i$의 값을 구하면 $-x_i$의 값도 구할 수 있다. 그리고 계수 표기법으로 된 다항식은 홀수차함수와 짝수차함수로 나누기 아주 쉽다.

이렇게 나눠진 다항식은
$P(x)=P_e(x^2)+xP_o(x^2)$처럼 나타내면 될 것이다.

그런데, 각각의 $P_e$와 $P_o$에 대해서도 짝수차와 홀수차로 분할할 수 있을 것이다.(분할 정복)

이걸 계속 반복한다... 그러면 최종적으로 $\log N$번 분할하게 될 것이다.
분할해서 나온 결과들은 계속해서 $±$의 관계이므로 곱해서 취합하고, 다시 취합하고, 다시 취합하면 된다.

그런데, 분할한 부분을 곱해가며 취합하는 과정에서 계속해서 $±$ 쌍을 만들기 위해 음의 제곱근이 유도되는 부분이 있다. 이를 어떻게 해야할까?
앞서 나온 $w=e^{-i\frac{2\pi}{N}}$를 다시 살펴보자.

오일러 공식에 따라 $e^{-2\pi i}=1$이다. 즉 $w^n=1$이다. 이런 숫자를 1의 거듭제곱근(Root of unity)이라고 부르며, $n$제곱을 하면 1로 돌아오는 주기성을 가지고 있다.

이런 성질 덕분에 행렬에 $w$를 넣는다면 중간에 어떤 연산을 거치든 최종적으로 return할 때 1로 돌아오게 된다.

FFT의 역변환인 IFFT같은 경우에는 $e^{i\theta}$에 의한 회전 방향이 반대이므로 $w$ 대신 $w^{-1}$을 사용하고, 데이터 길이가 $N$임을 감안해 $N$으로 나눠줘서 할 수 있다.

자, 이제 충분히 만족스러운 $O(NlogN)$ 안에 FFT와 IFFT 모두를 할 수 있게 되었다. 그렇다면 공간복잡도는 어떨까? 통상적인 방법을 사용한다면 다항식을 분할할 때마다 배열을 재지정해야 하므로 메모리 사용량이 상당히 클 것이다.

만약 처음부터 배열이 $[0, 4, 2, 6, 1, 5, 3, 7]$처럼 정렬되어 있다면 굳이 배열을 재지정하지 않고 재사용할 수 있을 것이다.
이 때 정렬되는 인덱스는 이진수의 역순 배열이다. $(100\rightarrow 001)$

이 방법으로 공간복잡도 역시 만족스럽게 줄일 수 있다.

...물론 복소수 연산을 해야 한다는 크나큰 문제가 남아있다.

가장 간단하고 직관적이고 비효율적인 재귀 코드

import cmath

def fft(x, inv=False):
    N = len(x)
    if N <= 1:
        return x
    even = fft([x[i] for i in range(0, N, 2)])
    odd = fft([x[i] for i in range(1, N, 2)])
    T = [cmath.exp(complex(0, (2 if inv else -2) * cmath.pi * k / N)) * odd[k] for k in range(N // 2)] # odd에 x값, 즉 w^k를 곱해준다.
    rst = [even[k] + T[k] for k in range(N // 2)] + [even[k] - T[k] for k in range(N // 2)]
    if inv:
        rst = [round((i/N).real) for i in rst] # IFFT의 경우 N으로 나눈다.
    return rst