이산 푸리에 변환 (Discrete Fourier Transform)

Photo by Mitchell Y on Unsplash

글을 적는 현재 티스토리 에디터가 어색하게 느껴질 정도로 오랜만에 포스팅을 해 봅니다. 그동안 졸업 준비를 하느라 바빴습니다. 졸업 논문을 적고 심사 발표를 준비하는 일은 여간 힘든 일이 아니더군요. 다행히 예심은 잘 마쳤고, 본심까지 부족한 부분을 잘 보완할 일만 남았습니다. 그 사이 약간의 짬이 나서 글을 적습니다.

 

이번 포스팅의 주제는 이산 푸리에 변환(discrete Fourier transform, DFT)입니다. 약간 뜬금없게도 이 주제를 공부하게 된 계기는 쇼어의 알고리즘(Shor's algorithm) 때문입니다. 요새 학교에서 세미나 모임을 만들어 양자 컴퓨팅(qunatum computing)을 공부하고 있습니다. 여기서 여러 유명한 결과들을 살펴 보고 있는데, 그 중에서도 가장 유명한 결과는 단연 합성수를 소인수분해하는 쇼어의 알고리즘이라 생각합니다. 그런데 웬걸, 그 알고리즘을 제대로 이해하기 위해서는 푸리에 변환, 특히 DFT를 알아야 하더군요. 이름은 익히 들었지만 내용을 잘 몰랐기에 이번 기회에 한번 공부해 보자는 생각이 들었습니다. 공부하면서 저 스스로 이해한 내용을 공유해 봅니다.

허수와 복소수

DFT를 이해하기 위해서는 먼저 1의 거듭제곱근이 무엇인지 알아야 합니다. 그리고 1의 거듭제곱근이 무엇인지 알기 위해서는 허수 및 복소수의 세계에 대해서 공부해야 합니다.

 

허수 \(i\)는 고등학생 시절 처음 접했던 것으로 기억합니다. 그때는 무작정 무언가를 제곱해서 \(-1\)을 얻고 싶어서 "존재하지도 않는 수" \(i\)를 만들었다는 식으로 배웠던 것 같습니다. 이제야 느끼는 것이지만, 저러한 방식으로 허수 및 복소수를 가르치는 것은 그것들의 진정한 의미를 제대로 전달하지 못하는 처사라고 생각합니다.

 

이제부터 제가 따로 공부하면서 이해한 허수 및 복소수의 개념을 설명해 보고자 합니다. 사실 이를 위해서는 수론(number theory)이나 대수(algebra) 등 다양한 수학적 배경이 필요합니다만, 아무래도 저는 수학을 전공하지는 않았어서 그 이해가 조악합니다. 다만, 짧은 지식으로나마 저 스스로 직관적으로 이해한 방법을 소개합니다. 이 주제와 관련하여 더 좋은 자료들이 많으니 직접 찾아 보시는 것도 추천합니다. 특히 아래 영상은 제가 이 글을 작성할 때 도움이 된 영상입니다.

 

DMT PARK.  세상에서 가장 아름다운 수식을 이해해보자 (이과용). 2021.

우리가 익숙한 실수(real number)는 일차원적입니다. 모든 실수를 직선 위에 하나씩 순서를 지키며 옮겨 놓을 수 있죠. 이제 실수들의 곱셈을 생각해 보겠습니다. 아래의 예시는 우리가 상식적으로 알고 있는 내용들입니다. 편의 상 정수만 고려했습니다.

\[ 2 \times 3 = 6, \quad\quad 2 \times -3 = -6, \quad\quad -2 \times -3 = 6. \]

여기서 한 가지 의문이 생깁니다. \(2 \times 3\)은 \(2\)를 세 번 더한 것이니 \(6\)이 당연해 보입니다. 그렇다면 \(2 \times -3\)이나 \(-2 \times -3\)은 무슨 의미일까요? 전자는 순서를 바꾸면 \(-3\)을 두 번 더한 것과 같으니 \(-6\)이라고 억지를 부릴 수는 있겠습니다만, 후자는 그러한 꼼수도 통하지 않습니다.

 

정리하자면, 양수끼리의 곱셈은 직관적으로 받아들일 수 있겠지만, 곱셈에 음수가 참여하는 순간 우리의 상식만으로는 그 해석이 모호해지게 됩니다. 이 의문은 곧 \(-1\)의 제곱근이 무엇인지에 대한 질문으로도 직결됩니다. 동일한 실수를 두 번 곱했을 때 \(-1\)을 얻을 방법이 전무하기 때문입니다. 이러한 문제를 어떻게 해결할 수 있을까요?

 

이 모든 문제는 우리가 수를 일차원으로만 생각했기 때문에 발생한 것입니다. 하지만 수가 일차원이어야 한다는 법은 없습니다. 그 또한 고정관념일 뿐이죠. 따라서 기존 실수 직선 축과 직교하는 새로운 축을 생각해 보겠습니다. 그제야 우리는 동일한 무언가를 두 번 행하여 \(-1\)을 얻는 방법을 찾을 수 있습니다.

그림 1. -1의 제곱근을 얻는 방법

바로 \(1\)에서 출발해서 반시계 방향으로 \(90^\circ = \frac{\pi}{2}\) 회전하는 것이죠. 이때, 한 번 회전했을 때 도달하는 지점은 새로운 축 위에서 중심으로부터 \(1\)만큼 떨어진 곳입니다. 이를 기존의 축과 구분하기 위해 \(1i\) 혹은 그냥 \(i\)로 표현하겠습니다. 그리고 이것이 곧 우리가 익히 배운 허수 \(i\)의 의미가 되겠습니다. 자연스럽게 복소수는 기존의 실수 축과 새로운 "허수" 축으로 이루어진 좌표계에 대응합니다. 이를 복소 좌표계라고 부릅니다.

 

그렇다면, 이차원의 복소수는 음수가 참여하는 곱을 잘 설명할 수 있을까요? 그렇습니다. \(-1\)을 곱하는 것을 원점 대칭, 혹은 \(180^\circ = \pi\) 회전하는 것으로 생각하면 되기 때문입니다. 이는 음수와 음수의 곱이 양수가 되는 현상도 직관적으로 설명해 줍니다.

 

개인적으로 복소수를 이렇게 이해를 하고 나니 새로운 축에 "허황된 수", 허수라는 이름이 붙은 것이 약간은 불공평하다는 생각이 들었습니다. 불완전한 실수를 보완하여 주는 이 축 또한 실재하기 때문이죠. 그저 우리는 예전부터 실수의 세계가 익숙했고, 실수의 세계에서 허수의 세계를 상상하는 것은 "차원이 다른" 이야기였을 뿐인데 말입니다.

1의 거듭제곱근

허수 및 복소수의 개념을 이 정도까지 이해했으면 1의 거듭제곱근(roots of unity)도 직관적으로 이해할 수 있으리라 생각합니다. 어떠한 자연수 \(n\)에 대해, 1의 \(n\)차 거듭제곱근은 \[ x^n = 1 \]인 \(x\)들을 의미합니다. 실수의 세계만 생각을 한다면, \(n\)이 짝수일 때는 \(1\)과 \(-1\)이, 홀수일 때는 오직 \(1\)만이 위 식을 만족할 것입니다.

 

하지만 앞에서 우리는 허수를 도입하면서 허수가 관여하는 곱셈은 반시계 방향으로 회전시키는 것에 대응한다는 것을 보았습니다. 따라서 복소수의 세계에서는 더 많은 수들이 위 식을 만족하게 될 것을 예상할 수 있습니다. 그리고 정확히 복소 좌표계에서 중심이 원점이고 반지름이 1인 원을 그렸을 때, 실수 축의 1을 시작으로 원주를 정확히 \(n\) 등분 하였을 때의 각 지점들이 해당 식을 만족하게 됩니다.

그림 2. 1의 8차 거듭제곱근

1에서 반시계 방향으로 돌았을 때 처음으로 만나는 수를 특별히 \(\omega_n\)이라고 부르겠습니다. 경우에 따라 문맥 상 명확하면 \(n\)을 생략하고 \(\omega\)라고만 적을 수도 있습니다. 그러면 직관적으로 두 번째 수는 \(\omega\)가 이루는 각도로 두 번 회전했다는 의미이므로 \(\omega^2\)이 됨을 알 수 있습니다. 같은 이치로 모든 점은

\[\omega^0 = 1, \; \omega,  \; \omega^2, \; \ldots, \; \omega^{n - 1}\]

에 대응하게 됩니다.

 

이렇게 얻은 수들이 \(x^n = 1\)을 만족한다는 것을 직관적으로 보이겠습니다. 먼저, 원주를 \(n\) 등분 하였으므로

\[ \omega = \cos \frac{2\pi}{n} + i \sin \frac{2\pi}{n} \]

이 된다는 것을 확인하시기 바랍니다. 따라서 모든 \(j = 0, \ldots, n - 1\)에 대해,

\[ \omega^j = \cos \frac{2 \pi j}{n} + i \sin \frac{2\pi j}{n} \tag{1} \]

을 만족합니다. \(\omega^j\)를 \(n\) 제곱했다는 의미는 그것이 이루는 각도로 \(n\) 번 돌렸다는 뜻이므로

\[ (\omega^j)^n = \cos \left(\frac{2 \pi j}{n}\cdot n \right)  + i \sin \left( \frac{2\pi j}{n} \cdot n \right) = 1\]

이 된다는 것을 알 수 있습니다.

 

지수에 음수가 들어가는 경우도 생각해 봅시다. 일례로 \(\omega^{-1}\)을 고려해 보죠. 직관적으로는

\[ \omega \cdot \omega^{-1} = 1 \]

이 되어야 할 것입니다. 앞에서 \(\omega\)를 곱한다는 것은 반시계 방향으로 \(\frac{2\pi}{n}\) 돌리는 것에 대응한다고 하였습니다. 그렇다면 과연 \(\omega^{-1}\)을 곱하는 것을 어떻게 해석해야 다시 1로 돌아오게 만들 수 있을까요? 바로 반대 방향, 즉 시계 방향으로 회전시키는 것입니다. 즉, \(\omega^{-1}\)은 1에서 시계 방향으로 \(\frac{2\pi}{n}\) 만큼 돌린 수에 대응한다는 것을 알 수 있습니다. 잘 생각해 보면, 이는 \[ \omega^{-1} = \omega^{n - 1} \]과 같음을 유추할 수 있습니다.

그림 3. 1의 8차 거듭제곱근에서 \(\omega^{-1}\)의 의미.

위 논의를 확장하면, \[ \omega^{-1}, \; \omega^{-2}, \; \ldots, \; \omega^{-(n - 1)}  \]도 모두 1의 거듭제곱근이 된다는 것을 알 수 있습니다.

이산 푸리에 변환

이제부터 본론인 이산 푸리에 변환, 줄여서 DFT를 알아 보겠습니다. 참고로 벡터와 행렬을 표현할 때 0부터 시작하는 인덱스 방법을 따를 것입니다. 예를 들어, 어떤 \(n\) 차원 벡터 \(v\)의 각 원소는

\[ v = \begin{bmatrix} v_0 \\ v_1 \\ v_2 \\ \vdots \\ v_{n - 1} \end{bmatrix} \]

로 표현됩니다. 어떤 \(m \times n\) 행렬 또한 행은 0부터 \(m - 1\), 열은 0부터 \(n - 1\)의 인덱스를 갖습니다.

 

어떤 자연수 \(n\)에 대하여 \(n\)차 DFT는 다음과 같은 \(n \times n\) 행렬로 정의됩니다:

\[\mathcal{DFT}^{(n)} := \frac{1}{\sqrt{n}} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega & \omega^2 & \cdots  & \omega^{n-1} \\ 1 & \omega^2 & \omega^4 & \cdots & \omega^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{n-1} & \omega^{2(n-1)} & \cdots & \omega^{(n-1)(n-1)}\end{bmatrix}.\]

여기서 \(\omega := \omega_n\)입니다. 문맥 상 명확하면 \((n)\)을 제거하고 \(\mathcal{DFT}\)로만 적을 수도 있습니다. DFT가 행렬로 정의된 것에서 곧바로 알 수 있듯이, \(n\) 차원 복소수 벡터를 입력 받아 \(n\) 차원 복소수 벡터를 출력하는 함수 \(\mathcal{DFT}:\mathbb{C}^n \to \mathbb{C}^n\)으로도 생각할 수 있습니다.

역행렬

\(\mathcal{DFT}\)는 역행렬(inverse matrix)이 존재합니다. 이는 다음과 같습니다.^[ㄱ]

\[\mathcal{DFT}^\dagger := \frac{1}{\sqrt{n}} \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega^{-1} & \omega^{-2} & \cdots  & \omega^{-(n-1)} \\ 1 & \omega^{-2} & \omega^{-4} & \cdots & \omega^{-2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega^{-(n-1)} & \omega^{-2(n-1)} & \cdots & \omega^{-(n-1)(n-1)}\end{bmatrix}.\]

이것이 역행렬이 됨을 증명해 봅시다. 우리는

\[ \mathcal{DFT}^\dagger \mathcal{DFT} = I^{(n)} \]

이 되는 것을 보이고 싶습니다.(여기서 \(I^{(n)}\)은 \(n \times n\) 단위 행렬입니다.) 어떤 \(j = 0, \ldots, n - 1\)와 \(k = 0, \ldots, n - 1\)에 대해, \(\mathcal{DFT}^\dagger \mathcal{DFT} \)의 \(j\) 행 \(k\) 열의 원소의 값은 \( \mathcal{DFT}^\dagger \)의 \(j\)행과 \(\mathcal{DFT}\)의 \(k\) 열의 곱으로 정해집니다. 이는 곧

\[ \frac{1}{n} \begin{bmatrix} 1 & \omega^{-j} & \cdots & \omega^{-j(n-1)} \end{bmatrix} \begin{bmatrix} 1 \\ \omega^{k} \\ \vdots \\ \omega^{k(n-1)} \end{bmatrix} = \frac{1}{n} \sum_{\ell = 0}^{n - 1} \omega^{(k - j)\ell} \]

이 된다는 것을 알 수 있습니다.

정리 1. 합의 결과

---

\(\displaystyle \sum_{\ell = 0}^{n - 1} \omega^{(k - j)\ell} = \begin{cases} n & \text{ if } j = k \\ 0 & \text{ if } j \neq k  \end{cases}\)

[증명] 만약 \(j = k\)라면 \(\omega^0 = 1\)이므로 곧장 정리가 성립함을 알 수 있습니다. 이제 \( j \neq k \)일 때를 생각해 보겠습니다. 이전 절에서 \(\omega^{k - j}\)도 1의 거듭제곱근이라는 사실을 알 수 있습니다. 즉,

\[ x^n - 1 = (x - 1)(x^{n - 1} + x^{n - 2} + \dots + 1) = 0 \]

의 해가 되어야 합니다. 이때, \(j \neq k\)이므로 \(\omega^{k - j} \neq 1\)입니다. 따라서 반드시 두 번째 항은

\[ \omega^{(k - j)(n-1)} + \omega^{(k - j)(n-2)} + \ldots + 1 = \sum_{\ell = 0}^{n-1} \omega^{(k-j)\ell} = 0 \]

이 되어야 함을 이끌어 낼 수 있습니다. ■

 

이 정리를 통해 \(\mathcal{DFT}^\dagger \mathcal{DFT}\)는 정확히 대각선에 위치한 원소만 1을 갖고, 나머지는 모두 0인 단위 행렬이 된다는 것을 얻을 수 있으며, 이로써 \( \mathcal{DFT}^\dagger \)가 \(\mathcal{DFT}\)의 역행렬이 된다는 것을 보이게 됩니다.

 

어떤 정사각 행렬이 역행렬을 갖는다는 것은 아주 중요한 성질입니다. \(\mathcal{DFT}\)를 정의역과 공역의 크기가 같은 함수로 생각할 때, 정의역 \(\mathbb{C}^n\)의 모든 벡터가 공역 \(\mathbb{C}^n\)의 모든 벡터와 일대일 대응(one-to-one correspondence, bijection)을 이룬다는 의미입니다. 좀더 풀어서 설명하자면, 두 개 이상의 서로 다른 벡터가 같은 함숫값을 가질 수 없으며, 공역의 모든 벡터들도 정의역에 정확히 하나의 원상(preimage)을 갖는다는 뜻입니다. 이 성질은 아래의 논의에서 매우 중요한 지위를 갖습니다.

주기와 주파수

DFT는 왜 중요할까요? 다양한 이유가 있겠지만, 그중에서도 가장 유명한 이유를 꼽자면, 어떤 주기(period)를 띠는 벡터를 주파수(frequency) 영역으로 옮겨 준다는 것입니다. 이를 예시와 함께 설명해 보겠습니다. 다음과 같은 80 차원 벡터 \(v\)를 생각해 보겠습니다. 편의 상 Python의 list 형식으로 표기하였습니다.

 

v = [2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6]

 

이를 잘 보면, [2, 0, 3, 5, 7, 4, 1, 6]의 여덟 개 원소가 정확히 반복되고 있습니다. 따라서 주기는 8이 됩니다. 이제 이 벡터에 대해 DFT를 진행해 보겠습니다. 그 결과는 복소수이므로, 각 원소의 절댓값^[ㄴ]을 찍어 보겠습니다.

그림 4. 주기가 8인 \(v\)의 DFT 적용 후 각 원소의 절댓값.

흥미롭게도 정확히 10의 배수에 해당하는 원소에서만 0보다 큰 절댓값을 갖는다는 것을 알 수 있습니다. 10은 무슨 의미일까요? 원래 \(v\)는 80 차원에서 8의 주기를 갖고 있었습니다. 주파수는 주어진 시간 동안 얼마나 자주 반복되는지를 나타내는 수치로 \(v\)의 주파수는 \(80/8 = 10\)으로 정의할 수 있습니다. 따라서 DFT를 수행하면 결과 벡터는 원래 벡터의 주파수의 배수에서 0보다 큰 절댓값을 갖는다고 유추해 볼 수 있습니다.

 

한 번 더 확인할 겸 같은 차원의 다른 주기를 갖는 벡터를 생각해 봅시다. 이번 \(v\)는 아래와 같습니다.

 

v = [2, 0, 3, 8, 5, 7, 9, 4, 1, 6,
     2, 0, 3, 8, 5, 7, 9, 4, 1, 6,
     2, 0, 3, 8, 5, 7, 9, 4, 1, 6,
     2, 0, 3, 8, 5, 7, 9, 4, 1, 6,
     2, 0, 3, 8, 5, 7, 9, 4, 1, 6,
     2, 0, 3, 8, 5, 7, 9, 4, 1, 6,
     2, 0, 3, 8, 5, 7, 9, 4, 1, 6,
     2, 0, 3, 8, 5, 7, 9, 4, 1, 6]

 

이는 [2, 0, 3, 8, 5, 7, 9, 4, 1, 6]의 열 원소가 반복되고 있으므로 10의 주기를 가지며, 자연스럽게 주파수는 8이 됩니다. DFT를 적용한 후의 결과를 그려 보면 아래와 같이 우리가 예상한 대로 나온다는 것을 알 수 있습니다.

그림 5. 주기가 10인 \(v\)의 DFT 적용 후 각 원소의 절댓값.

 

어떤 벡터가 완벽한 주기성을 갖는다는 것은 그 자체로 좋은 일입니다. 그리고 DFT는 이러한 벡터를 받아 그것의 주파수를 완벽히 복원해 냅니다. 하지만 언제나 세상은 그리 만만하지 않죠. 세상에는 주기를 띠는 것들로 넘쳐나지만 그렇다고 그 모두가 완벽한 주기성을 갖지는 않을 것입니다. DFT가 유용하려면 얼추 반복되는 벡터가 주어질 때, 그것의 주파수를 얼추 반환하여 주는 강건함(robustness)이 필요합니다. 이를 실험해 보기 위해서 이번에는 처음 주기가 8인 \(v\)에서 마지막 두 원소를 지워 보았습니다.

 

v = [2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4, 1, 6,
     2, 0, 3, 5, 7, 4]

 

더이상 \(v\)는 완벽한 주기를 갖지는 않습니다. 하지만 여전히 얼추 8의 주기를 갖고 있기는 하죠. 따라서 DFT를 적용한 후에도 얼추 10의 주파수를 얻기를 예상해 볼 수 있습니다. 실지로 아래 그림에서와 같이 10의 배수 언저리에서 상대적으로 큰 절댓값을 갖는다는 것을 확인할 수 있습니다.

그림 6. 주기가 약간 어그러진 \(v\)의 DFT 적용 후 각 원소의 절댓값.

푸리에 변환은 이러한 성질 덕분에 다양한 분야에서 널리 사용되고 있습니다. 특히 쇼어의 알고리즘에서도 이 성질은 매우 요긴하게 쓰입니다. 푸리에 변환의 보다 자세하고도 직관적인 설명은 아래 영상을 참고하시기 바랍니다.

 

3Blue1Brown. 푸리에 변환이 대체 뭘까요? 그려서 보여드리겠습니다. 2018.

다항식으로 해석하기

DFT를 이해하는 또 다른 좋은 방법은 다항식의 세계에서 해석하는 것입니다. 임의의 \(f \in \mathbb{C}^n\)를 \(x^j\)의 계수가 \(f_j\)인 다항식으로 연관지어 보겠습니다. 즉,

\[ f=\begin{bmatrix}f_0 \\ f_1 \\ \vdots\\ f_{n-1}\end{bmatrix} \iff f(x) = f_0 + f_1 x + \ldots + f_{n - 1} x^{n-1} = \sum_{\ell = 0}^{n - 1} f_\ell x^\ell \]

로 연관을 짓는 것이죠.

 

이제 \(\mathcal{DFT}\)를 적용해 봅시다. 계산해 보면, \( \mathcal{DFT} f \)의 \(j\) 인덱스 원소는

\[ [\mathcal{DFT}f]_j = \frac{1}{\sqrt{n}} \sum_{\ell = 0}^{n - 1} f_\ell \omega^{j \ell} = \frac{1}{\sqrt{n}} f(\omega^j) \]

가 된다는 것을 쉽게 알 수 있습니다. 다시 적으면,

\[ \mathcal{DFT} f = \frac{1}{\sqrt{n}} \begin{bmatrix} f(\omega^0) \\ f(\omega^1) \\ \vdots \\ f(\omega^{(n-1)})  \end{bmatrix} \]

과 같다는 것을 알 수 있습니다.

 

이를 통해 우리는 최대 차수가 \(n-1\)인 어떤 다항식의 계수를 표현하는 벡터를 \(\mathcal{DFT}\)에 넣어 변환하면 해당 다항식의 특정한 함숫값 \(n\)개를 얻을 수 있습니다. 이 자체는 사실 그리 놀라운 일은 아닙니다. 어떤 다항식의 계수들을 알고 있을 때 함숫값을 구하는 것은 그리 어려운 일이 아니니 말이죠.

 

하지만 DFT가 흥미로운 이유는 역연산이 가능하다는 데 있습니다. 즉, 특정한 함숫값 \(n\)개를 갖고 역행렬 \(\mathcal{DFT}^\dagger\)에 넣어 변환하면 해당 다항식의 계수를 곧바로 얻을 수 있다는 뜻입니다. 한번 \(f(1) = 10\), \(f(2) = -6\), \(f(3) = 7\), \(f(4) = 2\), \(f(5) = -8\)을 만족하는 4차 다항식 \(f\)를 계산해 보시기 바랍니다. 못할 이유도 없겠지만, 꽤 지난한 과정이 필요하다는 것을 알 수 있을 것입니다. 실지로, 임의의 함숫값 \(n\) 개에서 다항식을 얻는 것은 가우스 소거법(Gaussian elimination) 등을 통해 \(O(n^3)\) 시간에 해결이 가능하지만, DFT를 사용하면 기본적으로는 \(O(n^2)\)시간, (후술하겠지만) 더 잘 하면 \(O(n \log n)\) 시간에도 가능합니다. 물론 DFT를 사용하려면 \( f(\omega^0), f(\omega^1), \ldots, f(\omega^{n-1}) \)의 특정한 함숫값이 필요합니다.

합성곱

DFT를 다항식의 계수와 특정한 함숫값 사이의 변환으로 이해하면 DFT를 활용하여 쉽게 합성곱(convolution)을 구하는 방법을 얻을 수 있습니다. 합성곱은 다양한 분야에서 활용되고 있으며, 특히 신호 및 이미지/영상 처리 분야에서 요긴하게 쓰이는 개념입니다. 이와 관련하여 아래 영상을 시청하시는 것을 추천드립니다.

3Blue1Brown.  But what is a convolution?. 2022.

이제 합성곱을 정의해 보겠습니다. 어떤 \(m\) 차원 벡터 \(f \in \mathbb{C}^m\)과 \(n\) 차원 벡터 \(g \in \mathbb{C}^n\)에 대해, \(f\)와 \(g\)의 합성곱 \(f * g \in \mathbb{C}^{m + n -1}\)은 다음을 만족하는 \((m + n - 1)\) 차원의 벡터로 정의됩니다.

• 모든 \(j = 0, 1, \ldots, m + n - 2\)에 대해, \(f*g\)의 \(j\) 인덱스 원소는 \[ (f*g)_j = \sum_{k = 0}^{m + n - 2} f_k g_{j - k} \]이다. 이때 \(f_{-1}\), \(f_{m}\)이나 \(g_{-2}\), \(g_{n + 1}\) 등과 같이 원래 \(f\)와 \(g\)의 차원의 범위를 넘어서는 인덱스의 원소는 0으로 간주한다.

예를 들어, \(f \in \mathbb{C}^4\)이고 \(g \in \mathbb{C}^3\)라면, \(f*g \in \mathbb{C}^6\)는

\[\begin{array}{rl} (f*g)_0 & = f_0 g_0 \\ (f*g)_1 & = f_0 g_1 + f_1 g_0 \\ (f*g)_2 & = f_0 g_2 + f_1 g_1 + f_2 g_0 \\ (f*g)_3 & = f_1 g_2 + f_2 g_1 + f_3 g_0 \\ (f*g)_4 & = f_2 g_2 + f_3 g_1 \\ (f*g)_5 & = f_3 g_2 \end{array}\]

과 같습니다.

 

위 결과를 잘 살펴 봅시다. 앞에서 우리는 벡터 \(f\)와 \(g\)를 각각 \(m-1\)차 다항식 \(f(x)\)와 \(n-1\)차 다항식 \(g(x)\)의 계수에 대응시켰습니다. 그러면 \(f*g\)는 무엇에 대응할까요? 바로 두 다항식의 곱 \(f(x)g(x)\)의 계수에 대응합니다. 실지로 앞에서 봤던 예시에서, 3차 다항식 \(f(x)\)와 2차 다항식 \(g(x)\)의 곱 \(f(x)g(x)\)는

\[ f_0 g_0 + (f_0 g_1 + f_1 g_0) x + (f_0 g_2 + f_1 g_1 + f_2 g_0) x^2 + (f_1 g_2 + f_2 g_1 + f_3 g_0) x^3 + (f_2 g_2 + f_3 g_1) x^4 + f_3 g_3 x^5 \]

라는 5차 다항식이 되며, 각각의 계수는 \(f*g\)와 동일함을 확인할 수 있습니다.

 

따라서 \(f\)와 \(g\)의 합성곱 \(f*g\)를 얻기 위해서는 각각 다항식으로 표현했을 때 \(f(x)g(x)\)의 계수를 구하면 됩니다. \(f(x)g(x)\)의 최대 차수는 \(m + n - 2\)로, 총 \(m + n - 1\) 개의 함숫값을 안다면 \(f(x)g(x)\)를 특정할 수 있습니다. 앞에서 논의했던 대로, 아무 함숫값을 갖고 온다면 이를 빠르게 구하기 쉽지 않을 것입니다. 하지만 1의 \(m + n - 1\)차 거듭제곱근 \(\omega := \omega_{m + n - 1}\)에 대해,

\[ \begin{bmatrix} f(\omega^0)g(\omega^0) \\ f(\omega)g(\omega) \\ f(\omega^2) g(\omega^2) \\ \vdots \\ f(\omega^{m + n - 2})g(\omega^{m + n - 2}) \end{bmatrix} \]

를 어떤 방식으로든 구할 수 있다면, 이를 \([\mathcal{DFT}^{(m+n-1)}]^\dagger\)에 넣고 변환하여 \(f(x)g(x)\)의 계수를 구할 수 있을 것입니다.

 

그러면 해당 함숫값들은 어떻게 구할 수 있을까요? 어렵지 않습니다. 원래는 \(m\) 차원인 \(f\)에 \(n - 1\) 개의 0을 아래에 덧대어 줍니다. 같은 방식으로 \(g\)에는 \(m - 1\) 개의 0을 덧대어 줍니다. 그러면 \(f\)와 \(g\) 모두 \(m + n - 1\) 차원의 벡터이지만, 다항식으로 표현하여 보면 결국 기존의 \(f(x)\)와 \(g(x)\)와 동일합니다. 이후 각각에 대해 \(\mathcal{DFT}^{(m+n-1)}\)을 수행하여 준다면 우리는

\[ \frac{1}{\sqrt{m + n - 1}}\begin{bmatrix} f(\omega^0) \\ f(\omega) \\ f(\omega^2) \\ \vdots \\ f(\omega^{m + n - 2}) \end{bmatrix}, \quad \frac{1}{\sqrt{m + n - 1}}\begin{bmatrix} g(\omega^0)  \\ g(\omega) \\ g(\omega^2) \\ \vdots \\ g(\omega^{m + n - 2}) \end{bmatrix} \]

을 얻게 됩니다. 각 원소끼리 곱셈을 해 주면 우리가 원하는 벡터를 구할 수 있습니다.

 

위 논의를 통해 우리는 아래의 정리를 이끌어 낼 수 있습니다.

정리 2. DFT와 합성곱

---

어떤 \(f \in \mathbb{C}^m\), \(g \in \mathbb{C}^n\)에 대해,

\[ f*g = \sqrt{m + n - 1} \; \mathcal{DFT}^\dagger \; [ \mathcal{DFT}f \circ \mathcal{DFT}g ] \]

를 만족한다. 단, 위 식에서 \(f, g\)는 모두 \(m + n - 1\) 차원으로 0을 덧대어 확장된 벡터이며, \(\mathcal{DFT}\)와 \(\mathcal{DFT}^\dagger\) 모두 \(m + n - 1\)차 DFT의 행렬들이다. 같은 차원의 두 벡터 \(u, v\)에 대해 \(u \circ v\)는 각 원소끼리 곱하는 아다마르 곱(Hadamard product)이다.

고속 푸리에 변환

마지막으로 컴퓨터과학도로서 DFT의 효율성을 논해 보고자 합니다. 어떤 자연수 \(n\), \(f \in \mathbb{C}^n\)에 대하여, \( \mathcal{DFT} f \)는 얼마나 빠른 시간에 구할 수 있을까요? 가장 멍청하게는 \(O(n^2)\) 시간이 필요할 것입니다. 하지만 \(n\)이 2의 거듭제곱 꼴, 즉, 어떤 자연수 \(c \geq 1\)에 대해, \(n = 2^c\)인 경우라면 우리는 이보다 훨씬 빠른 속도인 \(O(n \log n)\) 시간에 \( \mathcal{DFT} f \)를 구할 수 있습니다. 이 방법을 고속 푸리에 변환(fast Fourier transform, FFT)이라고 일컫습니다. 참고로 아래의 논의를 확장하면 \( \mathcal{DFT}^\dagger f \)도 같은 시간에 계산할 수 있습니다.

 

FFT는 다항식 표현법으로 이해하는 것이 편합니다. 우리의 목표를 다시 적어 보겠습니다. 우리는 \(f(x) = \sum_{\ell = 0}^{n - 1} f_\ell x^\ell\)에 대해,

\[ \begin{array} {rl} f(\omega^0) & = \sum_{\ell = 0}^{n - 1} f_\ell (\omega^0)^\ell \\ f(\omega) & = \sum_{\ell = 0}^{n - 1} f_\ell \omega^\ell \\ f(\omega^2) & = \sum_{\ell = 0}^{n - 1} f_\ell (\omega^2)^\ell \\ & \vdots \\ f(\omega^{n - 1}) & = \sum_{\ell = 0}^{n - 1} f_\ell (\omega^{n - 1})^\ell \end{array} \]

을 모두 구해야 합니다. 여기서 \(\omega\)는 1의 \(n\)차 거듭제곱근 \(\omega_n\)입니다.

 

이 중에서 어떤 \(j\) 인덱스의 것을 갖고 오겠습니다.(단, \(j = 0, 1, \ldots, n - 1\)) 이를 풀어 적어 보면

\[ f(\omega^j) = f_0 + f_1 \omega^{j} + f_2 \omega^{2j} + \ldots + f_{n - 1} \omega^{(n-1)j} \]

이 됩니다. 위 식을 차수가 짝수인 부분과 홀수인 부분으로 쪼개 보겠습니다. 2의 거듭제곱 꼴이라는 가정에 의해 \(n\)은 반드시 짝수입니다.

\[ \begin{array}{rl} f(\omega^j) & = (f_0 + f_2 \omega^{2j} + \ldots + f_{n - 2} \omega^{(n-2)j}) + (f_1 \omega^{j} + f_3 \omega^{3j} + \ldots + f_{n - 1}\omega^{(n-1)j}) \\ & = (f_0 + f_2 \omega^{2j} + \ldots + f_{n - 2} \omega^{(n-2)j}) + \omega^j (f_1 + f_3 \omega^{2j} + \ldots + f_{n - 1} \omega^{(n -2)j}) \end{array} \tag{2} \]

두 번째 등식은 뒤의 항들에서 \(\omega^j\)를 빼낸 것입니다.

 

기존의 \(f\)에서 크기가 절반인 다음 두 벡터 \(f^\mathsf{even}, f^\mathsf{odd} \in \mathbb{C}^{n/2}\)를 만들어 봅시다.

\[ f^\mathsf{even} := \begin{bmatrix} f_0 \\ f_2 \\ \vdots \\ f_{n - 2} \end{bmatrix}, \quad f^\mathsf{odd} := \begin{bmatrix}  f_1 \\ f_3 \\ \vdots \\ f_{n - 1} \end{bmatrix}. \]

이들 각각도 다항식으로 나타내 보면,

\[  \begin{array}{rl} f^\mathsf{even}(x) & = f_0 + f_2 x + \ldots + f_{n - 2} x^{\frac{n}{2} - 1} \\ f^\mathsf{odd}(x) & = f_1 + f_3 x + \ldots + f_{n - 1} x^{\frac{n}{2} - 1} \end{array} \tag{3} \]

의 \(\frac{n}{2} - 1\)차 다항식이 된다는 것을 알 수 있습니다. 위의 식 (2)에 이를 적용하여 보면

\[ f(\omega^j) = f^\mathsf{even}(\omega^{2j}) + \omega^j f^\mathsf{odd}(\omega^{2j}) \tag{4} \]

와 동일하다는 것을 이끌어 낼 수 있습니다.

 

아래의 정리는 \(n\)차 거듭제곱근과 \(\frac{n}{2}\)차 거듭제곱근 사이의 관계를 조명합니다. 그림 1을 참고하세요.

정리 3. \(\omega^{2j}\)의 성질

---

어떤 짝수 \(n \geq 4\)와 \(j = 0, \ldots, n - 1\)에 대해, \[ \omega^{2j}_n = \omega^{(j \mod \frac{n}{2})}_{n/2} \]를 만족한다.

[증명] 식 (1)을 확장하여 \[ \omega^{2j}_n = \cos (\frac{2\pi}{n} \cdot 2j) + i \sin ( \frac{2\pi}{n} \cdot 2j ) = \cos (\frac{2 \pi}{n/2} j) + i \sin (\frac{2\pi}{n/2} j ) = \omega^j_{n/2} \]를 이끌어 낼 수 있습니다. 또한 \(j = \frac{n}{2} q + k \) 꼴로 나타내면, \[ \omega^j_{n/2} = \cos(2q\pi + \frac{2\pi}{n/2} k) + i \sin(2q\pi + \frac{2\pi}{n/2} k) = \omega^r_{n/2} \]가 됨을 바로 알 수 있습니다. ■

 

이 정리를 활용하여 식 (4)를 다시 적어 보면, \(k := j \mod \frac{n}{2}\)에 대해,

\[ f(\omega^j_n) = f^\mathsf{even}(\omega^k_{n/2}) + \omega^j_n f^\mathsf{odd}(\omega^k_{n/2}) \tag{5} \]

가 됩니다. 다시 말해, 우리가

\[ \{ f^\mathsf{even}(\omega^k_{n/2}) \mid k = 0, \ldots, \frac{n}{2} - 1 \} \cup \{ f^\mathsf{odd}(\omega^k_{n/2}) \mid k = 0, \ldots, \frac{n}{2} - 1\} \tag{6} \]

의 모든 값을 알고 있다면, 식 (5)를 각 \(j = 0, \ldots, n-1\)마다 적용하여 우리가 원하는 \(\{f(\omega^j_n) \mid j = 0, \ldots, n - 1\}\)을 모두 구할 수 있다는 뜻입니다.

 

그렇다면 식 (6)의 값들은 어떻게 구하면 될까요? 재귀적으로 해결하면 됩니다. 현재 우리가 해결하려는 문제를 구조화하면 다음과 같습니다.

• 어떤 2의 거듭제곱 \(n \geq 2\)에 대해, \(n - 1\)차 다항식 \(f(x)\)가 주어졌을 때 \( \{ f(\omega^j_n) \mid j = 0, \ldots, n-1 \} \)을 구하시오.

\(n\)은 2 이상의 2의 거듭제곱 꼴이라는 가정에 의해, \(n > 4\)인 경우에는 \(n/2\)도 2 이상의 2의 거듭제곱입니다. 따라서 식 (6)의 모든 값들은 아래의 두 입력

• \(n := n/2\), \(f := f^\mathsf{even}\)
• \(n := n/2\), \(f := f^\mathsf{odd}\)

을 재귀적으로 해결하면 구할 수 있게 됩니다. 

 

지금까지 논의한 내용을 정리하면 다음 알고리즘을 얻을 수 있습니다.

알고리즘 1. 고속 푸리에 변환(FFT)

---

입력: 어떤 2의 거듭제곱 수 \(n \geq 2\), \(f \in \mathbb{C}^n\) 혹은 \( n-1 \)차 다항식 \(f(x)\)

출력: \(\mathcal{DFT}f\) 혹은 \( \{ f(\omega^j_n) \mid j = 0, \ldots, n-1 \} \)

 

1. 만약 \(n = 2\)라면, \(f(1)\)과 \(f(-1)\)을 계산하여 출력한다.
2. 만약 \(n > 2\)라면, 아래를 수행한다.
   1. 식 (3)에 따라 \(\frac{n}{2} - 1\)차 다항식 \(f^\mathsf{even}(x)\)와 \(f^\mathsf{odd}(x)\)를 만든다.
   2. \(f^\mathsf{even}\)과 \(f^\mathsf{odd}\) 각각을 입력으로 재귀적으로 호출하여 \( \{f^\mathsf{even}(\omega^k_{n/2}) \mid k = 0, \ldots, \frac{n}{2} - 1 \}\)과 \( \{ f^\mathsf{odd}(\omega^k_{n/2}) \mid k = 0, \ldots, \frac{n}{2} - 1 \} \)을 구한다.
   3. 각 \(j = 0, \ldots, n - 1\)에 대해, 식 (5)에 따라 조합하고 이를 모두 출력한다.

이 알고리즘의 수행 시간을 분석해 보겠습니다.

정리 4. FFT 수행 시간

---

알고리즘 1은 \(O(n \log n)\) 시간에 동작한다.

[증명] 어떤 \(n \geq 2\)에 대해, 이 알고리즘이 \(n\) 차원 벡터 \(f\)를 입력받았을 때의 수행 시간을 \(T(n)\)이라고 부르겠습니다. 단계 1에서 \(T(2) = O(1)\)임을 알 수 있습니다. 이제 \( n > 2 \)의 경우를 생각해 봅시다. 단계 2-b를 제외하고는 모두 \(O(n)\) 시간에 수행이 가능합니다. 단계 2-b에서는 절반의 크기의 벡터 \(f^\mathsf{even}\)과 \(f^\mathsf{odd}\) 각각에 대해서 재귀적으로 알고리즘을 수행합니다. 따라서 우리는 다음의 점화식을 이끌어 낼 수 있습니다.

\[T(n) = \begin{cases} O(1), & \text{ if } n = 2, \\ 2 T(\frac{n}{2}) + O(n), & \text{ otherwise. } \end{cases}\]

이는 병합 정렬(merge sort)의 점화식과 동일하며, 같은 방식으로 \(T(n) = O(n \log n)\)임을 보일 수 있습니다. ■

마치며

이것으로 이산 푸리에 변환에 관한 포스팅을 모두 마칩니다. 서론에서도 밝힌 바와 같이 DFT를 공부하게 된 계기는 양자 컴퓨팅과 쇼어의 알고리즘 때문이었습니다. 푸리에 변환 자체를 언젠가는 공부해야지 생각만 하다가, 뜬금없지만 얼결에라도 공부를 하게 되어 기쁘게 생각합니다. 다른 분야에서는 푸리에 변환이 어떻게 활용되는지 궁금하기도 합니다.

 

FFT에서 우리는 \(n\)이 2의 거듭제곱 꼴이어야 한다는 가정을 했습니다. 만약 \(n\)이 2의 거듭제곱 꼴이 아니라면 재귀적으로 동작하는 중에 언젠가 홀수인 \(n\)을 처리해야 합니다. 문제는 \(n\)이 홀수 일 때는 \(\frac{n}{2}\)가 정수가 아니어서 정리 3이 더 이상 성립하지 않는다는 것입니다. 그러한 경우를 적절히 처리하는 방법도 알려져 있다고 합니다만, 제가 참고한 교과서에서는 해당 내용을 소개하지 않았습니다. 개인적으로 한번 풀어 보고, 포스팅함 직하면 정리해 보도록 하겠습니다.

 

양자 컴퓨팅을 열심히 공부하고 있습니다. 처음에는 진짜 하나도 이해가 되지 않았는데, 세미나도 하고 혼자서 공부도 하다 보니 어느새 기본적인 내용들은 조금씩 익숙해져 가는 것 같습니다. 다음에는 양자 컴퓨팅 주제로 포스팅을 하고자 합니다. 저희 학교에서 학생들끼리 모여 진행하는 세미나는 아래 링크의 사이트에 모두 아카이빙하고 있습니다. 궁금하신 분들은 아래 링크로 방문하시면 되겠습니다.

Yonsei CS Theory Student Group

 

읽어 주셔서 고맙습니다.

참고 문헌

[1] Jon Kleinberg and Eva Tardos. Algorithm design. Pearson Education, 2013.

 

[2] Michael Loceff. A course in quantum computing (for the community college). 2015.

각주

[ㄱ] \(\mathcal{DFT}^\dagger\)는 \(\mathcal{DFT}\)의 켤레 전치(conjugate transpose) 행렬입니다.

 

[ㄴ] 어떤 복소수 \(a + bi\)의 절댓값은 \(\sqrt{a^2 + b^2}\)으로 정의됩니다.