안녕하세요. 지난 포스팅의 디지털 영상 처리 - 선택적 필터링에서는 밴드차단 및 통과 필터, 그리고 노치 필터에 대해서 알아보았습니다. 오늘은 번외로 대부분의 프로그래밍 언어에서 구현되어 있는 고속 푸리에 변환(Fast Fourier Transform; FFT) 알고리즘을 보도록 하겠습니다. MATLAB에는 fft2, 파이썬의 넘파이에서는 np.fft.fft2로 구현이 되어있을 겁니다.
1. 2D DFT 분리성(Separability)
2D DFT는 각 차원에 따라서 수행하는 2개의 1D DFT로 분리될 수 있습니다. 일단, 다시 2D DFT를 상기하면 그 식은 아래와 같습니다.
$$F(\mu, \nu) = \sum_{x = 0}^{M - 1}\sum_{y = 0}^{N - 1} f(x, y)e^{-2j\pi\left(\frac{\mu x}{M} + \frac{\nu y}{N}\right)}$$
그런데 위 식을 저희는 아래와 같이 바꾸어볼 수 있습니다.
$$F(\mu, \nu) = \sum_{x = 0}^{M - 1}e^{-2j\pi\frac{\mu x}{M}}\sum_{y = 0}^{N - 1} f(x, y)e^{-2j\pi\frac{\nu y}{N}}$$
$$ = \sum_{x = 0}^{M - 1} F(x, \nu) e^{-2j\pi\frac{\mu x}{M}}$$
위 식의 변환 원리는 $y$에 대한 합을 수행할 때, $x$를 상수라고 생각하고 하면 $y$축에 대한 단순한 1D DFT임을 알 수 있습니다, 여기서 $F(x, \nu) = \sum_{\mu = 0}^{M - 1} F(\mu, \nu) e^{-2j\pi\frac{\nu y}{N}}$입니다. 따라서 2D DFT는 $f(x, y)$의 각 행에 대해 1D 변환을 계산한 후 그 계산 결과를각 열을 따라서 1D 변환을 계산하여 얻을 수 있습니다.
2. DFT 알고리즘을 이용한 IDFT 계산
이번에는 2D IDFT를 다시 상기해보도록 하겠습니다.
$$f(x, y) = \frac{1}{MN} \sum_{\mu = 0}^{M - 1} \sum_{\nu = 0}^{N - 1} F(\mu, \nu) e^{2j\pi\left(\frac{\mu x}{M} + \frac{\nu y}{N}\right)}$$
위 식에 저희가 복소 공액(complex conjugate)를 취하고 양변에 $MN$을 곱하면 아래의 식을 얻을 수 있습니다.
$$MNf^{*}(x, y) = \sum_{\mu = 0}^{M - 1} \sum_{\nu = 0}^{N - 1} F^{*}(\mu, nu)e^{-2j\pi\left(\frac{\mu x}{M} + \frac{\nu y}{N}\right)}$$
그런데 위 식에서 오른쪽은 $F^{*}(\mu, \nu)$의 2D DFT를 의미합니다. 따라서 저희는 2D DFT를 위한 알고리즘을 이용해서 $F^{*}(\mu, \nu)$를 대입한다면 $MNf^{*}(x, y)$를 얻을 수 있습니다. 그러므로 $f(x, y)$를 얻고 싶다면 위의 2D DFT 결과에 $MN$으로 나눈 뒤 복소 공액을 취하면 2D IDFT를 계산하지 않고도 $f(x, y)$를 얻을 수 있습니다.
3. 고속 푸리에 변환(Fast Fourier Transform;FFT)
저희가 순수한 2D DFT를 계산한다고 가정하면 영상의 크기, 즉 총 $(MN)^{2}$의 곱셈과 덧셈을 수행해야합니다. 만약 $M = N = 1024$라면 총 1조번의 연산량을 수행해야만 하죠. 하지만, 이제부터 소개해드릴 FFT를 이용하면 2천만번의 연산으로 큰 감소를 하게 됩니다. 일단 FFT 알고리즘의 기본 틀은 "연속 배증 방법(successive-doubling method)"에 그 뿌리를 두고 있습니다. 그리고 저희는 2D DFT의 분리성에 의해서 1D DFT에 대한 FFT만 알아보도록 하겠습니다.
원래 1D DFT의 수식은 $F(\mu) = \sum_{x = 0}^{M - 1} f(x)e^{-2j\pi\frac{\mu x}{M}}$이지만 FFT에 대한 이야기를 할 때는 $W^{\mu x}_{M} = e^{-2j\pi\frac{\mu x}{M}}$ 그리고 $M = 2^{n}$이라고 가정한 뒤 아래와 같이 적는 것이 정석입니다. 대부분의 영상의 크기가 256, 512, 1024와 같이 잡는 것을 생각하면 그리 이상한 가정은 아니라고 생각합니다.
$$F(\mu) = \sum_{x = 0}^{M - 1} f(x)W^{\mu x}_{M}$$
그런데, 여기에서 $n$은 양의 정수이기 때문에 양의 정수 $K$에 대해서 $M = 2K$라고도 적을 수 있습니다. 이 식을 위에 대입하면 아래와 같습니다.
$$F(\mu) = \sum_{x = 0}^{2K - 1} f(x)W^{\mu x}_{2K}$$
그리고 위 식을 전개해보도록 하겠습니다.
$$F(\mu) = f(0)e^{-2j\pi\frac{\mu 0}{2K}} + f(1)e^{-2j\pi\frac{\mu 1}{2K}} + f(2)e^{-2j\pi\frac{\mu 2}{2K}} +f(3)e^{-2j\pi\frac{\mu 3}{2K}} + \cdots $$
$$ = \left[f(2 \cdot 0)e^{-2j\pi\frac{\mu (2 \cdot 0)}{2K}} + f(2 \cdot 1)e^{-2j\pi\frac{\mu (2 \cdot 1)}{2K}} + \cdots \right] + \left[f(2 \cdot 0 + 1)e^{-2j\pi\frac{\mu (2 \cdot 0 + 1)}{2K}} + f(2 \cdot 1)e^{-2\pi\frac{\mu (2 \cdot 1 + 1)}{2K}} + \cdots \right]$$
$$ = \sum_{x = 0}^{K - 1}f(2x)W^{\mu (2x)}_{2K} + \sum_{x = 0}^{K - 1} f(2x + 1)W^{\mu (2x + 1)}_{2K}$$
또한 $W^{\mu x}_{K}$의 정의에 의해 저희는 아래의 식들을 얻을 수 있습니다.
$$W^{2\mu x}_{2K} = e^{-2j\pi\frac{2\mu x}{2K}} = e^{-2j\pi\frac{\mu x}{K}} = W^{\mu x}_{K}$$
$$W^{\mu (2x+1)}_{2K} = e^{-2j\pi\frac{\mu (2x + 1)}{2K}} = e^{-2j\pi\frac{\mu x}{K}}e^{-2j\pi\frac{\mu}{2K}} = W^{\mu x}_{K}W^{\mu}_{2K}$$
이 두 식들을 다시 $F(\mu)$에 대입하면 아래와 같은 식이 됩니다.
$$F(\mu) = \sum_{x = 0}^{K - 1} f(2x) W^{\mu x}_{K} + \sum_{x = 0}^{K - 1} f(2x + 1) W^{\mu x}_{K} W^{\mu}_{2K}$$
이번에는 이 식들을 쪼개서 각각 $F_{\text{even}}(\mu)$과 $F_{\text{odd}}(\mu)$를 정의해보도록 하겠습니다.
$$F_{\text{even}}(\mu) = \sum_{x = 0}^{K - 1} f(2x)W^{\mu x}_{K}$$
$$F_{\text{odd}}(\mu) = \sum_{x = 0}^{K - 1} f(2x + 1)W^{\mu x}_{K}$$
그러면 $F(\mu) = F_{\text{even}}(\mu) + F_{\text{odd}}(\mu)W^{\mu}_{2K}$로 줄어들게 됩니다. 그리고 $W^{\mu x}_{M}$과 관련된 2가지 식을 다시 활용합니다. 이 과정에서 오일러 공식을 활용하시면 됩니다.
$$W^{\mu + M}_{M} = e^{-2j\pi\frac{\mu + M}{M}} = e^{-2j\pi\frac{\mu}{M}}e^{-2j\pi} = W^{\mu}_{M}$$
$$W^{\mu + M}_{2M} = e^{-2j\pi\frac{\mu + M}{2M}} = e^{-2j\pi\frac{\mu}{2M}}e^{-j\pi} = -W^{\mu}_{2M}$$
따라서, $F(\mu + K) = F_{\text{even}}(\mu) - F_{\text{odd}}(\mu)W^{\mu}_{2K}$가 됩니다. 이제 분석에 필요한 식들을 하나로 보아보도록 하겠습니다.
$$F_{\text{even}}(\mu) = \sum_{x = 0}^{K - 1} f(2x)W^{\mu x}_{K}$$
$$F_{\text{odd}}(\mu) = \sum_{x = 0}^{K - 1} f(2x + 1)W^{\mu x}_{K}$$
$$F(\mu) = F_{\text{even}}(\mu) + F_{\text{odd}}(\mu)W^{\mu}_{2K}$$
$$F(\mu + K) = F_{\text{even}}(\mu) - F_{\text{odd}}(\mu)W^{\mu}_{2K}$$
이 식들을 이용하여 FFT를 이용하면 큰 계산성능을 향상할 수 있습니다.
