Discrete Fourier Transform

Introduction

Discrete Fourier transform is one of the most important linear transformation for signal processing. Although I have learned discrete Fourier transformation in college and graduate school, it was not systematic at that time. Because I did not have too much chance to actually use it without having to be dependent on third-party library or software, I only remember it is a transformation between the time domain and the frequency domain without being able to tell too much mathematical details.

In this blog post, I would like to discuss some of the details of discrete Fourier transform and its relationship with Fourier transform.

Fourier Transform

In mathematics, a Fourier transform (FT) is a mathematical transform that decomposes a function (often a function of time, or a signal) into its constituent frequencies, which is sometimes called frequency spectrum. The constituent frequencies could also be used to reconstruct the function back.

To decompose a function $f(x)$ to its frequency spectrum function $F(\xi)$ or $F(\omega)$, we use Fourier transform.

$$
F(\xi) = \int_{-\infty}^{\infty} f(x) e^{-2\pi i \xi x} dx
$$

where $\xi$ is any real-valued frequency.

Sometimes, if we define the radian-frequency $\omega = 2\pi \xi$,

$$
F(\omega) = \int_{-\infty}^{\infty} f(x) e^{-i \omega x} dx
$$

To reconstruct the function $f(x)$ from its frequency spectrum function $F(\xi)$ or $F(\omega)$, we use inverse Fourier transform.

$$
f(x) = \int_{-\infty}^{\infty} F(\xi) e^{2\pi i \xi x} d\xi
$$

or equivalently,

$$
f(x) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega x} d\omega
$$

Discrete Fourier Transform

Discrete Fourier transform replaces the infinite integral in the Fourier transform with a finite sum. Given a signal sequence $\{f(x_0), f(x_1), \cdots, f(x_{n-1})\}$, where the sampling interval between $x_i$ and $x_{i+1}$ is a constant, we defined the discrete Fourier transform, and got a frequency sequence $\{F(\omega_0), F(\omega_1), \cdots, F(\omega_{n-1})\}$.

$$
F(\omega_k) = \sum_{n = 0}^{N-1} f(x_n) e^{-i \omega_k x_n}
$$

where

$k = 0, 1, 2, \cdots, N-1$.

$N$ is the number of samples from $f$ and $F$.

$f(x_n)$ is the signal amplitude at $x_n$. It could be real or complex value.

$T$ is the sampling interval and $T = x_{i+1} - x_{i}$.

$x_n = nT$.

$F(\omega_k)$ is the spectrum of $f$ at radian-frequency $\omega_k$.

$f_0 = \frac{1}{NT}$ is the fundamental frequency, because we treat the $N$ signal examples as if they are from one period.

$\Omega = \frac{2\pi}{NT}$ is the radian-frequency sampling interval.

$\omega_k = k\Omega = \frac{2\pi k}{NT}$ is the $k$th radian-frequency sample.

$f_s = \frac{1}{T}$ is the sampling rate.

Similarly, the inverse discrete Fourier transform could be derived based on the definition of discrete Fourier transform.

$$
f(x_n) = \frac{1}{N} \sum_{k = 0}^{N-1} F(\omega_k) e^{i \omega_k x_n}
$$

Sometimes, the sampling interval $T$ is not important and we are interested in the magnitude of the signal and frequency samples. Given a signal sequence $\{f[0], f[1], \cdots, f[N-1]\}$ and we got a frequency sequence $\{F[0], F[1], \cdots, F[N-1]\}$ after fourier transformation. We have the more general form of discrete Fourier transform.

$$
\begin{align}
F[k] &= \sum_{n = 0}^{N-1} f[n] e^{-i \omega_k x_n} \\
&= \sum_{n = 0}^{N-1} f[n] e^{-i \frac{2\pi k}{NT} nT} \\
&= \sum_{n = 0}^{N-1} f[n] e^{-i \frac{2\pi k n}{N}} \\
\end{align}
$$

Similarly, we have the inverse discrete Fourier transform.

$$
\begin{align}
f[n] &= \frac{1}{N} \sum_{k = 0}^{N-1} F[k] e^{i \omega_k x_n} \\
&= \frac{1}{N} \sum_{k = 0}^{N-1} F[k] e^{i \frac{2\pi k}{NT} nT} \\
&= \frac{1}{N} \sum_{k = 0}^{N-1} F[k] e^{i \frac{2\pi k n}{N}} \\
\end{align}
$$

The reason why we only need $N$ frequency samples becomes obvious here. Because

$$
\begin{align}
F[k+N] &= \sum_{n = 0}^{N-1} f[n] e^{-i \frac{2\pi (k+N) n}{N}} \\
&= \sum_{n = 0}^{N-1} f[n] e^{-i \frac{2\pi k n}{N}} e^{-i 2\pi n} \\
&= \sum_{n = 0}^{N-1} f[n] e^{-i \frac{2\pi k n}{N}} \times 1\\
&= \sum_{n = 0}^{N-1} f[n] e^{-i \frac{2\pi k n}{N}}\\
&= F[k] \\
\end{align}
$$

Because the frequency has a “period” of $N$, only $N$ frequency examples are required.

The frequencies corresponding to the frequency samples $F[k]$ is $k f_0$, i.e., the frequencies are $0, \frac{1}{NT}, \frac{2}{NT}, \cdots, \frac{N-1}{NT}$ for $F[0], F[1], F[2], \cdots, F[n]$, respectively. Similarly the radient frequency corresponding to the frequency samples $F[k]$ is $2 \pi k f_0$.

Unitary Discrete Fourier Transform

Because discrete Fourier transform is a linear transformation, we could express such transformation using matrix multiplications.

Surprisingly, the $e^{i \frac{2\pi k}{N}}$ for $k \in [0,N-1]$ in the inverse discrete Fourier transform are actually the $N$th root of unity, $\omega_N^{0}, \omega_N^{1}, \cdots, \omega_N^{N-1}$, i.e., $\omega_N^{k} = e^{i\frac{2\pi k}{N}}$.

So we could write the inverse discrete Fourier transform using the following expression.

We define the inverse discrete Fourier transform Vandermonde matrix $\mathbf{DFT}^{\ast\dagger}$ as

$$
\begin{align}
\mathbf{DFT}^{\ast\dagger} &=
\begin{bmatrix}
\omega_N^{0 \cdot 0} & \omega_N^{0 \cdot 1} & \cdots & \omega_N^{0 \cdot N-1} \\
\omega_N^{1 \cdot 0} & \omega_N^{1 \cdot 1} & \cdots & \omega_N^{1 \cdot N-1} \\
\vdots & \vdots & \ddots & \vdots \\
\omega_N^{N-1 \cdot 0} & \omega_N^{N-1 \cdot 1} & \cdots & \omega_N^{N-1 \cdot N-1} \\
\end{bmatrix}
\end{align}
$$

where $\omega_N^{i \cdot j} = (\omega_N^{i})^j = \omega_N^{ij}$.

We further arrange the signal and frequency samples as column vectors.

$$
F_N =
\begin{bmatrix}
F[0] \\
F[1] \\
\vdots \\
F[N-1] \\
\end{bmatrix}
$$

and

$$
f_N =
\begin{bmatrix}
f[0] \\
f[1] \\
\vdots \\
f[N-1] \\
\end{bmatrix}
$$

Therefore,

$$
\begin{align}
f_N &= \frac{1}{N} \mathbf{DFT}^{\ast\dagger} F_N \\
F_N &= \mathbf{DFT}^{\ast} f_N \\
\end{align}
$$

However, $\mathbf{DFT}^{\ast\dagger}$ is not unitary.

$$
\mathbf{DFT}^{\ast\dagger}_{j,k} = \omega_N^{jk}
$$

$$
\mathbf{DFT}^{\ast}_{j,k} = \omega_N^{-kj}
$$

Therefore,

$$
\begin{align}
(\mathbf{DFT}^{\ast} \mathbf{DFT}^{\ast\dagger})_{j,k} &= \sum_{i=0}^{N-1} \omega_N^{ji} \omega_N^{-ik} \\
&= \sum_{i=0}^{N-1} \omega_N^{-i(k-j)} \\
&=
\begin{cases}
\sum_{i=0}^{N-1} 1 & \text{when $k = j$}\\
\omega_N^{j-k} \sum_{i=0}^{N-1} \omega_N^{i} & \text{when $k \neq j$}\\
\end{cases} \\
&=
\begin{cases}
N & \text{when $k = j$}\\
\omega_N^{j-k} \frac{1 (1 - \omega_N^N)}{1-\omega_N^1} & \text{when $k \neq j$}\\
\end{cases} \\
&=
\begin{cases}
N & \text{when $k = j$}\\
\omega_N^{j-k} \frac{1 (1 - 1)}{1-\omega_N^1} & \text{when $k \neq j$}\\
\end{cases} \\
&=
\begin{cases}
N & \text{when $k = j$}\\
0 & \text{when $k \neq j$}\\
\end{cases} \\
\end{align}
$$

It is not unitary, but we are close. We define a new matrix

$$
\mathbf{DFT}^{\dagger} = \frac{1}{\sqrt{N}} \mathbf{DFT}^{\ast\dagger}
$$

This time, it is easy to show that this time $\mathbf{DFT}^{\dagger}$ is a unitary matrix.

To conduct discrete Fourier transform and inverse discrete Fourier transform, we have

$$
\begin{align}
f_N &= \frac{1}{\sqrt{N}} \mathbf{DFT}^{\dagger} F_N \\
F_N &= \sqrt{N} \mathbf{DFT}^{} f_N \\
\end{align}
$$

The unitary discrete Fourier transform and the inverse unitary discrete Fourier transform accordingly are

$$
\begin{align}
f^{\prime}_N &= \mathbf{DFT}^{\dagger} F_N \\
F^{\prime}_N &= \mathbf{DFT}^{} f_N \\
\end{align}
$$

Examples

Here I implemented discrete Fourier transform, inverse discrete Fourier transform, unitary discrete Fourier transform, and inverse unitary discrete Fourier transform using unitary discrete Fourier transform matrix multiplications. The correctness is also verified using the Numpy Discrete Fourier Transform library.

dft.py

import numpy as np

def get_inverse_unitary_dft_vandermonde_matrix(N):

    def vandermonde_row(k):
        omega = np.exp(2 * np.pi * k / N * 1j)
        row = [1]
        for _ in range(1, N):
            row.append(omega * row[-1])
        return row

    matrix = []
    for row in range(N):
        matrix.append(vandermonde_row(k=row))

    matrix = np.matrix(matrix)
    matrix = 1 / np.sqrt(N) * matrix

    return matrix

def dft(f_N):

    N = len(f_N)
    matrix = get_inverse_unitary_dft_vandermonde_matrix(N)
    matrix = matrix.getH()
    F_N = np.sqrt(N) * np.matmul(matrix, f_N)

    return F_N

def inverse_dft(F_N):

    N = len(F_N)
    matrix = get_inverse_unitary_dft_vandermonde_matrix(N)
    f_N = 1 / np.sqrt(N) * np.matmul(matrix, F_N)

    return f_N

def unitary_dft(f_N):

    N = len(f_N)
    F_N = 1 / np.sqrt(N) * dft(f_N)

    return F_N

def unitary_inverse_dft(F_N):

    N = len(F_N)
    f_N = np.sqrt(N) * inverse_dft(F_N)

    return f_N

def main():

    random_seed = 0
    np.random.seed(random_seed)

    N = 4

    f_N = np.random.uniform(low=-10, high=10, size=N)
    F_N = np.random.uniform(low=-10, high=10, size=N)

    # Verify the correctness of DFT implementations
    assert np.allclose(dft(f_N=f_N), np.fft.fft(f_N)), "The DFT implementation result does not match the Numpy implementation result!"
    assert np.allclose(inverse_dft(F_N=F_N), np.fft.ifft(F_N)), "The inverse DFT implementation result does not match the Numpy implementation result!"
    # Unitary DFT
    # https://numpy.org/doc/stable/reference/routines.fft.html#normalization
    assert np.allclose(unitary_dft(f_N=f_N), np.fft.fft(f_N, norm="ortho")), "The unitary DFT implementation result does not match the Numpy implementation result!"
    assert np.allclose(unitary_inverse_dft(F_N=F_N), np.fft.ifft(F_N, norm="ortho")), "The unitary inverse DFT implementation result does not match the Numpy implementation result!"

if __name__ == "__main__":
    
    main()

References

• Complex Number Representations
• Introduction to the DFT