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

[경고] 아래 글을 읽지 않고 "이산 푸리에 변환"을 보면 바보로 느껴질 수 있습니다.

1. 푸리에 급수의 시작

2. 형보다 나은 아우: 푸리에 변환

3. 푸리에 사인 및 코사인 변환

4. 차분 방정식

이산 푸리에 변환(discrete Fourier transform, DFT)은 푸리에 변환(Fourier transform)을 컴퓨터 상에서 구현할 때 사용하는 중요 개념이다. 이산(離散, discrete) 푸리에 변환은 명칭에도 있듯이 연속적이지 않고 이산적으로 변하는 함수에 대한 푸리에 변환이다. 컴퓨터는 원론적으로 연속 함수(continuous function)를 다룰 수 없기 때문에, 수치 해석에서는 컴퓨터에 저장하기 쉬운 이산화(離散化, discretization)를 적극적으로 채택한다. 

[그림 1] 이산 푸리에 변환의 적용 예시(출처: wikipedia.org)

연속적으로 변하는 함수 $f(t)$를 이산적으로 바꾸려면, 연속 함수에 [그림 2]와 같은 표본화(sampling)를 적용한다. 표본화하는 주기를 $\Delta t$, 표본 개수를 $M$이라 하면, $f(t)$를 이산화한 전체 표본화 시간은 $T$ = $M \Delta t$가 된다. 하지만 표본화 시간인 $0 \le t < M\Delta t$에서는 함수 $f(t)$의 값을 알지만, 나머지 시간에서는 함수값을 알 수 없게 된다. 이때는 $f(t)$를 과감하게 주기 $T$를 가진 주기 함수(periodic function)라 가정한다. 그래서 표본화된 함수 $f_\text{sa}(t)$는 주기가 $T$이면서 $t$ = $m \Delta t$[$m$ = $0, \pm 1, \pm 2, \cdots$]에서만 값이 이산적으로 존재하는 주기 함수가 된다.

[그림 2] 연속 함수의 표본화 예시(출처: wikipedia.org)

DFT 관계식을 유도하기 위해, 먼저 복소 푸리에 급수(complex Fourier series)와 푸리에 변환(Fourier transform) 공식을 고려한다.

                  (1)

                  (2)

여기서 $\omega_0$ = $2 \pi/T$이다. 푸리에 변환은 $F(\omega)$ = $\lim_{T_\text{tot} \to \infty} \widetilde{F}_m T_\text{tot}$로 정의한다. 여기서 $T_\text{tot}$는 $f(t)$가 존재하는 전체 시간이다. 함수 $f(t)$가 전체 표본화 시간 $T$를 주기로 반복되고 이상적인 표본화는 디랙 델타 함수(Dirac delta function)라고 가정하면, 푸리에 변환보다는 식 (1)에 있는 푸리에 급수가 표본화된 함수 $f_\text{sa}(t)$의 해석에 더 적절하다. 그러면 $f_\text{sa}(t)$에 대해, 푸리에 변환에서 유추한 새로운 DFT의 푸리에 계수 $F_k$를 다음처럼 얻을 수 있다.

                  (3)

여기서 $\delta(t)$는 디랙 델타 함수, 표본화 계수 $f_m$ = $f(m \Delta t)$이다. 식 (3)에 $k$ = $k + M$을 대입하면, $F_k$ = $F_{k+M}$이 항상 성립한다. 이 결과는 $F_k$가 주기 $M$을 가지고 반복되는 계수임을 뜻한다. 푸리에 계수 $F_k$를 식 (1)에 대입해서 $F_k$의 무한 합으로 $f_\text{sa}(t)$를 다음과 같이 표현한다.

                  (4a: 푸리에 급수 이용한 정의)

                  (4b)

여기서 $T$ = $M \Delta t$, $f_\text{sa}(t)$의 정의역을 $0 \le t < T$로 제한, 무한한 $k$[ = $0, \pm1, \pm2, \cdots$]를 유한한 $0 \le k + mM < M$[$m$ = $0, \pm1, \pm2, \cdots$]로 바꾸며, $\delta(2 \pi t / \Delta t)$ = $\Delta t / (2 \pi) \delta(t)$이다. 식 (3)의 첫째식과 식 (4)의 마지막식을 비교해서 다음 결과를 유도한다.

                  (5)

따라서 이산 푸리에 변환쌍(DFT pair)을 다음처럼 표현할 수 있다.

                  (6)

이산 푸리에 역변환에 $1/M$이 나오는 이유는 적분의 이산화로 쉽게 설명할 수 있다. 미분소 $d \omega$의 이산화는 $\Delta \omega$ = $\omega_0$ = $2 \pi / T$ = $2 \pi / (M \Delta t)$이다. 이 결과를 식 (2)에 대입하고 $\Delta t$ = $1$로 두면, 최종적으로 $1/M$만 남는다. 식 (6)을 한 차원 더 확장해서 2차원 이산 푸리에 변환쌍을 정의한다.

                  (7)

식 (6)의 첫째식을 둘째식에 대입하거나 둘째식을 첫째식에 대입하면, 다음 항등식이 유도된다.

                  (8)

여기서 $\delta_{ml}$은 크로네커 델타(Kronecker delta)이다. 식 (8)의 최종식은 등비 급수(geometric series)를 적용해서 간단히 얻을 수도 있다.

                  (9)

이산 푸리에 변환을 기본으로 해서 더욱 다양한 이산 변환을 생성할 수 있다.

   1. 이산 코사인 변환(discrete cosine transform, DCT)   

[그림 1.1] 여러 가지 이산 코사인 변환(출처: wikipedia.org)

이산 코사인 변환(discrete cosine transform, DCT)은 DFT의 현실적인 변형이다. DFT는 푸리에 변환을 실제 수치 해석에 적용할 때 매우 유용하다. 하지만 DFT의 계수 $F_k$는 복소수라서 다루기 불편하다. 만약 DFT 대신 DCT를 쓰면, 변환 결과는 실수가 되기 때문에 아주 편해진다. 현실적인 신호 처리에 필수적인 DCT는 아메드Nasir Ahmed(1940–)가 1972년아메드 32세, 박정희 정부 시절에 제안했다. [그림 1.1]에 있는 여러 DCT 중에서 가장 쉬운 형태인 DCT-I을 보자. [그림 1.1]처럼 모든 DCT는 우함수(even function) 특성을 가지기 때문에, DFT 기준으로 보면 $f_m$ = $f_{-m}$이 항상 성립한다. 또한 표본화된 함수는 계속 반복되므로, DCT는 주기적인 특성이 나타난다. 예를 들어, DCT-I의 주기는 $2M-2$로 나온다.[그림 1.1에서 빨간색 표본이 $M$개, 이어지는 검정색 표본은 $M-2$개가 된다.] 하지만 DCT를 적용하는 함수가 반드시 우함수이어야 할까? 아니다. 함수 $f(t)$는 우함수일 필요는 없고, 어떠한 $f(t)$이든지 [그림 1.1]과 같이 대칭적으로 함수값을 배정하여 우함수를 강제적으로 만들 수 있다. 즉, [그림 1.1]에서 빨간색은 직접 측정해서 대상의 정보를 가진 실제 표본이고, 검정색은 각 DCT 방법에 맞도록 실제 표본으로 생성한 부차적인 계산 표본이다. [그림 1.1]의 DCT-I처럼 $0 \le m < M$에서 표본화한 함수값[그림 1.1에서 빨간색]을 $x_m$ = $f(m \Delta t)$라 한다. 그러면 표본화한 범위를 넘는 $M \le m < 2M-2$에 대해, 우함수 특성인 $f_m$을 다음처럼 표현할 수 있다.

                  (1.1)

여기서 $f_m$은 $m$ = $M-1$ 기준으로도 우함수이다. 의도적으로 주기 $2M-2$를 가진 우함수로 만든 DFT의 표본화 계수 $f_m$을 식 (6)에 대입해 본다.

                  (1.2)

여기서 $F_k$의 주기도 $2M-2$이다. DCT-I의 표본화 및 푸리에 계수를 각각 $x_m$ = $f_m$, $X_k = F_k /2$로 정의해 정리하면, DCT-I의 변환식을 얻을 수 있다.

                  (1.3a)

                  (1.3b)

                  (1.3c)

여기서 $\alpha_{m, M}$ = $1+\delta_{m0} + \delta_{m,M-1}$이다. DCT-I의 표본화 계수 $x_m$이 실수라서 $X_k$도 당연히 실수가 된다. 또한 식 (1.3)에 의해 $X_k$는 우함수이면서 $k$ = $M-1$ 기준으로도 우함수가 된다. 즉, $X_k$ = $X_{-k}$ = $X_{2M-2-k}$[$k$ = $0, 1, \cdots, M-1$]가 성립한다. 그러면 $k$의 범위를 $0 \le k < M$으로 제한해도 문제는 없다. 이 특성을 식 (6)에 대입해서 $x_m$을 위한 DCT-I의 역변환식도 얻는다.

                  (1.4a)

                  (1.4b)

                  (1.4c)

여기서 $0 \le m < M$이다. 이산 코사인 역변환에 등장한 $2/(M-1)$은 푸리에 코사인 역변환(inverse Fourier cosine transform, IFCT)으로 쉽게 설명한다. IFCT에서 $\pi/2 \cdot d \omega$는 $2/\pi \cdot 2 \pi /[(2M-2)\Delta t]$ = $2/(M-1)$으로 이산화되어서 식 (1.4b)와 같은 상수가 나온다. 식 (1.3)와 (1.4)에서 새로운 DCT-I 공식을 얻었지만, DCT의 표현식이 간결해 보이지는 않는다. 그래서 DFT의 표본화 계수 $f_m$을 다음처럼 다시 정의해서 DCT-I과는 다른 공식 전개를 해본다.

                  (1.5)

여기서 $f_m$의 주기는 $4M$이다. 식 (1.2)와 비슷하게 식 (1.5)를 식 (6)에 대입한다.

                  (1.6)

여기서 $F_k$의 주기도 $4M$이다. DFT의 푸리에 계수를 바꾸어 $X_k = F_k /2$로 정의하면, 식 (1.6)에 의해 [그림 1.1]에 있는 DCT-II의 푸리에 계수를 얻는다.

                  (1.7)

여기서 $0 \le k < M$이다. 식 (1.3)과 비교하면 식 (1.7)은 매우 깔끔한 DCT-II의 표현식이다. 그래서 보통 DCT라고 하면 DCT-II를 의미한다. DCT-II 푸리에 계수 $X_k$의 주기 특성을 알기 위해 $X_{2M \pm k}$를 식 (1.7)에 대입하면, [그림 1.1]에 있는 DCT-III와 비슷하게 $X_{2M \pm k}$ = $-X_k$를 얻을 수 있다. 푸리에 계수 $X_k$의 주기는 분명 $4M$이지만, $X_k$와 $X_{2M \pm k}$는 서로 기함수 관계이기도 하다. 그러면 식 (1.7)에서 $k$의 범위를 $0 \le k < M$로 제한해도 문제없다. 이를 이용해 DCT-II의 역변환식을 구해보자.

                  (1.8)

여기서 $X_{2M}$ = $-X_0$. 기함수 관계 $X_{2M \pm k}$ = $-X_k$를 다시 식 (1.8)에 적용해서 간결하게 정리한다.

                  (1.9)

여기서 $0 \le m < M$, $X_M$ = $0$이다. 노이만 수(Neumann number) $\varepsilon_k$를 이용해 식 (1.9)을 더 간단하게 표현한다.

                  (1.10)

여기서 $\varepsilon_m$ = $2 - \delta_{m0}$, $\delta_{ml}$은 크로네커 델타(Kronecker delta)이다. 식 (1.3)와 (1.10)에서 유추해서 [그림 1.1]에 있는 DCT-III의 변환식도 정의할 수 있다.

                  (1.11)

여기서 $0 \le k < M$. DCT-II와 DCT-III은 비율 상수가 약간 다르지만 서로 역변환 관계를 가진다. 식 (7)에 있는 2차원 DFT와 비슷하게 2차원 이산 코사인 변환쌍을 정의할 수 있다. 예를 들어, DCT-II의 변환쌍은 다음과 같다.

                  (1.12)

                  (1.13)

[그림 1.1]에 4가지 종류의 DCT를 소개했지만, DCT의 표본화 계수 $x_m$에 기반을 둔 DFT의 표본화 계수 $f_m$을 다양하게 정의해서 다채로운 우리만의 DCT를 정의할 수도 있다. DCT를 어떻게 정의하든지 우함수 혹은 기함수의 대칭을 사용하므로, 표본화된 함수의 주기는 원래 경우보다 배 이상 늘어난다. 주기 $T$가 늘어나면, $\omega_0$ = $2 \pi/T$의 관계에 의해 주파수 영역의 계수 간격이 줄어든다. 따라서 DFT보다는 DCT의 푸리에 계수가 저주파에 더 많이 모이므로, 고주파 성분을 자르는 손실 있는 영상 압축 등에 DCT를 많이 사용한다.

DCT-I으로 얻어지는 풍성한 결과 중의 하나는 다음과 같은 항등식이다. 식 (1.4b)에 식 (1.3b)를 넣고 간단히 정리하면 곧 바로 증명이 된다.

                  (1.14a)

                  (1.14b)

여기서 $m$과 $m'$가 변하는 범위는 $0, 1, \cdots, M-1$이다. DCT-I을 쓰는 대신 식 (9)에 실수부를 선택하고 첨자 $k$를 두 영역으로 정리해도 유도가 가능하다.

                  (1.15a)

                  (1.15b)

                  (1.16)

여기서 $k$의 영역은 $k_1$ = $0, 1, \cdots, M-2$와 $k_2$ = $M-1, M, \cdots, 2M-3$으로 분해된다. 식 (1.15a)에서 $m$이 짝수라면, $k$와 $k + M-1$ 항은 서로 같은 값이어서[$\because$ $\pi m / (M-1) \cdot (M-1)$ = $m \pi$] $k_1$과 $k_2$를 가진 급수로 나누어진다. 다만 $m$이 홀수인 경우는 약간 다르게 접근한다. 새로운 첨자 $k_1$과 $k_2$로 만든 급수 합은 부호가 반대라서[$\because$ 이 두 급수를 합한 값은 자동적으로 0], $k_1$ 전체[= $0, 1, \cdots, M-2$]와 $M-1$을 가진 급수를 다시 $k_1$과 $M - 1 - k_1$인 항으로 정리해서 식 (1.15b)를 증명한다. 이때 $m$ = $0$인 급수는 합이 $M$이 되어 식 (1.15b)와 다르므로, $k$ = $0$과 $M-1$의 값을 반으로 줄여서 식 (1.15b)가 되게 한다.

   2. 이산 사인 변환(discrete sine transform, DST)   

[그림 2.1] 여러 가지 이산 사인 변환(출처: wikipedia.org)

이산 코사인 변환과 동일한 방법으로 푸리에 사인 변환(Fourier sine transform)의 이산화인 이산 사인 변환(discrete sine transform DST)을 정의할 수 있다. 표본화 계수가 $m$ = $0$과 $M-1$을 기준으로 기함수라 가정해서 DFT를 적용한다.

                  (2.1)

                  (2.2)

식 (2.1)의 둘째식에 강제로 $m$ = $M-1$을 넣으면, $f_{M-1}$ = $-f_{M-1}$이 되어서 $f_{M-1}$ = $0$이 되어야 한다. 비슷하게 $m$ = $2M-2$인 경우에 $f_{2M-2}$ = $f_0$ = $-f_{0}$으로 인해 $f_0$ = $0$으로 나온다. 그래서 식 (1.3)처럼 $x_m$ = $f_m$, $X_k$ = $i F_k/2$로 놓고 식 (2.2)를 깔끔히 정리한다.

                  (2.3)

여기서 $x_0$ = $x_{M-1}$ = $0$이다. 식 (1.3)과 같은 형태로 식 (2.3)을 다시 기술한 공식은 DST-I이라 부른다.

                  (2.4)

표본화 계수 $x_m$과 유사하게 식 (2.4)에 의해 $X_0$ = $X_{M-1}$ = $0$이다. 또한 $X_k$는 $k$ = $0$과 $M-1$에 대해 기함수를 이룬다. 즉, $x_m$처럼 $X_k$ = $-X_{2M-2-k}$가 성립한다. 식 (2.4)에 대응하는 이산 사인 역변환(inverse discrete sine transform, IDST)은 다음과 같다.

                  (2.5a)

                  (2.5b)

이산 사인 변환인 식 (2.4)에서 $x_0$과 $x_{M-1}$은 0이라서 의미가 없으므로, 이 두 값을 제외해서 전체 표본수가 $M$이 되도록 식 (2.4)와 (2.5)의 $M$을 $M+1$로 늘린다.

                  (2.6)

                  (2.7)

여기서 $m$ = $0, 1, \cdots, M-1$, $k$ = $0, 1, \cdots, M-1$이다. 식 (2.6)을 DST-I이라 이름 붙이기도 한다.

식 (2.4)를 식 (2.5b)에 대입해서 DCT와 비슷한 항등식을 만들 수 있다.

                  (2.8)

여기서 $m$과 $m'$의 범위는 $1, 2, \cdots, M-2$로 한정한다. 혹은 식 (2.8)의 좌변에 삼각 함수 항등식을 적용하고 식 (1.15b)를 넣어서 더 쉽게 증명하기도 한다.

                  (2.9)

여기서 계산을 쉽게 하기 위해 급수 범위를 $0 \le k \le M-1$로 늘린다.

   3. 이산 혼합 푸리에 변환(discrete mixed Fourier transform, DMFT)   

이산 혼합 푸리에 변환(discrete mixed Fourier transform, DMFT)은 이산 사인 및 코사인 변환을 섞어서 정의한다[1]. DMFT를 이산화하기 전의 적분 변환은 바로 혼합 푸리에 변환(mixed Fourier transform, MFT)이다.

                  (3.1)

여기서 $x_m$ = $m \Delta x$, $p_k$ = $k \pi / [(M-1)\Delta x]$ = $k \Delta p$, $\Delta p\Delta x$ = $\pi /(M-1)$, $\alpha$는 $x_0$ 위치의 경계 조건을 결정하는 복소 상수이다. 조건 $\Re[\alpha] \le 0$에 대한 DMFT의 역변환은 다음처럼 공식화된다.

             (3.2)

혼합 푸리에 역변환(inverse mixed Fourier transform, IMFT)의 특성으로 인해, $\Re[\alpha] > 0$인 역변환은 다소 복잡해진다.

             (3.3)

만약 $\alpha$ = $0$이면, DMFT는 전형적인 이산 코사인 변환으로 간략화된다. 하지만 혼합 푸리에 변환을 만드는 적분을 단순히 이산화해서 표현하기 때문에, 아래 적분과 다르게 $e^{-\alpha x_m}$과 $\alpha \sin(p_k x_m) - p_k \cos (p_k x_m)$의 급수에는 직교성이 성립하지 않는다.

                  (3.4)

                  (3.5)

결국 식 (3.2)와 (3.3)에 변형을 가해서 이산화한 두 함수 사이에 직교성이 항상 생기도록 한다.

                  (3.6)

여기서 $\operatorname{Sa}(x)$ = $\sin x \mathbin{/} x$, $M$이 충분히 클 때에 $\operatorname{Sa}(p_k \Delta x)$ $\approx$ $1$, $\operatorname{Sa}(x)$를 도입함으로써 $p_k$와 $\cos(p_k x_m)$의 곱이 $\sin (p_k x_m)$의 차분(差分, difference)으로 바뀐다. 식 (3.6)에 등장한 차분은 중심 $x_m$에 대해 대칭적으로 차이를 만들어서 대칭 차분(symmetric difference) $\Delta_s$가 된다. 여기에 곱의 차분 공식을 적용해서 차분 관계를 $f_m$으로 바꾼다.

             (3.7)

식 (3.7)은 근사식이지만 $M$이 충분히 크면 항등식에 가까워진다. 근사이지만 타당한 식 (3.7)의 최종 결과를 식 (3.6)에 대입해서 $f_m$에 대한 2계 차분 방정식(the second order difference equation)이 생기도록 한다.

                  (3.8)

만약 $M$을 충분히 키우면, 근사식 (3.8)을 등식처럼 생각할 수 있다. 또한 식 (3.8)은 정확히 혼합 푸리에 변환의 이산화가 된다. 또 다른 관점으로 식 (3.8)은 이산 사인 변환이라서, 식 (2.5b)에 따라 $g_m \mathbin{/} (2 \Delta x)$의 역변환이 바로 얻어진다.

                  (3.9)

여기서 $g_0$ = $g_{M-1}$ = $0$, $m$ = $1,2,\cdots,M-2$이다. 미분 방정식의 특수해와 일반해 개념을 써서, 식 (3.9)를 $f_m$에 대해 풀어서 나오는 해는 $f_m$ = $f_{m,p} + f_{m,g}$로 둔다. 여기서 $f_{m,p}$와 $f_{m,g}$는 각각 차분 방정식의 특수해와 일반해이다. 식 (3.9)를 만족시키는 특수해 $f_{m,p}$는 식 (3.6)을 참고해서 설정한다.

                  (3.10)

여기서 $A_k$는 결정해야 하는 상수이다. 식 (3.10)을 식 (3.9) 넣고 좌변과 우변을 비교해서 $A_k$를 구한다.

                  (3.11)

따라서 $A_k$는 식 (3.11)의 마지막식에 나오므로 $f_{m,p}$를 완전히 공식화한다.

                  (3.12)

식 (3.9)에 있는 2계 차분 방정식의 일반해 $f_{m,g}$는 특성 방정식을 풀어서 정한다.

                  (3.13)

여기서 $\Delta_s (r^m) \mathbin{/} (2\Delta x)$ = $-\alpha r^m$, $c_1, c_2$는 경계 조건으로 계산하는 상수이다. 식 (3.12)와 (3.13)을 합쳐서 이산 혼합 푸리에 역변환(inverse discrete mixed Fourier transform, IDMFT)을 확정한다.

                  (3.14)

여기서 일반해를 표현할 때에 $(-r_1)^{-m}$ 대신 대칭적인 $(-r_1)^{M-m-1}$을 쓰기도 한다.

[참고문헌]

[1] J. R. Kuttler and R. Janaswamy, "Improved Fourier transform methods for solving the parabolic wave equation," Radio Sci., vol. 37, no. 2, pp. 5-1–5-11, Apr. 2002.

[다음 읽을거리]

1. 삼각 함수 보간

2. Z 변환

3. 고속 푸리에 변환