Direct implementation of the discrete Fourier Transform

Q6.7.1Direct implementation of the discrete Fourier Transform

Compare the speed of execution of NumPy's np.fft.fft algorithm and that of the direct implementation of the equation

F_k = \sum_{m=0}^{n-1}f_m\exp\left( -\frac{2\pi i m k}{n} \right), \quad k=0,1,2,\cdots, n-1

Hints: treat the direct equation as a matrix multiplication (dot product) of an array of $$n$$ function values (random ones will do) with the $n \times n$ array with entries $\exp(-2\pi i m k / n)$ ( $m, k = 0, 1, \cdots, n-1$ ). Use IPython's %timeit magic.

Show solution

The two methods for calculating the DFT can be timed using the IPython %timeit magic function:

In [x]: import numpy as np
In [x]: n = 512
In [x]: # Our input function is just random numbers
In [x]: f = np.random.rand(n)

In [x]: # Time the NumPy (Cooley-Tukey) DFT algorithm
In [x]: %timeit np.fft.fft(f)
3.67 μs ± 21 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [x]: # Now calculate the DFT by direct summation
In [x]: k = np.arange(n)
In [x]: m = k.reshape((n, 1))
In [x]: w = np.exp(-2j * np.pi * m * k / n)
In [x]: %timeit np.dot(w, f)
76.4 μs ± 1.04 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [x]: # Check the two methods produce the same result
In [x]: ftfast = np.fft.fft(f)
In [x]: ftslow = np.dot(w, f)
In [x]: np.allclose(ftfast, ftslow)
Out[x]: True

The Cooley-Tukey algorithm implemented by np.fft.fft is found to be almost 30 times faster than the direct method. In fact, this algorithm can be shown to scale as $\mathcal{O}(n\log n)$ compared with $\mathcal{O}(n^2)$ for direct summation.