푸리에 급수
푸리에 변환에 대해 알아보기 전에 매우 친숙한 푸리에 급수에 대해서 먼저 알아볼 필요가 있다.
먼저, 푸리에 급수의 정의에 대해 알아보자.
그리고 이를 식으로 표현하면 주기가 T인 함수 f의 푸리에 급수를 다음과 같이 나타낼 수 있다.
\[\begin{align*} \lim_{N\rightarrow\infty}{f_N\left(t\right)} &= a_0+\sum_{n=1}^{\infty}\left(a_n\cos\frac{2\pi n}{T}t+b_n\sin\frac{2\pi n}{T}t\right) \\ &= a_0+\sum_{n=1}^{\infty}\left(a_n\cos nwt+b_n\sin nwt\right) \end{align*}\]이때 $a_n,b_n$ 은 푸리에 계수라고 정의한다. 그렇다면 주어진 복잡한 주기함수에서 푸리에 계수들은 어떻게 계산할 수 있을까?
핵심 아이디어는 적분이다. 삼각함수 sin,cos의 경우 한 주기에 대한 정적분의 값은 항상 0이 된다는 사실과, 주어진 함수의 한 주기를 잘라서 적분한 값은 선형 결합된 삼각함수를 모두 적분하여 합한 값과 같다는 사실을 이용하면 다음과 같은 식을 이끌어 낼 수 있다.
\[\int_{0}^{T}{f\left(t\right)dt=T\cdot a_0+0+\cdots}+0\]이로써 $a_0$의 값을 쉽게 계산해낼 수 있다.
\[a_0=\frac{1}{T}\int_{0}^{T}f\left(t\right)dt\]나머지 $a_n,b_n$의 값들도 같은 아이디어로 계산해낼 수 있다. 양변을 적분하여 구하고자 하는 상수만 포함된 식이 남게 만들면 되는 것이다.
식을 어떻게 변형하면 하나의 상수만 남기고 모두 제거할 수 있을까? 바로 각 항에 곱해진 삼각함수를 양변에 곱한 뒤 적분을 해주면 된다.
즉 $a_1$을 구하고 싶다면 원래 식에서 $a_1$에 곱해져 있는 $\cos wt$를 다시 양변에 곱해주는 것이다. 이 방법의 삼각함수의 곱셈 공식에 의한 것인데, $\sin$과 $\cos$에 관련된 삼각함수의 곱셈 공식은 다음과 같다.
\[sinA\cdot cosB=\frac{1}{2}(sin(A-B)+sin(A+B))\\ sinA\cdot sinB=\frac{1}{2}(cos(A-B)-cos(A+B))\\ cosA\cdot cosB=\frac{1}{2}(cos(A-B)+cos(A+B))\\\]이를 바탕으로 푸리에 급수 식에서 삼각함수를 곱한 뒤 양변을 적분한 결과를 살펴보자.
여기서도 중요한 포인트는 삼각함수의 경우 한 주기 내에서 적분을 한 값은 $0$이 된다는 것이다. 이를 잘 이용하면 $\sin$과 $\cos$ 함수가 곱해져 있는 경우에는 주파수가 같지 않거나, 또는 같더라도 $\sin 0=0$ 이므로 전체 적분 값은 언제나 $0$이 된다는 사실을 알 수 있다.
그렇다면 $\sin$과 $\sin$, $\cos$와 $\cos$ 함수가 곱해져 있는 형태는 어떨까? 만약 주파수가 다르다면 마찬가지로 전체 적분 값은 $0$이 될 것이고, 주파수가 같다면 $\cos(A-B)=1$이 되어 적분을 하게 되면 다음과 같은 값을 가지게 될 것이다.
\[\int_{0}^{T}{A\cdot \sin\ wt\cdot \sin\ wt\ dt=}\int_{0}^{T}{A\cdot\frac{1}{2}\left(1-cos\ 2wt\right)dt=}\frac{AT}{2}\]그렇기에 양변에 현재 구하고자 하는 상수에 곱해져 있는 삼각함수를 곱해주고 적분한다면 상수가 포함된 값만 남게 되는 것이다! 이를 수식으로 다시 정리한다면 다음과 같다.
\[\begin{split} \int_{0}^{T}{f\left(t\right)\times\cos xwt\ dt=}\int_{0}^{T}{\left[a_0\cos xwt+a_1\cos wt\times \cos xwt+\cdots \\ +a_x\cos xwt\times \cos wxt+\cdots+b_1\sin wt\times \cos xwt+\cdots\right]dt}=\frac{T}{2}\cdot a_x \end{split}\]따라서 우리는 $a_x$의 값을 다음과 같은 적분 식을 통해 구할 수 있다.
\[a_x=\frac{2}{T}\int_{0}^{T}{f\left(t\right)\times \cos xwt\ dt}\]마찬가지로, b_x의 값도 같은 방법으로 구할 수 있다.
\[b_x=\frac{2}{T}\int_{0}^{T}{f\left(t\right)\times \sin xwt\ dt}\]하지만 여기서 멈추지 않고 더 나아가서 복소 지수를 사용하여 식을 표현할 수도 있다. 바로 오일러 공식을 활용하면 된다.
\[e^{i\theta}=\cos\theta+i\cdot\sin\theta\ \ \ \ e^{-i\theta}=\cos\theta-i\cdot \sin\theta\]위 두 식을 더하고 뺌으로써 sin,\ cos 함수를 나타낼 수 있다.
\[\sin nwt=\frac{e^{inwt}-e^{-inwt}}{2i}\ \ \ \cos nwt=\frac{e^{inwt}+e^{-inwt}}{2}\]이를 본 식에 다시 대입하여 지수 함수를 기준으로 묶어준 결과는 다음과 같다.
\[a_0+\sum_{n=1}^{\infty}\left(a_n\frac{e^{inwt}+e^{-inwt}}{2}+b_n\frac{e^{inwt}-e^{-inwt}}{2i}\right)=a_0+\sum_{n=1}^{\infty}\left(\frac{1}{2}\left(a_n-ib_n\right)e^{inwt}+\frac{1}{2}\left(a_n+ib_n\right)e^{-inwt}\right)\]여기서 지수 함수의 계수들을 묶어 본다면
\[c_0=a_0,c_n=\frac{1}{2}\left(a_n-ib_n\right),c_{-n}=\frac{1}{2}\left(a_n+ib_n\right)\]라고 할 수 있다. 그리고 여기서 $a_n,b_n$을 미리 계산하였던 적분 꼴에서 삼각함수 부분을 복소 지수 함수 꼴로 치환해보자.
\[a_n=\frac{2}{T}\int_{0}^{T}{f\left(t\right)\times c o s\ nwt\ dt}=\frac{2}{T}\int_{0}^{T}f\left(t\right)\frac{1}{2}{(e}^{inwt}+e^{-inwt})dt=\frac{1}{T}\int_{0}^{T}{f}\left({t}\right){e}^{inwt}{dt}+\frac{1}{T}\int_{0}^{T}{f}\left({t}\right){e}^{-{inwt}}{dt}\] \[{ib}_n=\frac{2i}{T}\int_{0}^{T}{f\left(t\right)\times s i n\ nwt\ dt}=\frac{2i}{T}\int_{0}^{T}f\left(t\right)\frac{1}{2i}{(e}^{inwt}-e^{-inwt})dt=\frac{1}{T}\int_{0}^{T}{f}\left(t\right){e}^{inwt}{dt}-\frac{1}{T}\int_{0}^{T}{f}\left(t\right){e}^{-inwt}{dt}\]이를 다시 $c_n$에 대입한다면
\[c_n=\frac{1}{T}\int_{0}^{T}f\left(t\right)e^{-inwt}dt\]이 되고, 이를 다시 본 식에 대입한다면
\[f\left(t\right)=\sum_{n=-\infty}^{\infty}{c_ne^{inwt}}\]이와 같은 깔끔한 형태의 완성된 식을 얻을 수 있다.
그런데 이처럼 복소 지수를 활용하여 표현함으로써 가지는 의미는 무엇일까? 그 열쇠는 바로 오일러의 공식에 있다.
복소 평면 위에서 각도에 따른 원의 식은 $e^{i\vartheta}$로 나타낼 수 있는데, 푸리에 급수 식에 등장하는 $c_ne^{inwt}$는 반지름이 $c_n$이고, 각속도가 $nw$인 원들을 나타내는 것이다! 그리고 이 원들을 평면에 사영함으로써 파동의 형태를 표현할 수 있는 것이다.
이렇게 푸리에 급수에 대해서 어느정도 알아보았다. 그런데, 푸리에 급수의 한계는 그 정의에서 알 수 있듯이 바로 주기를 가지는 함수에만 적용이 가능하다는 것이다. 하지만 만약 주기가 없는 함수를 주기가 무한인 함수로 보는 참신한 아이디어를 바탕으로 확장한다면 이 또한 삼각함수의 급수 전개로 나타낼 수 있을 것이다. 그리고 이것이 바로 푸리에 변환이다!
푸리에 변환
그렇다면 푸리에 변환을 표현하기 위해 어떻게 식을 확장 시켜야 할까? 앞서 푸리에 급수를 통해 주어진 주기 함수를 다양한 주파수를 가진 삼각함수로 분해할 수 있었다.
이때 다양한 주파수는 주어진 함수의 기본 주파수의 배수들로 구성되어 있었다. 그런데 푸리에 변환은 비주기 함수들을 주기가 무한인 함수로 보는 것이고, 주파수는 주기의 역수이므로 기본 주파수는 0에 수렴하고, 이는 곧 분해된 함수들의 주기가 연속적으로 표현될 것임을 뜻한다.
그렇기에 푸리에 급수 식에서 주기를 무한대로 보냄으로써 푸리에 변환 식을 구할 수 있다!
\[f\left(x\right)=\int_{-\infty}^{\infty}f\left(t\right)e^{-iwt}dt\] \[f\left(t\right)=\int_{-\infty}^{\infty}f\left(x\right)e^{iwt}dx\]푸리에 변환의 또 다른 의미는 원래 시간 영역에서 정의된 함수를 주파수 함수로 변환하는 것이기에 첫번째 식을 푸리에 변환, 주파수 함수로부터 시간 영역의 함수를 복원하는 두번째 식을 푸리에 역변환이라고 한다.
이제 연속적인 시간에서 이산적인 시간으로 넘어갈 차례이다. 앞서 진행한 푸리에 변환을 연속적인 시간과 주파수를 기준으로 진행했다면 이산적인 시간, 또는 이산적인 주파수에서 푸리에 변환을 진행하는 것에 대해 알아볼 필요가 있다. 이는 현실적인 문제와도 연관 되어있다. 아무리 연속된 값이 주어진다 하더라도 유한한 성능을 지닌 컴퓨터가 처리할 때는 항상 유한한 시간으로 잘라서, 즉 이산적인 데이터로 변환하여 처리해야 하기 때문이다.
먼저 신호 시간이 이산시간으로 주어진 상황에서의 푸리에 변환 DTFT, Discrete-Time Fourier Transform 이산시간 푸리에 변환에 대해 알아보자. 이산시간이란 무엇일까?
이산의 정의에 따르면 시간과 시간 사이에 간격이 있음을 의미한다. 즉 어떠한 신호, 파장에 대해 시간 간격을 두면서 측정한 신호인 것이다. 그리고 이 과정을 샘플링이라고 부른다. 주어진 연속 신호 $f(t)$를 간격 $T$로 샘플링을 진행하면 $f(t)=f(nT)$, 즉 $n$에 관한 식으로 쓸 수 있고, 이로써 주파수는 불연속적인 값을 가지게 되므로 적분 대신 시그마를 사용하여 DTFT의 정의를 완성할 수 있다.
\[f(x)=\sum_{n=-\infty}^{\infty}{f(n)e^{-iwn}}\]DTFT의 역연산 또한 마찬가지로 정의할 수 있다.
\[f\left(n\right)=\int_{-\infty}^{\infty}f\left(x\right)e^{iwn}dx\]이때, 이산시간 주기 신호는 삼각함수의 성질에 의해 다음과 같은 식이 성립한다.
\[\cos wn=\cos 2\pi fn=\cos 2\pi\left(f+1\right)n=\cdots\]즉, 주파수는 정수만큼 변하면 같은 신호가 됨을 알 수 있다. 따라서 동일한 주파수에 대해 중복으로 더해지는 것을 막기 위해 적분 범위를 0부터 1로(주기가 1이므로 상한과 하한의 차이가 1이 나도록 범위를 잡으면 된다) 변환해 준다.
\[f\left(n\right)=\int_{0}^{1}f\left(x\right)e^{iwn}dx\]하지만 여전히 시간과 마찬가지로 주파수 또한 연속적인 값을 측정하여 계산할 수 없기에 주파수에 대해서도 연속이 아닌 이산 주파수의 형태로 변환하여 계산해야 한다. 그리고 이를 DFT Discrete Fourier Transform이라고 부른다.
주파수 $f$를 $\frac{1}{N}$ 간격으로 샘플링 하게 되면, 주기가 곧 $N$이 되고, 주파수가 $f=\frac{k}{N},k\in\mathbb{Z}$인 곳에서만 값을 가지게 된다. 따라서 DTFT 식을 변형하면 DFT의 식을 이끌어낼 수 있다.
\[f\left(k\right)=\sum_{n=0}^{N-1}{f\left(n\right)e^{-i\frac{2\pi}{N}kn},0\le k<n}\]그렇다면 여기서 거꾸로 $f\left(n\right)$은 어떻게 계산, 즉 IDFT의 식은 어떻게 구할 수 있을까? 핵심 아이디어는 내적이다. 두 주파수가 있을 때, 이 둘의 한 주기 구간 내에서 내적 값은 주파수가 같다면 그 주기가 적분 값이 되고, 주파수가 다르다면 0이 된다. 이를 이용하면 다음과 같은 식을 이끌어 낼 수 있다.
\[f\left(n\right)=\frac{1}{N}\sum_{k=0}^{N-1}{f\left(k\right)e^{i\frac{2\pi}{N}kn},0\le n<N}\]빠른 푸리에 변환
이제 최종 목표, 빠른 푸리에 변환에 도달했다. 앞서 설명한 것처럼 DFT 연산을 빠르게 처리하여 시간 복잡도를 $O\left(N^2\right)$에서 $O\left(NlogN\right)$으로 줄이는 마법과도 같은 알고리즘이다.
현재 알려진 모든 FFT 알고리즘의 시간 복잡도는 $O\left(NlogN\right)$이며, 이보다 작은 시간 복잡도가 불가능하다는 것에 대해서는 아직 증명이 되지 않았다. 이 글에서는 가장 많이 사용되고 있는 Cooley-Tukey FFT Algorithm에 대해서 알아보고자 한다.
CTFT의 핵심 아이디어는 분할 정복이다. N개의 데이터를 반으로 반복적으로 쪼개고 계산함으로써 계산량을 줄여 나가는 것이다. 이때, 데이터를 어떠한 방식으로 쪼개는 것이 가장 효율적으로 시간 복잡도를 단축시키는데 도움을 주게 될까?
바로 짝수항과 홀수항으로 나누는 것이다. 이로써 기존의 계산량이 $N^2$이었다면 반으로 나눈 뒤 FFT를 해주면 $\left(\frac{N}{2}\right)^2+\left(\frac{N}{2}\right)^2=\frac{N^2}{2}$ 로 절반이 되고, 이를 계속 수행하기 위해서는 $N$이 항상 $2^x$꼴이어야 한다는 조건이 붙는다. 이제 다시 이를 이용해서 DFT 식을 다시 적어보자.
\[\begin{align*} f\left(k\right) &=\sum_{n=0}^{N-1}{f\left(n\right)e^{-i\frac{2\pi}{N}kn}} \\ &=\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n\right)e^{-i\frac{2\pi}{N}k\cdot2n}+}\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n+1\right)e^{-i\frac{2\pi}{N}k\cdot\left(2n+1\right)}} \\ &=\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n\right)e^{-i\frac{2\pi}{\frac{N}{2}}kn}+}e^{-i\frac{2\pi}{N}k}\cdot\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n+1\right)e^{-i\frac{2\pi}{\frac{N}{2}}k\cdot n}} \\ &=E_k+e^{-i\frac{2\pi}{N}k}O_k \end{align*}\]즉, 짝수항에 대한 DFT와 홀수항에 대한 DFT 연산 결과를 더하는 식으로 표현할 수 있기에 $E_k,\ O_k$의 조합만으로 DFT 연산 결과를 만들어 낼 수 있는 것이다. 그런데 여기서 더 나아가 복소 지수의 주기성에 의해 계산량을 또 줄일 수 있다.
\[e^{-i\frac{2\pi}{N}(k+\frac{N}{2})}=e^{-i\frac{2\pi}{N}k}\cdot e^{-i\frac{2\pi}{N}\frac{N}{2}}=e^{-i\frac{2\pi}{N}k}\cdot e^{-i\pi}=-e^{-i\frac{2\pi}{N}k}\]위 식에 의해 다음이 성립하고,
\[f\left(k+\frac{N}{2}\right)=E_k-e^{-i\frac{2\pi}{N}k}O_k\]이로써 총 계산량이 절반이 됨으로써 전체 시간 복잡도 $O\left(\frac{N}{2}logN\right)$를 완성할 수 있다!
그런데, 이를 어떻게 큰 수의 빠른 곱셈에 활용할 수 있는 것일까? 여기에는 또 하나의 트릭이 필요하다! 우선 다항식의 곱셈에 FFT를 어떻게 활용할 수 있는 지부터 알아보자
빠른 다항식의 곱셈
우리는 다항식을 표현할 때 주로 계수와 차수를 이용해 표현한다.
\[f\left(x\right)=a_0+a_1x+a_2x^2+\cdots a_nx^n\]그리고 이를 우리는 Coefficient Representation이라고 부른다. 하지만 우리는 이 방법 대신 $n+1$개의 $\lbrace{x}_i,f\left(x_i\right)\rbrace$ 쌍을 이용하여 표현하는 방법, Point-Value Representation을 사용해 보고자 한다.
이를 사용하면 두 다항식 $f\left(x\right),g\left(x\right)$을 곱할 때 더 높은 차수(더 낮은 차수의 경우 계수가 $0$인 항들이 빈자리를 채워준다고 생각하면 된다)만큼 고른 적절한 임의의 x값들을 준비하고, $f\left(x_i\right)\times g\left(x_i\right)=h\left(x_i\right)$를 이용하여 두 다항식의 곱 $h\left(x\right)$를 구할 수 있다! 그런데 여기서 만약 이 $x$값들과 함숫값 사이의 변환을 빠르게 할 수 있다면 다항식의 곱을 매우 빠르게 계산할 수 있을 것이다.
\[f\left(k\right)=\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n\right)e^{-i\frac{2\pi}{\frac{N}{2}}kn}+}e^{-i\frac{2\pi}{N}k}\cdot\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n+1\right)e^{-i\frac{2\pi}{\frac{N}{2}}k\cdot n}}\]위는 우리가 구한 DFT의 식이다. 이때, 복소 지수 부분을 잘 살펴보면 $e^{-i\frac{2\pi}{N}k}$ 꼴이 반복되고 있음을 확인할 수 있다. 이를 $W^{-k}$로 치환해보자. 그러면
\[f\left(k\right)=\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n\right)W^{-2nk}+}W^{-k}\cdot\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n+1\right)W^{-2nk}}\] \[f\left(k+\frac{N}{2}\right)=\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n\right)W^{-2nk}-}W^{-k}\cdot\sum_{n=0}^{\frac{N}{2}-1}{f\left(2n+1\right)W^{-2nk}}\]가 됨은 자명하다!
그런데 여기서 또 하나의 테크닉 Root of unity를 사용해보자. 이는 N제곱을 했을 때 $1$이 되는 복소수들의 집합을 의미하는데, 이는 복소 평면에 나타냈을 때 원을 $N$등분하여 찍은 점의 좌표와 동치이다. 이를 활용하여 $e^{\frac{2\pi ik}{n}}=W^k=cos\frac{2k\pi}{n}+i\cdot sin\frac{2k\pi}{n}$ 로 쌍을 뽑게 되면 정말 놀랍게도 지수가 상쇄됨을 확인해볼 수 있다!! 이런 아이디어는 어떻게 얻은 것인지 정말 아름답다고 표현할 수 밖에 없을 것 같다.
이로써 다항식을 쌍들을 이용하여 빠르게 곱할 수 있게 되었다. 그러면 이제 새롭게 생성된 쌍들로부터 다항식을 어떻게 다시 복원할 수 있을까? 이는 앞서 설명했던 IDFT, DFT의 역연산을 이용하면 된다.
\[f\left(n\right)=\frac{1}{N}\sum_{k=0}^{N-1}{f\left(k\right)e^{i\frac{2\pi}{N}kn},0\le n<N}\]IDFT 식을 잘 살펴보면 DFT 식에서 복소 지수 부분만 다른 것을 확인할 수 있는데, 복소 지수의 부호가 반대인 경우 오일러 공식을 통해 분해하여 나온 결과가 켤레 복소수가 됨을 확인할 수 있다. 이 덕분에 IDFT는 DFT를 돌리는데, 이때 앞서 설정한 $x_j$의 켤레 복소수를 넣어주면 된다!
이를 통해 큰 수의 곱셈을 각 자릿수를 계수로 가지는 다항식으로 취급하여 계산할 수 있다.
구현은 큰 수 곱셈 2 로 대체한다.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include <bits/stdc++.h>
#define fastio cin.tie(0),cout.tie(0),ios::sync_with_stdio(0)
using namespace std;
const double PI = acos(-1);
typedef complex<double> cpx;
string x,y;
vector<cpx> a,b;
void FFT(vector<cpx> &f, cpx w)
{
int n = f.size();
if(n==1) return;
vector<cpx> e(n/2), o(n/2);
for(int i=0;i<n;i++) (i%2?o:e)[i/2]=f[i];
FFT(e,w*w),FFT(o,w*w);
cpx wp(1,0);
for(int i=0;i<n/2;i++)
{
f[i]=e[i]+wp*o[i];
f[i+n/2]=e[i]-wp*o[i];
wp *= w;
}
}
void mul(vector<cpx> &a, vector<cpx> &b)
{
int n=1,x=a.size()+b.size()-1;
while(n<=a.size()||n<=b.size()) n<<=1;
n<<=1;
a.resize(n),b.resize(n);
vector<cpx> c(n);
cpx w(cos(2*PI/n),sin(2*PI/n));
FFT(a,w),FFT(b,w);
for(int i=0;i<n;i++) c[i]=a[i]*b[i];
FFT(c, cpx(1,0)/w);
int t; vector<int> ans(x+1);
for(int i=x-1;i>=0;i--) t=round((c[i]/cpx(n,0)).real()), ans[i]+=(ans[i+1]+t)/10,ans[i+1]=(ans[i+1]+t)%10;
if(ans[0]) cout << ans[0];
for(int i=1;i<=x;i++) cout << ans[i];
}
int main()
{
fastio;
cin >> x >> y;
for(int i=0;i<x.length();i++) a.push_back(x[i]-'0');
for(int i=0;i<y.length();i++) b.push_back(y[i]-'0');
if(x[0]=='0'||y[0]=='0') cout << 0;
else mul(a,b);
}
To be continued
더 나아가서 FFT의 비재귀 구현, 시간 단축, 실수 오차 줄이기 등의 테크닉에 대해서는 추후에 더 작성할 예정이다.
또한 실수 대신 정수를 사용하는 NTT에 대해서도 글을 작성할 예정이다.