2020년 9월 20일 일요일

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

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


이산 푸리에 변환(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$로 제한, $\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)을 다음처럼 표현할 수 있다.

                  (6)

이산 푸리에 역변환(inverse discrete Fourier transform, IDFT)에 $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)의 첫째식을 둘째식에 대입하거나 둘째식을 첫째식에 대입하면, 다음 항등식이 유도된다.

                  (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$이다. 길쌈의 DFT는 푸리에 변환과 마찬가지로 각 DFT의 곱이 된다.

                  (11a)

                  (11b)

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

                  (12a)

                  (12b)

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

                  (13)

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


   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}$는 특성 방정식을 풀어서 정한다.

                  (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 기저 함수의 직교성]

                          (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.

[다음 읽을거리]

댓글 16개 :

  1. 식 16에서 M이 샘플의 개수인가요? 원래 신호의 샘플의 개수가 N개 이면 M=N이 되는건가요?
    아니면.. 원래 신호를 대칭으로 만들어서 M이 2N-1이 되어야 하나요? ㅠㅠ

    답글삭제
    답글
    1. 1. 네, 맞습니다. 실제 얻은 표본의 개수입니다. [그림 3]에서는 빨간색으로 표시되어 있어요.

      2. 본문에도 있듯이, M개의 실제 표본으로 DCT에 맞는 계산 표본을 적절히 만들어요. [그림 3]에 있는 검정색이 계산 표본입니다.

      삭제
    2. 답글 감사합니다 거북이님! 한가지 더 궁금한게 있는데요. 식 15를 이용하여 F(k) 를 구하면 F(0)은 DC 값을 의미하는건 알겠는데..
      F(1), F(2), ... F(M-1)은 무엇을 의미하나요? 그림 그려서 보면.. 각 파수에 대한 계수가 아닌거 같아서요..

      삭제
    3. 식 (14)와 같은 표본화된 신호 $f_m$에 대한 DFT입니다. 이걸 깔끔한 DCT-II 공식으로 바꾸는 방법을 설명하고 있어요.

      삭제
  2. 감사합니다 천천히 읽어보겠습니다

    답글삭제
  3. 안녕하십니까 선생님. 궁금한 점이 있어 댓글을 달게 되었습니다.
    (4)번 수식의 두번째 줄에서 m에 대한 sigma를 '1의 CTFS 역변환'으로 해석해서 보려고 하는데, 2pi가 어디에서 튀어나온 것인지 궁금합니다.
    더불어 마지막 줄에서 T = M * delta t일텐데, delta t를 1로 둔 이유와 2pi는 어디로 사라졌는지 궁금합니다.

    답글삭제
    답글
    1. 선생님 말고, 편하게 전파거북이로 불러주세요.

      식 (4) 밑에 내용을 추가했어요. 핵심은 디랙 델타 함수를 푸리에 급수로 바꾼 겁니다.

      삭제
    2. 감사합니다! 전파거북이님!

      삭제
  4. 안녕하세요, 전파거북이님. 저는 Imperial College London에서 전자공학을 전공중인 학부생입니다. 전파거북이님 덕분에 Signals & Systems 강의에 한번도 출석하지 않고도 A*를 받았습니다. 진심으로 감사합니다.

    답글삭제
    답글
    1. 축하드리고 대단하네요, 익명님. 출석 않고 A 받기는 매우 어려운데요 ㅋㅋ

      삭제
    2. 출석보다 이 글 읽는게 더 도움 되네요 ㅎㅎ

      삭제
  5. 전파거북이님, 4b 식 1번줄에서 2번줄로 넘어갈때 F_k의 periodicity를 고려하고 sigma k from -inf to inf k/M = sigma n from -inf to inf (n) + sigma k from 0 to M-1 (k/M)이라는걸 고려했을때 algebraic하게 정리되는것으로 이해하고 따라가고 있는데요, 이 외에 다른 시각적인/푸리에적인 관점으로 이 algebraic manipulation을 바라볼 수 있을까요?

    답글삭제
    답글
    1. DFT 관계식은 유니터리 행렬을 구성해서, 식 (9) 밑에 추가 설명을 붙였습니다.

      삭제
    2. 4b식 전개에 대한 제 이해는 틀리지 않았을까요? F_k는 M-1 이후부터 반복되니 저렇게 써주고, exp(k/M) infinite series가 두 series로 분해된다고 보면 될까요? k/M = (k + n/M) where k ]-oo,oo[ , n[0,M-1]

      삭제
    3. 같은 지표(index) $k$를 써서 헷갈릴 수 있는데요, 나눗셈 원리에 따라 $k_1$ = $mM + k_2$로 분해해서 계산했어요. 여기서 $k_1$, $k_2$는 각각 첫째와 둘째 줄의 지표입니다.

      삭제

욕설이나 스팸글은 삭제될 수 있습니다. [전파거북이]는 선플운동의 아름다운 인터넷을 지지합니다.