Discrete Fourier Transform

# Discrete Fourier Transform The discrete Fourier transform (DFT) is a mathematical procedure used to determine the harmonic or frequency content in a discrete time signal sequence. The DFT originates from the continuous Fourier transform, $X(f)$, defined as follows: $$X(f) = \int_{-\infty}^{\infty} x(t) e^{-j 2 \pi f t}\,dt \ $$ where $x(t)$ is some continuous time-domain signal.

The DFT is defined as the discrete frequency-domain sequence, $X(m)$, defined as follows: $$\begin{aligned} X(m) & = \sum_{n=0}^{N-1} x(n) e^{-j \, 2 \pi \frac{n \, m}{N}} \quad \quad &\text{Exponential Form} \\ \\ X(m) & = \sum_{n=0}^{N-1} x(n) \left[cos \left( 2 \pi \, \frac{n \, m}{N} \right) \, - \, j \, sin \left(2 \pi \, \frac{n \, m}{N} \right) \right] \quad \quad &\text{Rectangular Form} \end{aligned}$$ where $x(n)$ is a discrete sequence of time-domain sampled values of the continuous variable $x(t)$, $e$ is th base of natural logarithms and $j=\sqrt{-1}$. The value of $N$ in the DFT equation determines the following things: 1. The number of input samples that are needed 1. The resolution of the frequency-domain results 1. The amount of processing time necessary to compute an N-point DFT

### DFT Concrete Example In this concrete DFT example we will let $N=4$. The rectangular form of this DFT is: $$X(m) = \sum_{n=0}^{3} x(n) \left[cos \left(2 \pi \, \frac{n \, m}{4} \right) \, - \, j \, sin \left(2 \pi \, \frac{n \, m}{4} \right) \right]$$ $$\bigg\downarrow$$ $$\begin{aligned} X(0) = \quad & x(0) \, cos\left(2 \pi \, \frac{0 \cdot 0}{4}\right) \, - \, j \, x(0) \, sin\left(2 \pi \, \frac{0 \cdot 0}{4}\right) \\ + \, & x(1) \, cos\left(2 \pi \, \frac{1 \cdot 0}{4}\right) \, - \, j \, x(1) \, sin\left(2 \pi \, \frac{1 \cdot 0}{4}\right) \\ + \, & x(2) \, cos\left(2 \pi \, \frac{2 \cdot 0}{4}\right) \, - \, j \, x(2) \, sin\left(2 \pi \, \frac{2 \cdot 0}{4}\right) \\ + \, & x(3) \, cos\left(2 \pi \, \frac{3 \cdot 0}{4}\right) \, - \, j \, x(3) \, sin\left(2 \pi \, \frac{3 \cdot 0}{4}\right) \\ \end{aligned}$$ $$\vdots$$ $$\begin{aligned} X(1) = \quad & x(0) \, cos\left(2 \pi \, \frac{0 \cdot 1}{4}\right) \, - \, j \, x(0) \, sin\left(2 \pi \, \frac{0 \cdot 1}{4}\right) \\ + \, & x(1) \, cos\left(2 \pi \, \frac{1 \cdot 1}{4}\right) \, - \, j \, x(1) \, sin\left(2 \pi \, \frac{1 \cdot 1}{4}\right) \\ + \, & x(2) \, cos\left(2 \pi \, \frac{2 \cdot 1}{4}\right) \, - \, j \, x(2) \, sin\left(2 \pi \, \frac{2 \cdot 1}{4}\right) \\ + \, & x(3) \, cos\left(2 \pi \, \frac{3 \cdot 1}{4}\right) \, - \, j \, x(3) \, sin\left(2 \pi \, \frac{3 \cdot 1}{4}\right) \\ \end{aligned}$$ $$\vdots$$ $$\begin{aligned} X(2) = \quad & x(0) \, cos\left(2 \pi \, \frac{0 \cdot 2}{4}\right) \, - \, j \, x(0) \, sin\left(2 \pi \, \frac{0 \cdot 2}{4}\right) \\ + \, & x(1) \, cos\left(2 \pi \, \frac{1 \cdot 2}{4}\right) \, - \, j \, x(1) \, sin\left(2 \pi \, \frac{1 \cdot 2}{4}\right) \\ + \, & x(2) \, cos\left(2 \pi \, \frac{2 \cdot 2}{4}\right) \, - \, j \, x(2) \, sin\left(2 \pi \, \frac{2 \cdot 2}{4}\right) \\ + \, & x(3) \, cos\left(2 \pi \, \frac{3 \cdot 2}{4}\right) \, - \, j \, x(3) \, sin\left(2 \pi \, \frac{3 \cdot 2}{4}\right) \\ \end{aligned}$$ $$\vdots$$ $$\begin{aligned} X(3) = \quad & x(0) \, cos\left(2 \pi \, \frac{0 \cdot 3}{4}\right) \, - \, j \, x(0) \, sin\left(2 \pi \, \frac{0 \cdot 3}{4}\right) \\ + \, & x(1) \, cos\left(2 \pi \, \frac{1 \cdot 3}{4}\right) \, - \, j \, x(1) \, sin\left(2 \pi \, \frac{1 \cdot 3}{4}\right) \\ + \, & x(2) \, cos\left(2 \pi \, \frac{2 \cdot 3}{4}\right) \, - \, j \, x(2) \, sin\left(2 \pi \, \frac{2 \cdot 3}{4}\right) \\ + \, & x(3) \, cos\left(2 \pi \, \frac{3 \cdot 3}{4}\right) \, - \, j \, x(3) \, sin\left(2 \pi \, \frac{3 \cdot 3}{4}\right) \\ \end{aligned}$$ This process is illustrated in more detail in the code below. Feel free to adjust the input to see how the DFT calculation changes.

import numpy as np
import matplotlib.pyplot as plt

# Sample rate and period
Fs = 8e3
Ts = 1/Fs

# Number of samples to use in DFT calculation
N = 8
n = np.arange(0, N)

# Calculate the stop time
t_stop = N * Ts

# Signal parameters
A_1 = 1
f_1 = 1000
theta_1 = 0
A_2 = 0.5
f_2 = 2000
theta_2 = np.pi * 3 / 4

# Discrete signals
x_1 = A_1 * np.sin(2 * np.pi * f_1 * n * Ts + theta_1)
x_2 = A_2 * np.sin(2 * np.pi * f_2 * n * Ts + theta_2)
x_in = x_1 + x_2

# Simulate a continuous signal by setting the sample rate
# sufficiently higher than the discrete sample rate
Fs_cont = 1000 * Fs
Ts_cont = 1/Fs_cont
t_cont = np.arange(0, t_stop, Ts_cont)

# Continuous signals
x_1_cont = A_1 * np.sin(2 * np.pi * f_1 * t_cont + theta_1)
x_2_cont = A_2 * np.sin(2 * np.pi * f_2 * t_cont + theta_2)
x_in_cont = x_1_cont + x_2_cont

# Plot the composite input data as well as each part of the input data
fig_comp, ax_comp = plt.subplots()
ax_comp.grid()
ax_comp.plot(t_cont/Ts, x_in_cont, label="$x_{in}$")
ax_comp.plot(n, x_in, "ro", label="$x_{in}$ sampled")
ax_comp.plot(t_cont/Ts, x_1_cont, ":", c=(0.7, 0, 0, 0.7),
             label=f"${A_1} \cdot sin(2\pi\,{f_1}\,t + {theta_1:.4g})$")
ax_comp.plot(t_cont/Ts, x_2_cont, "--", c=(0.1, 0.5, 0.3, 0.7),
             label=f"${A_2} \cdot sin(2\pi\,{f_2}\,t + {theta_2:.4g})$")

# Label and format the plot
ax_comp.set_ylim((-1.5, 2.5))
ax_comp.set_xlabel("n")
ax_comp.set_ylabel("Amplitude")
ax_comp.set_title(f"Composite Signal\n$x_{{in}} = {A_1} \cdot sin(2\pi\,{f_1}\,t + {theta_1:.4g}) + {A_2} \cdot sin(2\pi\,{f_2}\,t + {theta_2:.4g})$")
ax_comp.legend(loc="upper right")

# Create a plot of each complex sinusoid that is used to 
# calculate the DFT and perform the DFT calculations
Xm_arr = np.empty(N, dtype=complex)
for m in range(N):
    # Generate the complex sinusoids
    cos_m = np.cos(2 * np.pi * n * m / N)
    sin_m = np.sin(2 * np.pi * n * m / N)
    cos_m_cont = np.cos(2 * np.pi * m * Fs / N * t_cont)
    sin_m_cont = np.sin(2 * np.pi * m * Fs / N * t_cont)

# Plot the complex sinusoids along with the composite input singal
    fig, ax = plt.subplots()
    ax.grid()
    ax.plot(t_cont/Ts, x_in_cont, label="$x_{in}$")
    ax.plot(t_cont/Ts, cos_m_cont, ":", c=(0.7, 0, 0, 0.7),
            label=f"$cos(2 \pi \, {round(m*Fs/N, 3):.4g} \,t)$")
    ax.plot(t_cont/Ts, sin_m_cont, "--", c=(0.1, 0.5, 0.3, 0.7),
            label=f"$sin(2 \pi \, {round(m*Fs/N, 3):.4g} \,t)$")
    ax.plot(n, x_in, "ro")
    ax.plot(n, cos_m, ".", c=(0.7, 0, 0, 0.7))
    ax.plot(n, sin_m, ".", c=(0.1, 0.5, 0.3, 0.7))

# Label and format the plot
    ax.set_ylim((-1.5, 2.5))
    ax.set_xlabel("n")
    ax.set_ylabel("Amplitude")
    ax.set_title(f"m = {m}")
    ax.legend(loc="upper right")

plt.show()

# Step through each of the DFT calculations
    for k in range(N):
        real_str = f"({round(x_in[k], 3)} * {round(cos_m[k], 3)})"
        imag_str = f"j ({round(x_in[k], 3)} * {round(sin_m[k], 3)})"
        if k == 0:
            print(f"X({k}) = ", end="")
            print(f"{real_str:18}-{imag_str}")
        else:
            print(f"     + {real_str:18}-{imag_str}")

print("\n    ........................................................")
    for k in range(N):
        real_str = f"{round(x_in[k] * cos_m[k], 3)}"
        imag_str = f"j {round(x_in[k] * sin_m[k], 3)}"
        if k == 0:
            print("     = ", end="")
            print(f"{real_str:8}-{imag_str}")
        else:
            print(f"     + {real_str:8}-{imag_str}")

print("\n    ........................................................")
    Xm_real = round(np.sum(x_in * cos_m), 3)
    Xm_imag = -round(np.sum(x_in * sin_m), 3)
    # Round -0.0 to 0.0 to ensure phase is calculated correctly
    if Xm_real == -0.0: Xm_real = 0.0 
    if Xm_imag == -0.0: Xm_imag = 0.0 
    print(f"     = {Xm_real} - j {Xm_imag}")

print("\n    ........................................................")
    Xm = complex(Xm_real, Xm_imag)
    Xm_arr[m] = Xm
    print(f"     = {round(np.abs(Xm), 3)} ∠ {round(np.angle(Xm, deg=True), 3)}°")

# Plot X(m) magnitude
fig_mag, ax_mag = plt.subplots()
ax_mag.grid()
ax_mag.stem(n, np.abs(Xm_arr), basefmt="C0:")
ax_mag.set_xlabel("m")
ax_mag.set_ylabel("Magnitude")
ax_mag.set_title("Magnitude of X(m)")

# Plot X(m) phase
fig_phase, ax_phase = plt.subplots()
ax_phase.grid()
ax_phase.stem(n, np.angle(Xm_arr, deg=True), basefmt="C0:")
ax_phase.set_xlabel("m")
ax_phase.set_ylabel("Angle (°)")
ax_phase.set_title("Phase of X(m) in degrees")
ax_phase.set_ylim((-180, 180))

# Plot X(m) real part
fig_real, ax_real = plt.subplots()
ax_real.grid()
ax_real.stem(n, np.real(Xm_arr), basefmt="C0:")
ax_real.set_xlabel("m")
ax_real.set_ylabel("Magnitude")
ax_real.set_title("Real Part of X(m)")

# Plot X(m) imaginary part
fig_imag, ax_imag = plt.subplots()
ax_imag.grid()
ax_imag.stem(n, np.imag(Xm_arr), basefmt="C0:")
ax_imag.set_xlabel("m")
ax_imag.set_ylabel("Magnitude")
ax_imag.set_title("Imaginary Part of X(m)")
pass

## DFT Symmetry For a real valued input sequence, the complex DFT outputs will be symmetric about $m=\frac{N}{2}$. The real parts of the DFT output are *even symmetric* and the imaginary parts of the DFT output are *odd symmetric*. This means that the real DFT values on either side of $m=\frac{N}{2}$ have the same magnitude, but the imaginary DFT values on either side of $m=\frac{N}{2}$ have their signs flipped. The DFT magnitude is also *even symmetric* about $m=\frac{N}{2}$ and the phase is *odd symmetric* about $m=\frac{N}{2}$.

## DFT Linearity The DFT has a property called *linearity* that states that the DFT of the sum of two signals is equal to the sum of the transforms of each signal. This allows us to take the DFT of signals that contain multiple different frequency components and we will arrive at the same result as if we took the DFT of each individual signal.

## DFT Magnitudes The DFT output magnitudes are scaled by a factor of $N$, where $N$ is the length of the DFT. For a real input signal, the magnitude of the DFT output, $M_r$, will be equal to: $$M_r = A_0 \, \frac{N}{2}$$ where $A_0$ is the amplitude of the real input signal. For a complex input signal, the magnitude of the DFT output, $M_c$, will be equal to: $$M_c = A_0 \, N$$ where $A_0$ is the magnitude of the complex input signal ($A_0 \, e^{j \, 2 \pi \, f \, n \, T_s}$).

## DFT Frequency Axis The frequency that is associated with each value in $X(m)$ is a function of the sample rate, $F_s$. The formula to calculate the frequency for a given $m$ is: $$f(m) = \frac{m \, F_s}{N}$$

## DFT Shifting Theorem An important property of the DFT is known as the *shifting theorem*. This states that a shift in time of a periodic input sequencne manifests itself as a constant phase shift in the angles associated with the DFT results.

## Inverse DFT The inverse discrete Fourier transform (IDFT) is a process that allows us to convert frequency-domain data back to the time-domain representation. The IDFT is defined by the following equations: $$\begin{aligned} x(n) & = \frac{1}{N} \sum_{m=0}^{N-1} X(m) e^{j \, 2 \pi \frac{m \, n}{N}} \quad \quad &\text{Exponential Form} \\ \\ x(n) & = \frac{1}{N} \sum_{m=0}^{N-1} X(m) \left[cos \left( 2 \pi \, \frac{m \, n}{N} \right) \, + \, j \, sin \left(2 \pi \, \frac{m \, n}{N} \right) \right] \quad \quad &\text{Rectangular Form} \end{aligned}$$ Notice that the two things that changed between the DFT and the IDFT equations are that: 1. The IDFT equation has a $\frac{1}{N}$ scale factor to normalize the gain incurred in the DFT. 1. The IDFT has the sign flipped in the exponent of the exponential form, and on the imaginary component in the rectangular form.