푸리에 분석 - 고속 푸리에 변환 증명하기

1. 푸리에 분석 - 왜 할까요
2. 푸리에 분석 - 고속 푸리에 변환 증명하기
3. 푸리에 분석 - 고속 푸리에 변환 구현하기

이 페이지는 최병선 교수님의 Wavelet 해석을 바탕으로 작성되었습니다. 그러나 책과는 달리 해석적으로 엄밀한 내용을 추구하지는 않습니다. 이산 푸리에 변환의 구현은 여기를 참고 바랍니다.

푸리에 변환

함수의 푸리에 급수는 다음과 같이 정의됩니다. \begin{align} f(t) = \sum_{k=0}^{\infty} a_k cos(k \pi t) + b_k sin(k \pi t) \cdots (1) \end{align}

그런데 사실 사인과 코사인 함수로 따로따로 전개하는 건 좀 비효율적이긴 합니다. 왜냐하면 오일러의 등식인 $e^{ix} = \cos x + i \sin x$에 의해 사인과 코사인 함수를 하나로 묶을 수 있거든요. 그래서 푸리에 급수를 지수함수의 형태로 나타내면 다음과 같습니다. \begin{align} f(t) &= \sum_{k=-\infty}^{\infty} \hat{f}(k) e^{k i \pi t} \cdots (2) \\ \hat{f}(k) &= \frac{1}{2}\int_{-1}^{1} f(t) e^{-k i \pi t} \cdots (3) \\ \end{align} 표기는 깔끔해졌지만, 대신 복소수가 난입하게 되었습니다.

이산신호의 푸리에 변환

해석학 교과서나 공학수학 교과에서 만나는 푸리에 변환의 문제는 보통 연속함수에 대한 것입니다. 그 중에서도 성질이 좋은 함수들뿐입니다. 하지만 우리가 컴퓨터를 통해 마주하게 되는 자료는 대부분 이산적인 자료입니다. 즉, 값이 띄엄띄엄 떨어져 있다는 것입니다. 여기에서 값이 띄엄띄엄 떨어져 있다는 것은 수학의 관점에서는 하나, 둘, 셋, … 이렇게 셀 수 있다는 뜻이라고 봐도 됩니다. 그러니까, 두 값이 아무리 가까이 붙어 있어도 하나하나 셀 수 있다면 이산적이고, 셀 수 없다면 연속적이라고 생각하면 됩니다(예를 들어서, 0 다음의 실수를 생각할 수 있을까요?).

이산적인 데이터 역시 푸리에 급수를 생각할 수 있습니다. 이산 푸리에 변환의 목적 역시, 이산 신호를 주기함수의 합으로 보자는 것입니다. 달리 생각하면, 삼각함수를 이용해 이산 신호를 보간하는 것이라고 볼 수도 있겠네요. 어떤 이산신호가 벡터 $v = [v_0, \cdots , v_{N-1}]^T(\in C^N)$으로 주어졌다고 가정합시다. 그러면 벡터 $v$의 이산 푸리에 변환 $\hat{v} = [\hat{v_0}, \cdots , \hat{v_{N-1}}]^T(\in C^N)$은 다음과 같이 정의됩니다.

\begin{align} \hat{v_n} = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1}v_k e^{-\frac{2 \pi i n}{N} k} \cdots (4) \end{align}

여기에서 $\hat{v}_n$은 $\hat{v}$의 $n$번째 원소를 의미합니다. 식 (3)은 행렬로 표기할 수 있는데, 이는 다음과 같습니다. \begin{align} \hat{v} = Fv \end{align} 여기에서 F는 다음과 같이 정의됩니다. \begin{align} \omega_{n} &= e^{-\frac{2 \pi i n}{N} k} (n=0,1,\cdots,N-1) \end{align} \begin{align} F = \frac{1}{\sqrt{N}} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_1 & \omega_1^2 & \cdots & \omega_1^{N-1} \\ 1 & \omega_2 & \omega_2^2 & \cdots & \omega_2^{N-1} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & \omega_{N-1} & \omega_{N-1}^2 & \cdots & \omega_{N-1}^{N-1} \\ \end{bmatrix} \end{align}

거꾸로, 푸리에 변환을 알고 있는 상태에서 원래 신호를 복원하려면 $F^{-1}$를 구하면 됩니다. $\hat{v} = Fv$이므로 $v = F^{-1}\hat{v}$거든요. 그런데 역행렬을 계산하는 건 매우 시간이 오래 걸리고 귀찮은 작업입니다. 다행히도, 이산 푸리에 변환을 구하는 행렬은 매우 좋은 성질을 가지고 있어서, 다음과 같이 계산됩니다. \[ FF^{*} = I\] 쉽게 말해, $F$의 역행렬을 구하려면 $F$의 모든 원소의 허수 부분의 부호만 반대로 바꿔주면 됩니다!

고속 푸리에 변환

이렇게 구한 이산 푸리에 변환은 계산이 오래 걸립니다. 복소수를 한 번 곱하고 더하는 것을 하나의 연산이라고 한다면, 우리는 총 $N^2$번의 계산을 해야 합니다. 즉, $N$이 2배 커지면 계산량은 4배, 3배 커지면 9배, 10배 커지면 계산량은 100배로 늘어납니다. 이런 알고리즘을 흔히 $O(N^2)$이라고 합니다.

이건 너무 비효율적입니다. 그러다보니 어떻게 하면 더 빠르게 계산할 수 없을까 하는 연구가 이어졌고, 결국 우리의 선배들은 답을 찾았습니다. $O(N^2)$인 알고리즘을 $O(N \log N)$으로 낮추는 획기적인(!!) 알고리즘이 등장하였고, 이 알고리즘을 Cooley-Tukey 알고리즘이라고 부릅니다. 이렇게 이산 푸리에 변환을 빠르게 계산하는 알고리즘들을 우리는 고속 푸리에 변환(Fast Fourier Transform)이라고 부르고, 약자로는 FFT라고 많이들 씁니다. 여기에서는 Cooley-Tukey 알고리즘의 증명을 소개하도록 하겠습니다.

Divide and Conquer

고속 푸리에 변환의 핵심은 분할과 정복입니다. 즉, 복잡한 문제를 단순한 여러 개의 문제로 쪼개는 것입니다. 이산 푸리에 변환 행렬 $F$는 전혀 성기지 않은 행렬입니다. 행렬의 모든 칸에 숫자가 들어차 있고, 그 숫자들마저도 그렇게 작지도 않습니다. 그렇기 때문에 어느 숫자 하나 무시할 수 없죠. 그런데 이 행렬은 놀랍게도 여러 개의 성긴 행렬(Sparse Matrix)들로 분해할 수 있습니다. 성긴 행렬이라는 것은 행렬의 원소들의 대부분이 0으로 채워져 있는 행렬입니다. 만약 $F$가 성긴 행렬 여러 개로 분해될 수 있다면, 정공법으로 뚫는 것보다 훨씬 적은 계산만으로도 원래의 목적을 달성할 수 있습니다.

이하의 내용은 고속 푸리에 변환의 증명입니다. 선형대수학이나 해석학, 그리고 수학의 증명법에 익숙하지 않으신 분들은 넘어가셔도 고속 푸리에 변환을 응용하는 데에는 전혀 무리가 없습니다. 또한 아래의 증명에서는 표기 및 구현의 편의상 행렬과 벡터의 인덱스는 0부터 시작한다고 가정하겠습니다.

먼저, 다음 보조정리를 증명하면서 FFT의 증명을 시작하겠습니다.

Danielson-Lanczos Theorem

행렬 $P_M^e$, $P_M^o$를 차원이 $M \times 2M$인 행렬이라고 하고, $P_M^e = \left[ \delta_{2m,n}\right]$, $P_M^o = \left[ \delta_{2m+1,n}\right]$ 이라고 하자. 또한 $I_M$은 $M \times M$ 항등행렬(Identity Matrix)로 정의하자. 그리고 $M \times M$ 대각행렬 $\Omega_M$을 다음과 같이 정의하자. \begin{align} \Omega_M &= \left[ \omega_{jk} \right] \\ \end{align} \begin{align} \omega_{jk} &= \begin{cases} e^{\frac{-i \pi k}{M}} &, j = k \\ 0 &, j \neq k \\ \end{cases} \end{align}

다음과 같은 $2M \times 2M$ 행렬들을 정의하자. \begin{align} P_{2M} = \begin{bmatrix} .P_M^e \\ P_M^o \\ \end{bmatrix} \end{align}

\begin{align} E_{2M} = \begin{bmatrix} I_M & \Omega_M \\ I_M & -\Omega_M \\ \end{bmatrix} \end{align}

그러면 이산 푸리에 행렬 $F_{2M}$은 다음과 같이 분해된다. \[ F_{2M} = \frac{1}{\sqrt{2}} E_{2M} ( F_M \bigoplus F_M) P_{2M} \] 여기에서 \[ F_M \bigoplus F_M = \left[ \begin{array} .F_M & O \\ O & F_M \\ \end{array}\right]\] 이다.

Proof

만일 $n=0,1,\cdots,M-1$이면 벡터 $\sqrt{2M}F_{2M} \cdot v$의 $n$번째 원소는 다음과 같다: \begin{align} \left[ \sqrt{2M}F_{2M}v\right]_{n} &= \sum_{k=0}^{2M-1} v_k e^{\frac{-2 \pi i k n}{2M}} \\ &= \sum_{k=0}^{M-1} v_{2k} e^{\frac{-2 \pi i (2k) n}{2M}} + \sum_{k=0}^{M-1} v_{2k+1} e^{\frac{-2 \pi i (2k+1) n}{2M}} \\ &= \sum_{k=0}^{M-1} v_{2k} e^{\frac{-2 \pi i (2k) n}{2M}} + e^{\frac{-2 \pi i n}{2M}}\sum_{k=0}^{M-1} v_{2k+1} e^{\frac{-2 \pi i (2k) n}{2M}} \\ &= \sum_{k=0}^{M-1} v_{2k} e^{\frac{-2 \pi i k n}{M}} + e^{\frac{-\pi i n}{M}} \sum_{k=0}^{M-1} v_{2k+1} e^{\frac{-2 \pi i k n}{M}} \\ &= \sqrt{M} [ F_{M} v^{e}]_{n} + \sqrt{M} e^{\frac{-\pi i n}{M}} [ F_{M} v^{o}]_{n} \\ \end{align} 여기에서 $v^{e}=[v_0, v_2, \cdots, v_{2M-2}]$, $v^{o}=[v_1, v_3, \cdots, v_{2M-1}]$ 이다. 또한 $v^{e} = P_{M}^{e} v$, $v^{o} = P_{M}^{o} v$이다.

만일 $n=M,M+1,\cdots,2M-1$이면, 다음 등식이 성립한다: \begin{align} e^{\frac{-2 \pi i (2k) n}{2M}} &= e^{\frac{-2 \pi i}{M} k (n - M)} \\ e^{\frac{-2 \pi i (2k+1) n}{2M}} &= e^{\frac{-\pi i}{M}(n - M)} \cdot e^{\frac{-2 \pi i}{M} k (n - M)} \\ \end{align} \[ (\because e^{\frac{-2 \pi i (2k) n}{2M}} \cdot 1 = e^{\frac{-2 \pi i}{M} \cdot k n} \cdot e^{\frac{2 \pi i}{M} \cdot k M} )\] $\therefore$ 만일 $n=M,M+1,\cdots,2M-1$이면 \[ [\sqrt{2M} F_{2M} v]_{n} = \sqrt{M} [F_{M} v^e]_{n-M} - \sqrt{M} e^{\frac{-\pi i}{M} (n - M)} [F_{M} v^o]_{n-M} \] \begin{align} (\because [\sqrt{2M} F_{2M} v]_{n} &= \sum_{k=0}^{M-1} v_{2k} e^{\frac{-2 \pi i (2k) n}{2M}} + \sum_{k=0}^{M-1} v_{2k+1} e^{\frac{-2 \pi i (2k+1) n}{2M}} \\ &= \sum_{k=0}^{M-1} v_{2k} e^{\frac{-2 \pi i}{M} k (n-M)} - e^{\frac{-\pi i}{M} k (n-M)} \sum_{k=0}^{M-1} v_{2k+1} e^{\frac{-2 \pi i}{M} k (n-M)} \\ &= \sqrt{M} [F_{M} v^e]_{n-M} - \sqrt{M} e^{\frac{-\pi i}{M} (n - M)} [F_{M} v^o]_{n-M}) \end{align}

\begin{align} \therefore [F_{2M} v]_{n} = \begin{cases} [\frac{1}{\sqrt{2}} (F_{M} P_{M}^{e} + \Omega_{M} F_{M} P_{M}^{o}) v]_{n}, &0 \leq n < M \\ [\frac{1}{\sqrt{2}} (F_{M} P_{M}^{e} - \Omega_{M} F_{M} P_{M}^{o}) v]_{n}, &M \leq n < 2M-1 \\ \end{cases} \end{align}

\begin{align} \therefore F_{2M} &= \frac{1}{\sqrt{2}} \begin{bmatrix} F_{M} P_{M}^{e} + \Omega_{M} F_{M} P_{M}^{o} \\ F_{M} P_{M}^{e} - \Omega_{M} F_{M} P_{M}^{o} \\ \end{bmatrix} \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} F_{M} & \Omega_{M} F_{M} \\ F_{M} & -\Omega_{M} F_{M} \\ \end{bmatrix} \begin{bmatrix} P_{M}^{e} \\ P_{M}^{o} \\ \end{bmatrix} \\ &= \frac{1}{\sqrt{2}} \begin{bmatrix} I_{M} & \Omega_{M} \\ I_{M} & -\Omega_{M} \\ \end{bmatrix} \begin{bmatrix} F_{M} & O \\ O & F_{M} \\ \end{bmatrix} \begin{bmatrix} P_{M}^{e} \\ P_{M}^{o} \\ \end{bmatrix} \end{align}

\[ \therefore F_{2M} = \frac{1}{\sqrt{2}} E_{2M} (F_{M} \bigoplus F_{M}) P_{2M}\]

FFT

이제 드디어 고속 푸리에 변환의 증명입니다. 위에서 고속 푸리에 변환은 Divide and Conquer, 즉 복잡한 행렬을 단순한 행렬의 곱으로 표현하는 것이라고 했습니다. 다음 명제가 바로 이 내용입니다.

만일 자연수 $q \in \mathbb{N}$에 대하여 $N = 2^{q}$이면 이산 푸리에 변환 행렬 $F_N$은 다음과 같다: \[ F_N = \frac{1}{\sqrt{N}} G_0^{N} G_1^{N} \cdots G_{q-1}^{N} Q_N \] 여기에서 $Q_N$은 $N \times N$ 순열 행렬(Permutation Matrix)로 다음과 같이 정의된다: \begin{align} & Q_{2} = I_{2} \\ & Q_{2M} = (Q_{M} \bigoplus Q_{M}) P_{2M}, \;\; M = 2^{1}, 2^{2}, \cdots , 2^{q-1} \\ \end{align} 또한 $G_{k}^{N}$은 $N \times N$ 성긴 행렬(Sparse Matrix)로 다음과 같이 정의된다: \[ G_{k}^{N} = \underbrace{E_{N/2^{k}} \bigoplus E_{N/2^{k}} \bigoplus \cdots \bigoplus E_{N/2^{k}}}_{2^{k}\text{개}}\]

Proof

수학적 귀납법으로 증명한다. 만일 $N = 2$라면, Danielson-Lanczos Theorem으로부터 다음 등식이 성립한다. \[ F_2 = \frac{1}{\sqrt{2}} G_0 P_2 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \\ \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ \end{bmatrix} \] \[(\because F_2 = \frac{1}{\sqrt{2}} E_2 (F_1 \bigoplus F_1) P_2 = \frac{1}{\sqrt{2}} \begin{bmatrix} I_1 & \Omega_1 \\ I_1 & -\Omega_1 \\ \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & 1 \\ \end{bmatrix} P_2\]

위 명제가 $N = 2^1,2^2,\cdots,2^q$$에 대하여 성립한다고 가정하자. 그러면 $N = 2^{q+1}$에 대하여 이 명제가 성립하는 것을 보이면 된다. 먼저, 다음 두 가지 등식이 성립함을 보이자.

• $G_{k+1}^{2M} = G_{k}^{M} \bigoplus G_{k}^{M}$ \begin{align} (\because G_{k+1}^{2M} &= \underbrace{E_{2M/2^{k+1}} \bigoplus \cdots \bigoplus E_{2M/2^{k+1}}}_{2^{k+1}개} \\ &= \underbrace{E_{M/2^{k}} \bigoplus \cdots \bigoplus E_{M/2^{k}}}_{2^{k+1}개} \\ &= \underbrace{\left(E_{M/2^{k}} \bigoplus \cdots \bigoplus E_{M/2^{k}} \right)}_{2^{k}개} \bigoplus \underbrace{\left(E_{M/2^{k}} \bigoplus \cdots \bigoplus E_{M/2^{k}} \right)}_{2^{k}개} \\ &= G_{k}^{M} \bigoplus G_{k}^{M}) \end{align}
• $G_{0}^{N} = E_N$ \[(\because E_N = E_{N/2^{0}}, \; G_{0}^{N} = \underbrace{E_{N/2^{0}}}_{1개})\]

$\therefore$ $M = 2^{q}(\therefore N = 2M)$라 하면 Danielson-Lanczos Theorem으로부터 다음이 성립한다: \begin{align} F_{2M} &= \frac{1}{\sqrt{2}} E_{2M} (F_{M} \bigoplus F_{M}) P_{2M} \\ &= \frac{1}{\sqrt{2}} \cdot \frac{1}{\sqrt{M}} \cdot E_{2M} ((G_{0}^{M} G_{1}^{M} \cdots G_{q-1}^{M} Q_{M}) \bigoplus (G_{0}^{M} G_{1}^{M} \cdots G_{q-1}^{M} Q_{M})) P_{2M} \\ &= \frac{1}{\sqrt{2M}} E_{2M} ((G_{0}^{M} \bigoplus G_{0}^{M}) \cdots (G_{q-1}^{M} \bigoplus G_{q-1}^{M})(Q_{M} \bigoplus Q_{M})) P_{2M} \\ &= \frac{1}{\sqrt{2M}} E_{2M} G_{1}^{2M} \cdots G_{q}^{2M} (Q_{M} \bigoplus Q_{M})) P_{2M} \\ &= \frac{1}{\sqrt{2M}} G_{0}^{2M} G_{1}^{2M} \cdots G_{q}^{2M} (Q_{M} \bigoplus Q_{M})) P_{2M} \\ &= \frac{1}{\sqrt{2M}} G_{0}^{2M} G_{1}^{2M} \cdots G_{q}^{2M} Q_{2M} \\ \end{align} \[\therefore F_N = \frac{1}{\sqrt{N}} G_{0}^{N} G_{1}^{N} \cdots G_{q}^{N} Q_{N} \]

행렬의 의미

FFT의 증명 자체는 끝났습니다만, 이렇게 분해된 행렬이 대체 왜 빠른지는 조금 더 고민이 필요합니다. 먼저, 행렬 $G_{k}^{N}$을 보도록 하겠습니다. 위에서 정의한 것과 같이, $G_{k}^{N}$은 다음과 같이 정의됩니다. \[ G_{k}^{N} = \underbrace{E_{N/2^{k}} \bigoplus E_{N/2^{k}} \bigoplus \cdots \bigoplus E_{N/2^{k}}}_{2^{k}\text{개}}\] 결국 $G_{k}^{N}$을 알려면 $E_{N/2^{k}}$을 알아야 합니다. $E_{N/2^{k}}$은 다음과 같이 정의됩니다. \begin{align} E_{2M} = \begin{bmatrix} I_M & \Omega_M \\ I_M & -\Omega_M \\ \end{bmatrix} \end{align}

$E_4$일 때 한 번 보도록 하겠습니다. \begin{align} E_{4} = \begin{bmatrix} I_2 & \Omega_2 \\ I_2 & -\Omega_2 \\ \end{bmatrix} = \begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & \omega \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -\omega \\ \end{bmatrix} \end{align}

보시는 바와 같이, $E_4$의 대부분은 0이고, 각 행에는 0이 아닌 숫자가 2개씩밖에 없습니다. 이러한 관찰은 모든 $E_{2M}$에 대하여 성립합니다. 즉, 각 행에는 0이 아닌 숫자가 2개밖에 없습니다. 그렇기 때문에 행렬을 한 번 곱하는 것의 시간복잡도가 매우 낮아지게 됩니다. 벡터의 길이가 $N$이라고 하면, 행렬 $G_{k}^{N}$을 한 번 곱하는 것의 시간복잡도는 $O(N)$이 됩니다. 그리고 $G_{k}^{N}$의 $k$는 $k = 0, 1, \cdots , q-1$인데, $q = \log_{2} N$입니다. 그렇기 때문에 우리는 시간복잡도가 $O(N)$인 연산을 $\log_{2} N$ 번만 수행하면 됩니다.

그렇다면 $Q_{N}$은 어떤 행렬일까요?

비트를 뒤집자

$Q_{N}$은 순열 행렬이라고 했습니다. 순열 행렬은 벡터의 순서를 바꿔주는 행렬입니다. 예를 들어, $[a, b, c, d]^T$인 벡터가 있다고 하면, 순열 행렬을 통과하면 $[b, c, d, a]^T$가 될 수 있습니다($Q_N$이 이렇게 바꿔준다는 건 아닙니다).

$Q_{N}$은 순열 행렬 중에서도 아주 특별한 놈입니다. 이 녀석은 인덱스의 비트를 뒤집어버립니다. 주의할 것은, 벡터의 순서를 정반대로 뒤집는다는 뜻이 아닙니다. 벡터의 인덱스의 비트를 뒤집는다는 겁니다.

무슨 뜻인지 조금 더 알아보도록 하겠습니다. 예를 들어, 벡터의 길이가 8이라고 해보겠습니다. 그리고 벡터의 시작은 0번째부터 시작하도록 하겠습니다. 그러면 벡터는 다음과 같이 되어 있을 것입니다. \[ v = \begin{bmatrix} (0번째) \\ (1번째) \\ (2번째) \\ (3번째) \\ (4번째) \\ (5번째) \\ (6번째) \\ (7번째) \\ \end{bmatrix}\]

벡터의 인덱스를 비트로 표현하겠습니다. 비트로 표현한다는 것은 2진법으로 표기한다는 것과 같은 뜻입니다. \[ v = \begin{bmatrix} (000번째) \\ (001번째) \\ (010번째) \\ (011번째) \\ (100번째) \\ (101번째) \\ (110번째) \\ (111번째) \\ \end{bmatrix}\]

이제 $v$에 행렬 $Q$를 곱해주면 다음과 같이 됩니다. \[ Q \cdot v = \begin{bmatrix} (000번째) \\ (100번째) \\ (010번째) \\ (110번째) \\ (001번째) \\ (101번째) \\ (011번째) \\ (111번째) \\ \end{bmatrix} = \begin{bmatrix} (0번째) \\ (4번째) \\ (2번째) \\ (6번째) \\ (1번째) \\ (5번째) \\ (3번째) \\ (7번째) \\ \end{bmatrix}\]

벡터의 원소들의 순서가 바뀌었습니다. 벡터의 순서가 뒤집어졌기 때문에, 이 행렬을 Bit-Inversal Involution Matrix라고 부릅니다. 그리고 이 행렬을 곱하는 것 역시 실질적으로는 $O(N)$의 시간복잡도를 나타냅니다. 그렇기 때문에 FFT(정확하게는 Cooley-Tukey 알고리즘)의 시간복잡도는 $O(N \log N)$입니다.

참고문헌

최병선, Wavelet해석, 세경사, 2001