siware.dev

Discrete Fourier Transform

The idea of Fourier is that a function (signal) can be decomposed as infinite number of sine/cosine waves (harmonics). We can consider sinsin and coscos as the basis function and they form an orthogonal basis.

We can lookup on the Internet regarding the topic of Fourier and we will find there are extensive knowledge. Studying the terminology is important to avoid confusion. For starter, the concept of Fourier manifested into many variants:

The differences can be summarized in this table:

Series Input Periodicity Input Resolution Output Periodicity Output Resolution
CT-DFS Infinite Periodic Continuous Aperiodic Discrete
DT-DFS Infinite Periodic Discrete Periodic Discrete
CTFT Infinite Aperiodic Continuous Aperiodic Continuous
DTFT Infinite Aperiodic Discrete Periodic Continuous
DFT Finite Aperiodic Discrete Periodic Discrete

Note that FFT is an implmentation of DFT and practically what we use for computation. For more details on the differences, consult http://complextoreal.com/wp-content/uploads/2013/04/fft5.pdf.

Fourier Transform Pair

Forward Fourier Transform F{f(x)}\mathcal{F} \{f(x)\} is mathematical operation that transforms a function of time (or space) into a function of frequency. Backward (Inverse) Fourier Transform F1{f^(ω)}\mathcal{F}^{-1}\{\hat{f}(\omega)\} reconstructs the function of time from function of frequency. Together they form Fourier Transform pair.

Fourier Transform is commonly expressed in terms of complex exponentials (instead of sine/cosine). Recall Euler’s Formula is eix=cos(x)+i  sin(x)e^{ix} = cos(x) + i \; sin(x).

In addition, since:

eimx,einxdx=ππeimxeinxdx={0mn2πm=n\quad \langle e^{imx}, e^{-inx} dx \rangle = \displaystyle\int_{-\pi}^{\pi} e^{imx} e^{-inx} dx = \begin{cases} 0 & \quad m \neq n \\ 2\pi & \quad m = n \end{cases}

ei2πmx,ei2πnxdx=ππei2πmxei2πnxdx={0mn1m=n\quad \langle e^{i 2\pi mx}, e^{-i 2\pi nx} dx \rangle = \displaystyle\int_{-\pi}^{\pi} e^{i 2\pi mx} e^{-i 2\pi nx} dx = \begin{cases} 0 & \quad m \neq n \\ 1 & \quad m = n \end{cases}

This makes chosing ei2πξxe^{i 2\pi \xi x} as basis function convenient. Note that ξ\xi is used to denote an ordinary frequency (Hz) whereas ω=2πξ\omega = 2 \pi \xi is for angular frequency (rad/s).

Continuous Time Fourier Transform (CTFT) pair is given as:

In Discrete Time Fourier Transform (DTFT), the function f(x)f(x) is sampled with period TsT_s (freq 1ξs\frac{1}{\xi_s}) and this can be expressed as:

fs(x)=n=f(nTs)  δ(xnTs)\quad f_s(x) = \displaystyle\sum_{n=-\infty}^{\infty} f(nT_s) \; \delta(x - nT_s)

plugging fsf_s to CTFT forward transform yields:

f^(ξ)=n=f(nTs)  ei2πξnTs\quad \hat{f}(\xi) = \displaystyle\sum_{n=-\infty}^\infty f(nT_s) \; e^{-i 2\pi \xi nT_s}

and let the n-th sample fn=f(nTs)f_n = f(nT_s) and ξ^=ξ  Ts=ξ/ξs\hat{\xi} = \xi \; T_s = \xi / \xi_s, this simplifies to:

f^(ξ^)=n=fn  ei2πξ^n\quad \hat{f}(\hat{\xi}) = \displaystyle\sum_{n=-\infty}^{\infty} f_n \; e^{-i 2 \pi \hat{\xi} n}

Note that the output of CTFT and DTFT are continuous function of frequency. In addition DTFT requires infinite sum which is not suitable for computation. Discrete Fourier Transform (DFT) approximates the result by sampling the output function.

Given LL samples of a function ff and f^\hat{f} is discretized by NN, DFT can be expressed as:

fk^=n=0L1fn  ei2πNkn\quad \hat{f_k} = \displaystyle\sum_{n=0}^{L-1} f_n \; e^{-i \frac{2\pi}{N} k n}

Note: in most implementation L=NL = N.

To summarize, Fourier Transform pairs are:

Forward Transform Backward (Inverse) Transform
CTFT f^(ξ)=f(x)  ei2πξxdx\hat{f}(\xi) = \displaystyle\int_{-\infty}^\infty f(x) \; e^{-i 2\pi \xi x} dx f(x)=f^(ξ)  ei2πξxdξf(x) = \displaystyle\int_{-\infty}^\infty \hat{f}(\xi) \; e^{i 2\pi \xi x} d\xi
DFT fk^=n=0N1fn  ei2πNkn\hat{f_k} = \displaystyle\sum_{n=0}^{N-1} f_n \; e^{-i \frac{2\pi}{N} k n} fn=1Nk=0N1fk^  ei2πNknf_n = \frac{1}{N} \displaystyle\sum_{k=0}^{N-1} \hat{f_k} \; e^{i \frac{2\pi}{N} k n}

Properties of Fourier Transform

Spectral Derivative

The benefit of doing computation in Fourier space is that it is easy to compute derivative (a.k.a Spectral derivative).

    F{ddxf(x)}=dfdxeiωxdx\mathcal{F}\{\frac{d}{dx} f(x)\} = \displaystyle\int_{-\infty}^\infty \frac{df}{dx} e^{-i \omega x} dx

    =f(x)(iωeiωx)dx+[f(x)eiωx] = -\displaystyle\int_{-\infty}^\infty f(x) (-i \omega e^{-i \omega x})dx + \Bigl[ f(x) e^{-i \omega x}\Bigr]_{-\infty}^\infty

    =iωf(x)eiωxdx = i \omega \displaystyle\int_{-\infty}^\infty f(x) e^{-i \omega x} dx

    =iωF{f(x)} = i \omega \mathcal{F}\{f(x)\}

It turns out that computing FFT-based spectral derivative isn’t as straightforward as above. For complete discussion on this topic, please consult https://math.mit.edu/~stevenj/fft-deriv.pdf.

In summary, algorithm to compute first derivative:

  1. Given fnf_n for 0n<N0 \le n \lt N, use an FFT to compute fk^\hat{f_k} for 0k<N0 \le k \lt N
  2. Obtain fk^\hat{f_k}^\prime by multiplying {i2πNkk<N/2i2πN(kN)k>N/20k=N/2  (N  is  even)\begin{cases} i \frac{2 \pi}{N} k & \quad k \lt N/2 \\ i \frac{2 \pi}{N} (k - N) & \quad k \gt N/2 \\ 0 & \quad k = N/2 \; (N \; is \; even) \end{cases}
  3. Compute fnf_n^\prime from fk^\hat{f_k}^\prime via inverse FFT

To compute second derivative:

  1. Given fnf_n for 0n<N0 \le n \lt N, use an FFT to compute fk^\hat{f_k} for 0k<N0 \le k \lt N
  2. Obtain fk^\hat{f_k}^{\prime\prime} by multiplying {(2πNk)2kN/2(2πN(kN))2k>N/2\begin{cases} -(\frac{2 \pi}{N} k)^2 & \quad k \le N/2 \\ -(\frac{2 \pi}{N} (k-N))^2 & \quad k \gt N/2 \end{cases}
  3. Compute fnf_n^{\prime\prime} from fk^\hat{f_k}^{\prime\prime} via inverse FFT

Convolution

Convolution is defined as (fg)=f(xξ)g(ξ)dξ(f * g) = \displaystyle\int_{-\infty}^{\infty} f(x - \xi) g(\xi) d\xi

Fourier Transform property for convolution is that:

    F{(fg)}=F{f}F{g}=f^g^\mathcal{F}\{(f * g)\} = \mathcal{F}\{f\} \mathcal{F}\{g\} = \hat{f} \hat{g}

This basically means convolution equals to multiplication in Fourier space. This is another benefit of Fourier Transform.

Parseval’s Theorem

The energy in ff is proportional with energy in its fourier transform f^\hat{f}

    f(ω)^2dω=2πf(x)2dx\displaystyle\int_{-\infty}^{\infty} \lvert \hat{f(\omega)} \rvert^2 d\omega = 2 \pi \displaystyle\int_{-\infty}^{\infty} \lvert f(x) \rvert^2 dx

This theorem is useful because it verifies that if we truncate the small Fourier coefficients, we are still capturing most of function ff.

Linear Operator

It’s useful to know that Fourier Transform is a linear operator:

    F{αf(x)+βg(x)}=αF{f}+βF{g}\mathcal{F}\{\alpha f(x) + \beta g(x)\} = \alpha \mathcal{F}\{f\} + \beta \mathcal{F}\{g\}

References