Lec 13. Discrete Fourier Transform

본 게시물은 University of Arizona의 Prof. J. Scott Tyo의 OPTI512R 강의를 보고 정리한 내용입니다.
전공자가 아닌 만큼 부정확한 정보가 포함되어 있을 수 있습니다.
부정확한 정보가 있다면 알려주시면 감사하겠습니다.

---

The Discrete Fourier Transform

Space function이 continuous 하고 Aperiodic 할 때, 그 Fourier Transform 또한 continuous 하고 Aperiodic 하다. 만약 Space function이 periodic 하면, 우리는 Frequency domain에서 discrete 한 Foureir Series를 가지게 된다. 마찬가지로 function이 space에서 discrete 할 때, function의 Fourier Transform이 Continuous 하고, Periodic 하다. 이는 아래 Table에 정리되어 있다.

	Continuous Space	Discrete Space
Continuous Frequency	Fourier Xform, Aperiodic	DSFT, Periodic in Frequency
Discrete Frequency	Fourier Series, space periodic	DFT

 

마지막으로 Space와 Frequency 모두 Discrete 한 경우가 남았다. 일반적으로 Space와 Frequency가 모두 Continuous 한 경우 function이 Aperiodic 하기 때문에 우리는 결과가 periodic 할 것이라고 예상할 수 있을 것이다. 우리의 이론이 성립할지 알아보도록 하자. Discrete space와 Continuous Frequency의 DSFT(Discrete Space Fourier Transform)의 경우 우리는 index를 $-\infty ~ \infty$로 설정할  수 있기 때문에, Series가 Infinite Series가 된다는 점에 주목했었다. $f(x) = 0$이 interval $0 \leq x < L$의 밖에 존재한다는 사실을 알고 있다면 어떤 일이 발생하는지 생각해보자. 이 interval에서 N개의 equally-spaced sample이 존재한다고 가정한다면 다음이 성립한다.

 

$$x_z = {L \over N}, kx_s = {kL \over N}, k = 0,1, .....,N-1$$

 

우리는 Band limited 하고 Space limited 한 signal을 규정하였지만, 이 명제는 엄밀히 말하면 참 일 수 없는 난제가 발생한다. function $f(x)$가 space limited 하기 때문에, 우리는 아래의 조건이 만족되는 경우 Whittaker-Shannon sampling theorem을 통해 $F(\xi)$를 $n\triangle \xi$에서 sample로 재구성할 수 있다는 것을 안다.

 

$${1 \over \triangle\xi} \geq L$$

 

따라서 우리는 Sampling Theorem을 만족하는 Frequency domain의 minimum interval로 $\triangle \xi = {1 \over L}$을 선택할 것이다. 이제 Discrete space Fourier Transform을 다음과 같이 Sampling 하도록 하겠다.

 

$$F_s(\xi) = \sum^{N-1}_{k=0} f_k e^{-j2\pi \xi k x_s} \\ F_s(n \triangle \xi ) = F_n = \sum^{N-1}_{k=0} f_k e^{-j2\pi {L \over N}n{1 \over L}} = \sum^{N-1}_{k=0} f_k e^{-j2\pi {kn \over N}}$$

 

위 수식 중 두 번째 수식은 Sequence $f_k, k = 0,1,......,N-1$의 Discrete Fourier Transform(DFT)를 정의한다.

다른 Fourier Transform들과 마찬가지로 DFT 또한 Inverse Transform을 가지고 있는데 이는 다음과 같다.

 

$$f_k = {1 \over N} \sum^{N-1}_{n=0} F_n e^{j2\pi {kn \over N}}$$

---

Discrete Exponential Sequences

DFT는 특정한 Speical sequence인 Exponential sequence로 작성된다. 다음과 같은 Sqeuence가 있다고 하자.

 

$$\psi_k[n] = e^{j2\pi {kn \over N}}$$

 

이 Sequence는 $k = 0, 1, ....., N-1$과 $n = 0,1,....,N-1$에 대해 define 되어 있다. 우리는 이 Sequence를 N element vector로 생각할 수 있고, dot product를 적용하여 다음과 같은 결과를 얻을 수 있다 

 

$$<\psi_k[n], \psi_m^*[n]> = \sum^{N-1}_{n=0} e^{j2\pi {kn \over n}}e ^{-j2\pi {mn \over N}} \\ = \sum^{N-1}_{n=0} e^{j{2\pi \over N}(k-m)n}$$

 

위 수식은 아래 아래 수식의 Tabulated series이다.

 

$$\sum^{N-1}_{n=0} x^n = {1 - x^N \over 1 -x}$$

 

우리의 경우에는 $x = e^{j {2\pi \over N}(k-m)}$이고,

 

$$\psi_k[n] \cdot\psi_m^*[n] = {1 - e^{j2\pi (k-m)} \over 1 - e^{j{2\pi \over N} (k-m)}}$$

 

이다. $k$와 $m$이 모두 integer이므로, $k \neq m$이면, numerator는 정확히 0이 된다. $k=m. x = 1$이고, series가 정확히 N일 때 dot product를 다음과 같이 정의한다.

 

$$\psi_k[n] \cdot\psi_m^*[n] = N \delta_{km}$$

 

$\delta_{km}$ function은 Kroneker delta function으로 우리가 알고 있는 dirac delta function에 대한 Discrete analog이고, 다음과 같이 정의된다.

 

$$\delta_{km} = \begin{cases}1 & k = m\\0 & k \neq m\end{cases}$$

---

Properties of DFT

DFT는 Continuous 한 사촌 격인 Fourier Transform과 많은 동일한 Property를 가지고 있다. 하지만, Fourier Transform과는 구분되는 것들 또한 존재한다.

 

1. DFT는 Linear 하고, 이는 definition을 통해 쉽게 알 수 있다. 즉, $\mathcal{DFT} \{f_k + g_k\} = F_n + G_n$가 성립한다.

 

2. DFT는 periodic 하다.

 

$$F_{n+N} = \sum ^{N-1}_{k=0} f_k e^{-j2\pi {nk \over N}} e^{-j2\pi {Nk \over N}} \\ = \sum ^{N-1}_{k=0} f_k e^{-j2\pi {nk \over N}} \cdot 1 \\ = \sum ^{N-1}_{k=0} f_k e^{-j2\pi {nk \over N}} \\ = F_n$$

 

따라서 DFT는 Discrete 하면서 동시에 Periodic 하다는 것을 알 수 있다.

 

3. 마찬가지로  IDFT 또한 Periodic 하다.

 

$$F_{k+N} = {1 \over N}\sum ^{N-1}_{n=0} F_n e^{j2\pi {kn \over N}} e^{j2\pi {kN \over N}} \\ = e^{j2\pi {kN \over N}}\sum ^{N-1}_{n=0} F_n e^{j2\pi {kn \over N}} \\ = \sum ^{N-1}_{n=0}F_n e^{j2\pi {kn \over N}} \\ = f_k$$

 

Frequency domain sequence $F_n$ 뿐만 아니라 Space domain sequence $f_k$ 또한 Discrete 하고 Periodic 하다. $k <0$ 또는 $k \geq N$에 대하여 sequence $f_k = 0$이라고 가정하였지만, IFDT의 proerty가 Periodictity를 강조한다는 것을 알 수 있다. 이는 추후에 다루게 될 Convolution에 대해서 의미를 가지고 있다.

 

4. Modulated Signal의 DFT에 는 다음과 같다.

 

$$\mathcal{DFT} \{f_k e^{j2\pi {n_0 \over N}k}\} = \sum^{N-1}_{k=0} f_k e^{-j{2\pi \over N}k(n-n_0)} \\ = F_{n-n_0}$$

 

만약 $n-n_0$ 가 $0 \leq n < N-1$의 range를 벗어나는 경우 DFT의 periodic property가 사용 될 수 있다.

DFT의 Symmetry로 인해, 우리는 또한 Shifted signal의 DFT를 다음과 같이 정의할 수 있다.

 

$$\mathcal{DFT} \{f_{k-k_0}\} = e^{-j{2\pi \over N}k_0 n}F_n$$

 

5. Kroneker delta function의  DFT는 다음과 같다.

 

$$\mathcal{DFT} \{\delta_{km}\} = \sum^{N-1}_{k=0} \delta_{km} e^{-j2\pi {kn \over N}} \\ = e^{-j2\pi {m \over N}n}$$

 

Continuous 한 경우와 마찬가지로, Delta function의 DFT는 Complex exponential이라는 것 을 알 수 있다. DFT와 IDFT의 Symmetry 때문에 다음이 성립한다.

 

$$\mathcal{DFT} \{e^{j2\pi {n_0 \over N}k}\} = N\delta_{nn_0}$$

---

Evaluating Fourier Transform and 2D Fourier Transform in Matlab

지금부터는 Matlab에서 fft command를 활용하여 Foureir Transform을 수행하는 방법들에 대해 다룰 것이다. 우선 다음과 같은 function이 존재한다고 가정하자.

 

$$f(x) = \cos(2\pi x) + 3 \sin(6\pi x)$$

 

우리는 위 function에 대해 Analytic 하게 Fourier Transform을 적용할 수 있고, 다음을 얻을 수 있다는 것을 알고 있다.

 

$$F(\xi) = {1\over 2} [\delta(\xi - \xi_0) + \delta(\xi + \xi_0)] +{3 \over 2j}[\delta(\xi-3\xi_0) - \delta(\xi+3\xi_0)]$$

 

편의를 위해 위 function을 아래와 같이 정의하도록 하겠다.(Tistory codeblock에서 matlab을 지원하지 않아 python으로 설정 후 작성하였음)

 

x = linspace(0,2,1001);x=x(1:end-1);
f = cos(2*pi*x) + 3*sin(6*pi*x);

 

function $f(x)$는 아래 그림과 같이 Variable $x$의 범위인 0에서 2 직전까지 실행된다.

function의 length와 DFT의 length는 같고, DFT의 highest frequency가 ${1 \over \triangle x}$바로 아래이며, $\triangle x$는 Sampling interval이다. Frequency vector를 형성하기 위해 우리는 다음의 command를 사용한다.

 

dx = x(2) - x(1);
xi = linspace(0,1/dx, length(x)+1);xi = xi(1:end-1);

이 Command는 $\triangle \xi = {1 \over (\triangle xN)}$인 Frequency vector를 생성한다. 여기서 $N$은 sample의 개수이고, $\xi$ vector의 highest frequency는 $(1 \over \triangle x) - \trinagle \xi$이다.

 

이제 우리는 signal에 DFT를 취할 준비가 되었고, 이는 Matlab의 fft command를 사용하여 간단하게 이루어진다.

 

F = fft(f)*dx;
plot(xi,abs(F))

결과는 아래 그림과 같다. magnitude를 정확하게 맞추기 위해 Fourier Transform에 $dx$를 multiply 하였다(DFT의 정의가 Differential을 포함하지 않기 때문이다). 

 

 

위 그림에서 볼 수 있는 Fourier Transform의 결과는 우리에게 유용한 정보를 거의 담고 있지 않지만, 몇 가지 중요한 정보를 포함하고 있다. DFT는 Periodic 하며, $0$에서 ${1 \over (\triangle x 2)}$ 사이의 frequency만이 우리가 보던 것에 해당한다는 것이다(실제로는 그렇지 않다고 함). Matlab fft command는 Positive frequency에 대해서만 Fourier Transform을 계산한다. 우리는 ${1 \over (\triangle x2)}$에서 ${1\over \triangle x}$로 data를 가져와서 그들을 negative frequency로 전환하고자 한다. function F에 대해 이러한 작업을 처리하는 command는 fftshift이다. frequency variable에 대한 작업은 우리가 직접 해 주어야하기 때문에 다음의 command를 사용하였다.

plot([(xi(ceil(end/2)+1:end)-1/dx) xi(1:ceil(end/2))],fftshift(abs(F)))
set(gca,’xlim’,[-5 5])

위 Command는 첫 번째로, $\xi$ vector의 위쪽을 $1 \over \triangle x$ 만큼 아래로 shift 한 다음, 아래쪽 절반을 append 하는 vector를 만든다. command의 ceil(end/2) 부분은 {N \over 2}에서 round up 하여 command가 odd, even length vector 모두에 대해 작동하도록 보장한다. 마지막으로 우리는 $-5 \leq \xi \leq 5$로 axis를 limit 하여 우리가 관심 있는 정보를 확인할 수 있게 하였고, 그 결과는 아래 그림과 같다.

 

왼쪽의 그림은 magnitude를 나타내며, $\xi = \pm 1$과 $\xi = \pm 3$에서 4개의 $\delta$-function이 존재하는 것을 확인할 수 있다. $\delta-$function의 magnitude는 3:1의 ratio를 가질 것으로 예상되었다. 그러나 amplitude는 우리가 예상한 것보다 2배 큰 값을 가진다. 이러한 결과의 이유는 우리가 two period에 걸쳐 periodic function을 summation 했기 때문이다. 실제 Foureir Transform은 모든 period에 걸쳐 summation을 수행하며, magnitude는 infinite 하다. 또한 $\triangle \xi = {1 \over 2}$이기 때문에 이러한 impulse가 right area를 가지고 있음을 알 수 있다.

 

이제 오른쪽의 phase angle로 시선을 돌려보자. 우리는 곳곳에 angle이 존재하는 것을 알 수 있다. 그러나 이러한 phase anlge의 대부분은 $|F(\xi) \approx 0|$인 Frequency에 해당하므로 무의미하다. 우리가 신경 써야 할 4개의 Frequency에서 우리는 강조된 것과 같이 적절한 phase angle을 확인할 수 있다.

---

Frequency Leakage

우리가 학습한 DFT의 중요한 property 중 하나는 finite length sequence result의 DFT computation이 space $(k)$와 Frequency $(n)$ domain 모두에서 Periodicity를 enforce 한다는 것이다. 이러한 가정된 Periodicity는 중요한 결과를 가져온다.

x = linspace(0,3,101);
f = cos(2*pi*x) + .5*cos(2*pi*5*x);
dx = x(2) - x(1);
xi = linspace(0,1/dx,length(x)+1);xi = xi(1:end-1);
F = fft(f);

이 sequence의 DFT 결과는 아래 그림의 빨간색 선과 같다.

이 Sequence에는 $x=0$에서 $x=3$까지 실행되는 $N=101$개의 sample이 존재한다. Periodicity property를 고려하면 102번째 sample은 1번째 sample과 동일한 값을 가진다는 것을 알 수 있다. 이는 우리가 Fourier Transform을 적용하려 했던 original function이 아니다. 우리는 101번째 sample을 버려야 다음 sample이 1번째 sample $(x=0)$과 같은 original function을 얻을 수 있다. 이를 위해 아래의 command를 사용하였고, 그 결과는 위 그림의 파란색 선과 같다.

x = linspace(0,3,101);x = x(1:end-1);
f = cos(2*pi*x) + .5*cos(2*pi*5*x);
dx = x(2) - x(1);
xi = linspace(0,1/dx,length(x)+1);xi = xi(1:end-1);
F = fft(f);

---

Correcting for the Phase

DFT의 정의에 따르면 다음이 성립한다.

 

$$F_n = \sum^{N-1}_{k=0} f_k e^{-j{2\pi \over N}nk}$$

 

Sum은 $0$에서 $N-1$까지 이루어진다. 그러나 대부분의 경우 Fourier Transform을 적용하는 continuous function $f(x)$는 $-x_a < x < x_b$에서 정의된다. 즉, 첫 번째 sample이 $x=0$과 일치하지 않음을 의미한다. 종종 DFT compuation에서는 이를 해결하기 위해 shift를 수행한다. 

 

우리는 Lecture 8에서 $\mathcal{DFT} \{f[k-k_0]\} = e^{-j2\pi x_a \xi}F(\xi)$ 임을 배웠다. 따라서, $k=0$ sample이 $x=0$에 해당하지 않은 다는 사실을 설명하기 위해 Complex Exponential을 곱해야 한다.

 

x = linspace(-10,10,5001);x = x(5:end); %make selection asymmetric
dx = x(2)-x(1);
xi = linspace(1,1/dx,length(x)+1); xi = xi(1:end-1);
n = 0:(length(xi)-1); %these are the indices in the FD
f = sinc(x).^2;
k0 = x(1)/dx; %this finds the offset index in the SD
C = exp(-j*2*pi*k0*n/length(x)); %correction factor
F1 = fft(f);
F2 = fft(f).*C;

0에 대해 symmetric 하지 않은 $-2.45 \leq x \leq 2.5$ 범위에 대해 function을 정의하였다. complex eponential을 multiply 함으로 써 우리는 정확한 phase를 얻게 된다. F1과 F2의 real 및 imaginary part는 아래 그림에 표시되어 있다. 0 sample이 실제로 x vector에 포함되지 않은 경우에는 phase를 설명하기 어렵기에 주의해야 한다.

 

---

Defining Frequnecy Domain Filters

우리는 pure real function $h(x)$가 Hermitian symmetry, 즉 $H(-\xi) = H^*(\xi)$를 가지는 Fourier Transform을 가진다는 것을 배웠다. Matlab에서 filter를 frequency domain에 구성할 경우 이 symmetry가 유지되는지 확인해아 한다. 이를 확인하는 방법에는 여러 가지 방법이 존재한다.

 

예를 들어 $x_0=.25$이고 $h(x) = \delta(x-x_0)$인 filter를 생성한다고 할 때, 우리는 analytical 하게 이 function이 다음과 같은 Foureir Transform을 가진다는 것을 알고 있다.

 

$$H(\xi) = e^{-j2\pi \xi x_0}$$

 

하지만 우리의 상황에서는 이것을 우리의 function의 range에 있는 frequency에 대해 적절하게 구성해야 하고 이는 다음과 같은 command로 수행된다.

 

 H = exp(-j*2*pi*.25*[xi(1:ceil(end/2)) (xi(ceil(end/2)+1:end)-1/dx)]);

 

이 command는 symmetry를 올바르게 하기 위해 frequency의 second half를 shifted version으로 replace 한다. 아래 그림에서는 original function f와 shifted version real(ifft(H.*F))를 보여준다. 항상 약간의 residual imaginary part가 존재하기 때문에 real part에 ifft를 적용해야 한다. 이 경우 ifft의 imaginary part의 maximum value는 machine precision인 $10^{-15}$이다.