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

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

1. 푸리에 급수의 시작

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

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

4. 차분 방정식

5. 퇴플리츠 행렬

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

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

연속적으로 변하는 함수 $f(t)$를 이산적으로 바꾸려면, 연속 함수에 [그림 2]와 같은 표본화(sampling)를 적용한다. 표본화 주기 혹은 간격(sampling period or interval)을 $\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$로 제한, $\delta(2 \pi t / \Delta t)$ = $\Delta t / (2 \pi) \delta(t)$이다; 식 (4b)에서 첫째 줄의 무한한 $k$[= $0, \pm1, \pm2, \cdots$]를 나눗셈의 원리를 이용해 둘째 줄에서 유한한 $0 \le k + mM < M$[$m$ = $0, \pm1, \pm2, \cdots$]으로 바꾸며, 둘째 줄의 $m, k$는 각각 몫(quotient)과 나머지(remainder)가 된다. 식 (3)의 첫째식과 식 (4)의 마지막식을 비교해서 다음 결과를 유도한다.

                  (5)

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

                  (6a)

이산 푸리에 역변환(inverse discrete Fourier transform, IDFT)에 $1/M$이 나오는 이유는 적분의 이산화로 쉽게 설명할 수 있다. 미분소 $d \omega$의 이산화는 $\Delta \omega$ = $\omega_0$ = $2 \pi / T$ = $2 \pi \mathbin{/} (M \Delta t)$ = $2 \pi \Delta f$이다. 이 결과를 식 (2)에 대입하고 $\Delta t$ = $1$로 두면, 최종적으로 $1/M$만 남는다. 디랙 델타 함수를 쓰지 않고 푸리에 변환에 나오는 적분을 $\Delta t$와 $\Delta f$ = $1 \mathbin{/}(M \Delta t)$의 합으로 근사화해서 DFT 쌍을 만들 수도 있다.

                  (6b)

여기서 $\omega_k$ = $k \Delta \omega$ = $2 \pi k \mathbin{/}(M \Delta t)$ = $2 \pi k \Delta f$, $t_m$ = $m \Delta t$이다.

식 (6)을 한 차원 더 확장해서 2차원 이산 푸리에 변환쌍을 정의한다.

                  (7)

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

                  (8a)

                  (8b)

여기서 $\delta_{ml}$은 크로네커 델타(Kronecker delta)이다.

식 (8)의 최종식은 등비 급수(geometric series)를 적용해서 간단히 얻을 수도 있다.

                  (9)

또한 식 (9)는 인수 분해 공식 $z_m^M - 1$ = $(z_m - 1)$$(1 + z_m + z_m^2 + \cdots + z_m^{M-1})$에 $z_m$ = $e^{i2 \pi m / M}$, $z_m^M$ = $1$을 대입해 0이 됨을 보일 수 있다. 여기서 $m \ne 0$이다. 서로 다른 모든 기저 항(basis term)을 더하면 0이 된다는 뜻인 식 (9)는 기저 벡터(basis vector)의 선형 독립(linear independence)을 직관적으로 보여준다. 즉, $[1~z_m~z_m^2~\cdots~z_m^{M-1}]$인 기저 벡터는 모든 성분이 1인 기저 벡터에 직교한다. 식 (9) 대신 식 (8)의 첫째식을 고려하면, $[1~z_m~z_m^2~\cdots~z_m^{M-1}]$과 $[1~z_{m'}~z_{m'}^2~\cdots~z_{m'}^{M-1}]$는 켤레 전치 행렬(conjugate transpose) 연산에서 직교하므로, 이 기저 벡터로 구성한 행렬은 복소 영역의 직교 행렬인 유니터리 행렬(unitary matrix)이 된다. 여기서 $z_m^k$ = $z_k^m$이다. 따라서 식 (6)을 행렬 곱으로 바꾼 결과에는 유니터리 행렬이 등장하기 때문에, 이 행렬의 역행렬은 켤레 전치 행렬로 쉽게 구한다. 다시 말해 식 (9)는 식 (6)이 유니터리 관계임을 말하고 있어서, DFT를 보는 또 다른 관점은 유니터리 연산이다.

푸리에 변환(Fourier transform)의 대표 연산인 길쌈(convolution)을 이산화해서 DFT용 이산 길쌈(discrete convolution)을 정의할 수 있다.

                  (10)

여기서 $f_n, g_n$의 지표(index) $n$이 가진 주기는 $M$; $f,g$는 모두 $f_{n+M}$ = $f_n$인 순환 편이(circular shift)를 만족한다. 길쌈의 DFT는 푸리에 변환과 마찬가지로 각 DFT의 곱이 된다.

                  (11a)

                  (11b)

여기서 $p_k$ = $2 \pi k \mathbin{/} M$이다. 식 (8b)를 사용해서 함수 곱의 DFT는 각 DFT의 길쌈으로 표현됨도 유도한다.

                  (12a)

                  (12b)

식 (11b)와 (12b)를  종합해서, IDFT에 대한 길쌈 관계식도 추가로 얻는다.

                  (13)

DFT의 지표(index) $m$은 항상 0부터 시작하지만, 경우에 따라 수열의 중앙 지표(central index)가 0이 되도록 $m$을 음수부터 출발시키는 중앙화 이산 푸리에 변환(centered DFT, CDFT)도 자주 볼 수 있다. 예를 들어, 식 (6)에서 $M$ = $M_e + M_o$로 나누고 지표를 $m',k'$로 바꾸어 CDFT를 정의한다.

                  (14)

여기서 $M_e$ = $\lfloor M/2 \rfloor$, $M_o$ = $\lfloor (M+1)/2 \rfloor$; $\lfloor x \rfloor$은 $x$를 넘지 않는 최대 정수를 구하는 바닥 함수(floor function)이다. 지표 $m,k$의 시작점을 임의로 택하는 일반화 이산 푸리에 변환(generalized DFT, GDFT)도 만들 수 있다. 표본 개수 $M$이 짝수일 때, $m'$ = $0$ 혹은 $k'$ = $0$을 기준으로 양 및 음의 지표에 대한 대칭성이 깨진다. 그래서 짝수 판별식(discriminant of even number) $\Delta_e(M)$를 써서 $M$이 짝수이든 홀수이든 항상 대칭이 되게 하는 대칭 이산 푸리에 변환(symmetric DFT, SDFT)을 도입할 수 있다.

                  (15)

여기서 $\Delta_e(M)$ = $1 - \operatorname{rem}(M,2)$, $\operatorname{rem}(m, n)$은 $m$을 $n$으로 나눈 나머지(remainder)이다.

이상에서 논의한 기본적인 이산 푸리에 변환을 바탕으로 더욱 다양한 이산 변환을 생성할 수 있다.

   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}$, $\Sigma'$은 으뜸 합산(prime summation) 기호[1], [2]로서 처음과 끝 항은 $1/2$만 더하라는 합 연산이다. 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.4)

여기서 $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$, $\Sigma'$은 으뜸 합산(prime summation)이다. 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)가 되게 한다.

DCT의 길쌈 $(f * g)_m^c$은 식 (10)에 우함수 조건인 식 (1.1)을 대입해서 새롭게 정의한다[3], [4].

                  (1.17a)

                  (1.17b)

여기서 $(f*g)_m$은 DFT의 길쌈, $f, g$의 주기는 $2M-2$, 급수에 곱해진 계수 $1/2$은 $X_k = F_k /2$에서 나온다. 그러면 길쌈 신호의 DCT는 DFT와 동일하게 계산된다.

                  (1.18a)

                  (1.18b)

여기서 $F_k, G_k$는 각각 $f_m, g_m$의 DCT이다.

   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)

여기서 $\Sigma'$은 식 (1.3c)에 도입한 으뜸 합산(prime summation)이다[1], [2]. 표본화 계수 $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$로 늘린다.

식 (1.17)에 정의한 DCT의 길쌈과 비슷하게 DST의 길쌈 $(f * g)_m^s$도 정의한다.

                  (2.10)

여기서 $f_m, g_m$은 식 (2.1)처럼 기함수로 계산한다. 식 (2.10)으로 표현한 길쌈의 변환은 DST가 아닌 DCT를 선택한다.

                  (2.11a)

                  (2.11b)

여기서 $f_m, g_m$의 DST는 $F_k, G_k$이다. DST 연산에 뜬금없이 DCT가 나오는 이유는 길쌈 $(f * g)_m^s$가 기함수 아닌 우함수 특성을 가지기 때문이다.

                  (2.12)

여기서 식 (2.1)에 따라 $g_{-m}$ = $-g_m$을 사용한다.

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

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

                  (3.1a)

                  (3.1b)

여기서 $x_m$ = $m \Delta x$, $p_k$ = $k \pi / [(M-1)\Delta x]$ = $k \Delta p$, $\Delta p\Delta x$ = $\pi /(M-1)$, $\alpha_{m, M}$ = $1+\delta_{m0} + \delta_{m,M-1}$, $\Sigma'$은 으뜸 합산(prime summation)[1], [2], $\alpha$는 $x_0$ 위치의 경계 조건을 결정하는 복소수인 혼합 계수(mixed coefficient)이다. 조건 $\Re[\alpha] \le 0$에 대한 DMFT의 역변환은 다음처럼 공식화된다.

             (3.2a)

             (3.2b)

혼합 푸리에 역변환(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.3)에 변형을 가해서 $\Re[\alpha] > 0$인 경우에도 이산화한 두 함수 사이에 직교성이 생기도록 $F_k$를 재정의한다.

                  (3.6)

여기서 $\operatorname{Sa}(x)$ = $\sin x \mathbin{/} x$, $F_0$ = $F_{M-1}$ = $0$, $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)

여기서 중점 근사에 따라 $f_0 \sin(p_k x_0)$ = $0$ $\approx$ $[f_1 \sin(p_k x_1) + f_{-1} \sin(p_k x_{-1})]/2$, $f_{M-1} \sin(p_k x_{M-1})$ = $0$ $\approx$ $[f_M \sin(p_k x_M) + f_{M-2} \sin(p_k x_{M-2})]/2$이다. 식 (3.7)은 근사식이지만 $M$이 충분히 크면 항등식에 가까워진다. 근사이지만 타당한 식 (3.7)의 최종 결과를 식 (3.6)에 대입해서 $f_m$에 대한 2계 차분 방정식(the second order difference equation)이 생기도록 한다.

                  (3.8)

만약 $M$을 충분히 키우면, 근사식 (3.8)을 등식처럼 생각할 수 있다. 또한 식 (3.8)은 정확히 혼합 푸리에 변환[= $\int_0^\infty [df(x)/dx + \alpha f(x)]\sin(px)\,dx$]의 이산화가 된다. 상이한 관점으로 식 (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.9)에 나온 $f_m$에 식 (3.10)을 대입한 후, 그 결과의 좌변과 우변을 비교해서 $A_k$를 구한다.

                  (3.11)

여기서 $m$ = $1,2,\cdots,M-2$이다. 따라서 식 (3.11)에서 구한 $A_k$를 이용해서 $f_{m,p}$를 완전히 공식화한다.

             (3.12)

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

                  (3.13)

여기서 양끝점인 $m$ = $0$과 $M-1$에서 $\Delta_s (r^m) \mathbin{/} (2\Delta x)$ = $-\alpha r^m$을 만족하며, $c_1, c_2$는 $m$과 무관한 경계 조건으로 계산하는 상수[간단하게 $f_0$과 $f_{M-1}$로 $c_1, c_2$를 결정 가능]이다. 식 (3.12)와 (3.13)을 합쳐서 이산 혼합 푸리에 역변환(inverse discrete mixed Fourier transform, IDMFT)을 확정한다.

                  (3.14)

여기서 일반해를 표현할 때 식 (3.13)의 $(-r_1)^{-m}$ 대신 $m$에 대해 대칭적인 $(-r_1)^{M-m-1}$을 쓴다.[이때 일반해가 달라 식 (3.13)과 (3.14)의 $c_2$는 같지 않다.] 연립 방정식 개수 관점에서 $c_1, c_2$의 필요성을 이해할 수 있다. 식 (3.9)로 기술되는 차분 방정식의 개수는 $M-2$개밖에 되지 않아, $M$개인 미지수 $f_m$을 명확히 결정하기에는 방정식 2개가 부족하다. 이 결함은 $f_m$ 중의 어느 2개를 알고 있다는 가정으로 해결하며, 이런 타개책을 수학에서는 일반해라 명한다. 

상수 $\alpha \Delta x$가 매우 작은 경우에 $r_1, r_2$를 지수 함수로 근사할 수 있다.

                  (3.15)

식 (3.15)에 따라 $r_1$을 감소하는 지수 함수처럼 간주할 수 있어서 식 (3.5)처럼 급수 합을 정의한다.

                  (3.16)

신기하게도 식 (3.5)와 다르게, 식 (3.16)은 합산이 항상 0이 나오는 직교성(orthogonality)이 성립한다.

[DMFT 기저 함수(basis function)의 직교성]

                          (3.17)

여기서 $k$ = $0, 1, 2, \cdots, M-1$, $\Sigma'$은 식 (1.3c)에 나온 으뜸 합산(prime summation) 기호이다.

[증명]

기저 함수를 삼각 함수의 합차 공식으로 분해한다.

                  (3.18a)

식 (3.18a)를 합산하고 첨자 $m$을 맞추어 정리한다.

                  (3.18b)

여기서 $m$ = $1, 2, \cdots, M-2$이다. 처음 항인 $m$ = $0$ 경우는 식 (3.18b)가 성립하지 않지만, 으뜸 합산의 정의[처음 항은 반만 더함]로 인해 $m$ = $0$의 결과는 $m$ = $1$에 대한 식 (3.18b)의 둘째 항이 됨으로써 $m$ = $1$의 합에 포함된다.

                  (3.18c)

마찬가지로 식 (3.17)의 끝 항인 $m$ = $M-1$은 $m$ = $M-2$ 합산의 셋째 항으로 간주된다. 또한 $r_2$ = $(-r_1)^{-1}$은 식 (3.13)의 또 다른 근이므로, 식 (3.18b)에 의해 식 (3.17)의 둘째 줄도 유도된다.

______________________________

DMFT의 기저 함수 직교성을 식 (3.14)에 적용하여 식 (3.3)과 비슷한 방식으로 $c_1, c_2$를 도출할 수 있다[1].

                  (3.19a)

                  (3.19b)

                  (3.19c)

모든 $f_m$에 대한 합산이 번거로울 때는 $f_0, f_{M-1}$을 안다는 조건을 사용해서 $c_1, c_2$를 풀어도 식 (3.19c)와 동일한 결과를 얻을 수 있다. 어떤 경우든지 $F_k$만을 연산해서 모든 $f_m$을 일의적으로 결정하기는 불가능하다. 따라서 DMFT는 수학적으로 완전한 변환이 아니고 결점을 수정해서 써야 하는 불완전 변환이다.

[참고문헌]

[1] G. D. Dockery and J. R. Kuttler, "An improved impedance-boundary algorithm for Fourier split-step solutions of the parabolic wave equation," IEEE Trans. Antennas Propag., vol. 44, no. 12, pp. 1592–1599, Dec. 1996.

[2] 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.

[3] S. A. Martucci, "Symmetric convolution and the discrete sine and cosine transforms," IEEE Trans. Signal. Process., vol. 42, no. 5, pp. 1038–1051, May 1994.

[4] V. G. Reju, S. N. Koh, and I. Y. Soon, "Convolution using discrete sine and cosine transforms," IEEE Signal Process. Lett., vol. 14, no. 7, pp. 445–448, Jul. 2007.

[다음 읽을거리]

1. 삼각 함수 보간

2. Z 변환

3. 고속 푸리에 변환