An important concept for understanding the Fourier transform is the inner product, which is defined for two complex vectors $x, y \in \mathbb{C}^N$ of length $N \in \mathbb{N}$ as follows:
$$ \langle x | y \rangle := \sum_{n=0}^{N-1} x(n) \overline{y(n)}.$$The inner product may be interpreted as a measure of similarity between $x$ and $y$. The function np.vdot computes the dot product, but note that the complex conjugate is performed on the first argument! So, for computing $\langle x | y \rangle$ as defined above, one has to call np.vdot(y, x).
import numpy as np
x = np.array([1, 1j, 1 + 1j])
y = np.array([1.1 , 1j, 0.9 + 1.1j])
print('Vectors of high similarity:', np.vdot(y, x))
x = np.array([1, 1j, 1 + 1j])
y = np.array([1.1 , -1j, 0.1])
print('Vectors of low similarity:', np.vdot(y, x))
The discrete Fourier transform (DFT) of a signal $x(n)$ of length $N \in \mathbb{N}$ is given by
$$ X(k) := \sum_{n=0}^{N-1} x(n) \exp(-2 \pi i k n / N) $$for $k \in [0:N-1]$. By setting
$$\mathbf{u}_k(n) := \exp(2 \pi i k n / N)$$we can also interpret the DFT as an inner product:
$$ X(k) := \sum_{n=0}^{N-1} x(n) \overline{\mathbf{u}_k} = \langle x | \mathbf{u}_k \rangle.$$A noisy sinusoid shall serve as example signal to illustrate how this works:
%matplotlib inline
from matplotlib import pyplot as plt
N = 1024
n = np.arange(N)
k = 3
x = np.cos(2 * np.pi * k * n / N) + 0.2 * (np.random.rand(N) - 0.5)
plt.figure(figsize=(6, 6))
plt.plot(n, x, 'k')
plt.xlabel('$n$')
plt.ylabel('$x(n)$');
for this_k in [k-1, k, k+1]:
u = np.exp(2j * np.pi * this_k * n / N)
u_bar = np.conj(u)
dot_product = np.vdot(u, x)
plt.figure(figsize=(6*3, 6))
plt.suptitle('$k = %d$' % this_k)
plt.subplot(1,3,1)
plt.plot(n, x, 'k', label='$x(n)$')
plt.plot(n, np.real(u_bar), 'r', label='$\mathrm{Re}(\overline{\mathbf{u}_k(n)})$')
plt.fill_between(n, 0, x * np.real(u_bar), label='Area($x(n) \mathrm{Re}(\overline{\mathbf{u}_k(n)})$)', color='lightgray')
plt.legend(loc='lower left')
plt.xlabel('$n$')
plt.ylabel('$x(n)$')
plt.subplot(1,3,2)
plt.plot(n, x, 'k', label='$x(n)$')
plt.plot(n, np.imag(u_bar), 'r', label='$\mathrm{Im}(\overline{\mathbf{u}_k(n)})$')
plt.fill_between(n, 0, x * np.imag(u_bar), label='Area($x(n) \mathrm{Im}(\overline{\mathbf{u}_k(n)})$)', color='lightgray')
plt.legend(loc='lower left')
plt.xlabel('$n$')
plt.ylabel('$x(n)$')
plt.subplot(1,3,3)
plt.title(r'Bar Plot for $\langle x | \mathbf{u}_k \rangle$')
plt.bar([0, 1], [np.real(dot_product), np.imag(dot_product)], color='lightgray')
plt.xticks([0, 1], [r'$\mathrm{Re}(\langle x | \mathbf{u}_k \rangle)$', r'$\mathrm{Im}(\langle x | \mathbf{u}_k \rangle)$'])
plt.ylim([-600, 600])
plt.grid()
The $\mathrm{DFT}_N \in \mathbb{C}^{N\times N}$ matrix is given by
$$\mathrm{DFT}_N(n, k) = \mathrm{exp}(-2 \pi i k n / N) = \mathrm{exp}(-2 \pi i / N)^{kn}$$for $n\in[0:N-1]$ and $k\in[0:N-1]$. Defining $\Omega_N := \mathrm{exp}(-2 \pi i / N)$, this matrix looks like
$$ \mathrm{DFT}_N = \begin{pmatrix} 1 & 1 & 1 & \dots & 1 \\ 1 & \Omega_N & \Omega_N^2 & \dots & \Omega_N^{N-1} \\ 1 & \Omega_N^2 & \Omega_N^4 & \dots & \Omega_N^{2(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \Omega_N^{N-1} & \Omega_N^{2(N-1)} & \dots & \Omega_N^{(N-1)(N-1)} \\ \end{pmatrix} $$def generate_matrix_dft(N, K):
pass # your implementation here!
N = 32
dft_mat = generate_matrix_dft(N, N)
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.title('$\mathrm{Re}(\mathrm{DFT}_N)$')
plt.imshow(np.real(dft_mat), origin='lower', cmap='seismic', aspect='equal')
plt.colorbar()
plt.subplot(1,2,2)
plt.title('$\mathrm{Im}(\mathrm{DFT}_N)$')
plt.imshow(np.imag(dft_mat), origin='lower', cmap='seismic', aspect='equal')
plt.colorbar()
plt.tight_layout()
def dft(x):
pass # your implementation here!
N = 64
n = np.arange(N)
x = np.sin(2 * np.pi * 5 * n / N )
X = dft(x)
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.title('$x$')
plt.plot(x, 'k')
plt.subplot(1,2,2)
plt.title('$|X|$')
plt.plot(np.abs(X), 'k')
plt.tight_layout()
Now we want to approach the FFT as described in the book Chapter 2.4.3.
def twiddle(N):
pass # your implementation here!
One can split a DFT of size $N$ into two DFTs of size $N/2$. For this purpose, one has to compute the DFT of the even-indexed and the uneven-indexed entries of $x$:
\begin{align} (A(0), \dots, A(N/2-1)) &= \mathrm{DFT}_{N/2} \cdot (x(0), x(2), x(4), \dots, x(N-2))\\ (B(0), \dots, B(N/2-1)) &= \mathrm{DFT}_{N/2} \cdot (x(1), x(3), x(5), \dots, x(N-1)) \end{align}With those two DFTs of size $N/2$, one can compute the full DFT of size $N$:
\begin{align} X(k) &= A(k) + \Omega_N^k \cdot B(k)\\ X(N/2 + k) &= A(k) - \Omega_N^k \cdot B(k)\\ \end{align}for $k \in [0: N/2 - 1]$.
def dft_two_times_half_size(x):
pass # your implementation here!
N = 64
n = np.arange(N)
x = np.sin(2 * np.pi * 5 * n / N )
X = dft_two_times_half_size(x)
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.title('$x$')
plt.plot(x, 'k')
plt.subplot(1,2,2)
plt.title('$|X|$')
plt.plot(np.abs(X), 'k')
plt.tight_layout()
We already showed how to compute a DFT of size $N$ by splitting it into two DFTs of size $N/2$. By Recursively further splitting all DFTs, until we have only DFTs of length $1$, we reach the fast Fourier Transform (FFT).
def fft(x):
pass # your implementation here!
N = 64
n = np.arange(N)
x = np.sin(2 * np.pi * 5 * n / N )
X = fft(x)
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.title('$x$')
plt.plot(x, 'k')
plt.subplot(1,2,2)
plt.title('$|X|$')
plt.plot(np.abs(X), 'k')
plt.tight_layout()
We gained a lot!
N = 512
n = np.arange(N)
x = np.sin(2 * np.pi * 5 * n / N )
print('Timing for DFT: ', end='')
%timeit dft(x)
print('Timing for FFT: ', end='')
%timeit fft(x)
import timeit
Ns = [2 ** n for n in range(1,11)]
times_dft = []
times_fft = []
execuctions = 3
for N in Ns:
n = np.arange(N)
x = np.sin(2 * np.pi * 5 * n / N )
time_dft = timeit.timeit(lambda: dft(x), number=execuctions) / execuctions
time_fft = timeit.timeit(lambda: fft(x), number=execuctions) / execuctions
times_dft.append(time_dft)
times_fft.append(time_fft)
plt.figure(figsize=(16,5))
plt.plot(Ns, times_dft, '-xk', label='DFT')
plt.plot(Ns, times_fft, '-xr', label='FFT')
plt.xticks(Ns)
plt.legend()
plt.grid()
plt.xlabel('$N$')
plt.ylabel('Runtime (s)');