07-1. FFT 유도

 

이 페이지에서는 DFT를 고속 계산하게 해 주는 쿨리-튜키 알고리즘을 소개한다.

쿨리-튜키 알고리즘은 분할정복(Divide and Conquer)기반 알고리즘이다.

 

분할정복은 $O(N^2)$의 계산량을 $O(N\log{N})$으로 만들어주는 대표적인 알고리즘 유형이다. 

N개의 데이터를 $\frac{N}{2}$개씩의 데이터로 나누면서 계산해서 계산량을 줄이면서도, 원래의 결괏값이 나오도록 하는 알고리즘인데, 꽤 이해하기 힘든 알고리즘이다.

 

분할 정복 알고리즘

분할정복 알고리즘이 성립하려면, N개의 데이터를 $\frac{N}{2}$씩 2개로 나눠서 계산한 후 상수 시간 안에 합쳐서, 원래 N개일 때의 계산과 동일해야 한다.  또한, 데이터가 2개일 때의 계산이 간단해야 한다.

 

분할 정복 알고리즘을 적용할 수 있는 조건

• 조건 1: 1/2개의 데이터로 나눠서 계산 후 상수 시간 안에 합쳐서 원래 N개일 때의 계산 값으로 도출 가능
• 조건 2: 2개의 데이터에 대해서 상수 시간 안에 계산 가능

[그림 1] 분할정복 알고리즘 적용가능 조건

N개일 때의 계산량이, $\frac{N}{2}$개에 대해 각각 계산하게 되면 $\frac{N^2}{2}$으로 줄어들게 된다. (계산식이 이렇게 될 수 있는 조건이라야 한다. ) 이렇게 분할해서 바로 계산하는 것이 아니고, 분할한 데이터에 대해서 끝에 2개의 데이터가 될 때까지 계속 분할한다. 

 

1/2로 분할한 데이터를 다시 또 각각 1/2로 분할하는 것을 계속하면, 끝내는 2개의 데이터에 대해서 계산하는 단계까지 내려갈 것이고, 이 2개의 데이터에 대해서 계산한 것을 이제는 위쪽으로 보내서 합치면서 원래 N일 때의 계산을 수행하게 된다.

 이렇게 하게 되면, 분할하는데 $\log{N}$번 수행되고(2개씩 나눴으니깐), 분할한 후 합칠 때 N개에 대한 상수 시간의 계산이 소요되기에 한 레이어당 N개의 계산. 따라서 $N\log{N}$번의 계산으로 결괏값을 계산할 수 있다.

 

DFT에 대한 분할정복 알고리즘 적용

DFT 수식은 위의 분할정복 알고리즘을 적용할 수 있는 2가지 조건을 만족한다. 

 

• 조건 1: 1/2개의 데이터로 나눠서 계산 후 상수 시간 안에 합쳐서 원래 N개일 때의 계산 값으로 도출 가능
• 조건 2: 2개의 데이터에 대해서 상수 시간 안에 계산 가능

위 조건에 맞게끔 DFT의 식을 변형시켜 적용한 것이 쿨리-튜키 알고리즘이다.

 

DFT 식으로부터 분할정복 형태가 되도록 유도해보자.

 

먼저 길이 N인 y(n) 데이터를 2개로 나눠보자. 여기서 N은 짝수로 가정.

하나는 n이 짝수인 것, 그리고 다른 하나는 홀수로 나누면 정확하게 2개로 나눌 수 있을 것이다. (DFT 수식에 사용된 W 함수가 짝수, 홀수에 따라 구분되는 특성이 있기에 짝/홀수로 나누는 것임)

 

 

짝수로 구성된 것을 $p(n)$, 홀수로 구성된 것을 $q(n)$이라 하자.

 

$$ p(n)=y(2n), \; q(n)=y(2n+1), \; 0 \le n \le {\frac{N}{2}-1} \tag{1}$$

 

이제 N개의 원소를 가진 데이터 y(n)의 DFT식을 p(n)과 q(n)을 이용해서 표현하면,

 

$$ \begin{align} \\ Y(k) &= \sum _{n=0}^{N-1}y(n)e^{-i\frac{2\pi}{N}kn} \\ &= \sum_{n=even}^{N-1}{y(n)e^{-i\frac{2\pi}{N}kn}}+\sum_{n=odd}^{N-1}{y(n)e^{-i\frac{2\pi}{N}kn}} \\ &=\sum_{n=0}^{\frac{N}{2}-1}{y(2n)e^{-i\frac{2\pi}{N}k(2n)}} + \sum_{n=0}^{\frac{N}{2}-1}{y(2n+1)e^{-i\frac{2\pi}{N}k(2n+1)}}  \\ &=\sum_{n=0}^{\frac{N}{2}-1}{p(n)e^{-i\frac{2\pi}{\frac{N}{2}}kn}} + \sum_{n=0}^{\frac{N}{2}-1}{q(n)e^{-i\frac{2\pi}{N}k(2n)}\cdot e^{-i\frac{2\pi}{N}k}} \\ &=\sum_{n=0}^{\frac{N}{2}-1}{p(n)e^{-i\frac{2\pi}{\frac{N}{2}}kn}} +e^{-i\frac{2\pi}{N}k} \sum_{n=0}^{\frac{N}{2}-1}{q(n)e^{-i\frac{2\pi}{\frac{N}{2}}kn} }, \;  0\le k \le{n-1} \tag{2}\end{align} $$

 

식(2)에서 유도된 제일 마지막 식을 보면,  $\frac{N}{2}$개의 짝수항에 대한 DFT와, 홀수항에 대한 DFT식을 더하는 형태이다. 더 정확하게는, 짝수항 DFT와 어떤 상수값($e^{-i\frac{2\pi}{N}k}$)이 곱해진 홀수항 DFT의 합으로, 짝수항에 대한 DFT를 $P(k)$, 홀수항에 대한 DFT를 $Q(k)$라 하면 (식 2)는 다음과 같이 쓸 수 있다.

 

$$ Y(k)  = P(k) + e^{-i\frac{2\pi}{N}k}Q(k), 0\le k \le {N-1} \tag{3}$$

 

즉, N개 데이터에 대한 DFT식이, $\frac{N}{2}$개의 두 DFT 식의 합으로 표현된 것으로, 분할정복의 기본 조건 형태를 만족하게끔 변형이 된 것이다.

 

식의 유도과정을 좀 더 자세히 보자.

그런데, (식 2) 및 (식 3)이 왜 계산량을 적게 하는지, 그리고 왜 이렇게 나눈 것이 의미가 있는 걸까?

 

(식 3)에서 $W^{k}=e^{-i\frac{2\pi}{N}k}$로 놓고, k=0에서부터 하나씩 써보자.

 

$$ \begin{align}&Y(0)=P(0)+W^0Q(0) \\ &Y(1)=P(1)+W^1Q(1) \\ &... \\ &Y(N-1)=P(N-1)+W^{(N-1)}Q(N-1) \tag{4}\end{align}$$

 

(식 4)에서는 총 N번의 곱셈 연산과 N번의 덧셈 연산이 일어난다. Y(0) 계산을 위해 1번씩, Y(1)을 위해 1번... 이렇게 Y(7)까지 구하려면 N번의 곱셈과 N번의 덧셈을 하면 된다.

알고리즘의 연산량을 계산할 때 N+N은 그냥 $O(N)$이기에 그냥 N번 연산이 일어난다고 계산하자. (P와 Q에 대해서는, 이 단계에서는 이미 계산되어 있다고 가정하는 것이다.)

 

즉, N개에 원소에 대해서 2로 나눴을 때 연산량이 N회 발생.

 

이제 P와 Q를 구하기 위해서, P와 Q를 각각 절반으로 나눠서 계산하자. P를 절반으로 나는 것을 A, B라 하고, Q를 절반으로 나는 것을 C, D라고 하면,

 

$$ \begin{align}&P(0)=A(0)+W^0B(0) \\ &P(1)=A(1)+W^1B(1) \\ &... \\ &P(\frac{N}{2}-1)=A(\frac{N}{2}-1)+W^{(\frac{N}{2}-1)}Q(\frac{N}{2}-1) \end{align}$$

$$ \begin{align}&Q(0)=C(0)+W^0D(0) \\ &Q(1)=C(1)+W^1D(1) \\ &... \\ &Q(\frac{N}{2}-1)=C(\frac{N}{2}-1)+W^{(\frac{N}{2}-1)}D(\frac{N}{2}-1) \end{align}$$

 

P와 Q를 구하기 위해서도 역시 N번의 연산이 필요.

 

이러한 과정을 원소가 2개가 될 때까지 계속하게 되고, 원소가 2개씩 남았을 때도 역시 N번의 연산 필요.

 

N에서부터 시작해서 2개씩 나눠서 최종 원소가 2개 남을 때까지 수행되는 횟수는 $\log_2{N}$

 

따라서, 전체 연산 횟수는 $N\log{N}$이 된다. (밑 2는 생략해서 표시)

 

---

(식  3)에서, DFT값의 주기성과 복소지수함수의 성질을 이용하면, 좀 더 단순화시켜서 $\frac{N}{2} \log{N}$의 횟수에 계산이 되도록 할 수 있다.

 

(식 3)에서 P(k)와 Q(k)는 각각 $\frac{N}{2}$개의 원소에 대한 DFT이다. 따라서, $\frac{N}{2}$ 주기에 따라 같은 값을 가진다. (DFT라는 것이 N을 한 주기로 해서 계산되는 것이기에, Y(0) = Y(0+N)과 같이 N을 주기로 해서 같은 값을 가진다.)

따라서, $P(k) = P(k+\frac{N}{2})$이고, $Q(k)=Q(k+\frac{N}{2})$이다.

 

또한 복소지수함수의 성질에 의해서 $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}$가 된다.

($e^{-i\pi}=-1$이다. 계산은 $e^{-i\pi}=\cos{\pi}-i\sin{\pi}=-1$)

 

이제 (식 3)을 다음과 같이 쓸 수 있다.

$$ \begin{align}Y(k) = P(k) + e^{-i\frac{2\pi}{N}k}Q(k), \; &0 \le k \le {\frac{N}{2}-1} \\ Y(k+\frac{N}{2})=P(k)-e^{-i\frac{2\pi}{N}k}Q(k), \; &0 \le k \le {\frac{N}{2}-1} \tag{4} \end{align} $$

 

(식 4)에서 $e^{-i\frac{2\pi}{N}k}=W^k$로 놓으면,

$$ \begin{align}Y(k) = P(k) + W^kQ(k), \; &0 \le k \le {\frac{N}{2}-1} \\ Y(k+\frac{N}{2})=P(k)-W^kQ(k), \; &0 \le k \le {\frac{N}{2}-1} \tag{5} \end{align} $$

 

---

(식 5)에 대해서 N=8의 경우를 예로 들어보자.

 

$$ \begin{align} &Y(0)=P(0)+W^0Q(0) \\ &Y(1)=P(1)+W^1Q(1) \\ &Y(2)=P(2)+W^2Q(2) \\ &Y(3)=P(3)+W^3Q(3) \\ \\ &Y(4)=P(0)-W^0Q(0) \\ &Y(5)=P(1)+W^1Q(1) \\ &Y(6)=P(2)+W^2Q(2) \\ &Y(7)=P(3)+W^3Q(3) \end{align}$$

 

Y(0)~Y(3)까지는 (식 5)의 윗 식을, Y(4)~Y(7)까지는 (식 5)의 아랫 식을 이용한 것이다.

즉, Y(0)~Y(7)까지의 값을 4개의 P, W, Q값을 가지고 계산할 수 있는 것이다. 

 

위 식을 보면 대칭적인 규칙을 가지고 있다. 

그래프로 그려보면 그 규칙이 잘 보인다.

$ Y(0)=P(0)+W^0Q(0) $를 그래프로 표현한 것이다. 

$Y(1)=P(1)+W^1Q(1)$ 까지 넣어서 그려보자.

$Y(2)$와 $Y(3)$까지 그리면,

 

$Y(4)=P(0)-W^0Q(0)$로, 다시 $Y(0)$를 구할 때 사용했던 $P(0), W^0, Q(0)$를 이용하면서, 덧셈이 아니라 뺄셈을 수행한다.

 이제 남은 Y(5)~Y(7)까지를 다 그려 넣으면,

 

---

$P(0)$~$P(3)$와 $Q(0)$~$Q(3)$도 동일한 과정으로 구할 수 있다. 

 

$P(0)$~$P(3)$를 구하는데 사용되는 데이터를 $p(0)$~$p(3)$라고하고, $p(0)$~$p(3)$에서 짝수항을 $u$, 홀수항을 $v$로 표현하면, 

$$ \begin{align}&u(0)=p(0) \\ &u(1)=p(2) \\ &v(0)=p(1) \\ &v(1)=p(3) \end{align}$$

 

이제 $P(0)$~$P(3)$을 (식 5)에서 N=4를 넣고, $u$와 $v$에 대한 DFT값인 $U$와 $V$를 이용하여 표현하면, 

 

$$ \begin{align}&P(0)=U(0)+ W^0V(0) \\ &P(1)=U(1)+ W^1V(1) \\ &P(2)=U(0)- W^0V(0) \\ &P(3)=U(1)- W^1V(1) \end{align}$$ 

 

$P(0)$~$P(3)$를 구하는 것을 그래프로 그리면 다음과 같다.

 

마찬가지로 $Q(0)$~$Q(3)$에 대해도 $P(0)$~$P(3)$를 구하는 것과 유사하게 구할 수 있다.

$Q(0)$~$Q(3)$를 구하는데 사용될 데이터를 $q(0)$~$q(3)$라고 하고, 짝수항을 $r$, 홀수항을 $s$로 표현하면,

 

$$ \begin{align}&Q(0)=R(0)+ W^0S(0) \\ &Q(1)=R(1)+ W^1S(1) \\ &Q(2)=R(0)- W^0S(0) \\&Q(3)=R(1)- W^1S(1) \end{align}$$

 

---

이제 $U, V, R, S$를 구하는 것은, 각각 y(n)의 2개의 데이터를 이용해서 구하면 된다.

 

2개 데이터를 이용해서 DFT 값을 구하는 것은, 실제 DFT 식에 넣어보면 굉장히 단순한 실수연산으로 변환이 된다.

 

길이가 2인 x(n)에 대한 DFT 값 X(k)를 구하는 것을 생각해보자. DFT식은,

 

$$ \begin{align}X(k)&=\sum _{n=0} ^{1}x(n)e^{-i\frac{2\pi}{2}kn} \\ &=x(0)e^{-i\frac{2\pi}{2}k\cdot 0} + x(1)e^{-i\frac{2\pi}{2}k\cdot 1} &= x(0) + x(1)e^{-i\pi k}\end{align} $$

 

여기서 $k$의 값 {0,1}을 차례로 대입해서 계산해보면,

 

$ k=0$의 경우, $X(0) = x(0) + x(1)e^{-i\pi 0} = x(0) + x(1)$  --> $e^0=1$이기에,

$ k=1$의 경우, $X(1) = x(0) + x(1)e^{-i\pi 1} = x(0) -x(1)$  --> $e^{-i\pi}=-1$이기에

 

정리하면, 

 

$$ \begin{align} &X(0)=x(0)+x(1) \\ &X(1)=x(0)-x(1)\end{align}$$

 

2개 포인트에 대한 DFT 계산은, 이처럼 두 개 값의 덧셈과 뺄셈으로 간단히 계산될 수 있다. 

 

---

이처럼, 8개 데이터에 대한 FFT의 최종 값 Y(k)의 계산은,

• 3단계: Y(k) 값은, 4 포인트에 대한 FFT값인 P(k)와 Q(k)에 의해서 계산되고,
• 2단계: P(k)와 Q(k)는, 2 포인트에 대한 FFT 값인 U(k), V(k), R(k), S(k)에 의해 계산되고,
• 1단계: U(k), V(k), R(k), S(k)는, 원래 주어진 y(n)에 의해 계산된다.

실제 계산을 할 때 제일 아랫단에서의 y(n)에 의해 U(k), V(k), R(k), S(k)를 구하는 것이 제일 먼저 진행되기 때문에, 이 과정을 1단계로 표현한 것임

---

(식 5)가 FFT의 수식이고 알고리즘의 전부다.

 

식만 봐서는 어떻게 계산되는지 감이 잘 안 잡히고, 위에서 8개 포인트의 데이터에 대한 수식을 예로 들었으나, 아직 완벽하게 이해되지 않을 수 있다. 

 

실제 8개 포인트 데이터에 대한 FFT 계산을 손으로 한 단계씩 계산해보면, 어떻게 계산이 되는지 체득할 수 있을 것이다.  다음 페이지에서 그렇게 계산해볼 것이다.

 

-끝-

이전글: 07. FFT (Fast Fourier Transform, 고속 푸리에 변환)

다음글: 07-2. FFT 예제를 손으로 풀어보며 이해하기

다음다음글: 07-3. 엑셀에서 FFT 계산하기 (엑셀 자체 FFT 기능)