# Understanding the FFT Algorithm

https://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/

The Fast Fourier Transform (FFT) is one of the most important algorithms in signal processing and data analysis. I've used it for years, but having no formal computer science background, It occurred to me this week that I've never thought to ask how the FFT computes the discrete Fourier transform so quickly. I dusted off an old algorithms book and looked into it, and enjoyed reading about the deceptively simple computational trick that JW Cooley and John Tukey outlined in their classic 1965 paper introducing the subject.

The goal of this post is to dive into the Cooley-Tukey FFT algorithm, explaining the symmetries that lead to it, and to show some straightforward Python implementations putting the theory into practice. My hope is that this exploration will give data scientists like myself a more complete picture of what's going on in the background of the algorithms we use.

For simplicity, we'll concern ourself only with the forward transform, as the inverse transform can be implemented in a very similar manner. Taking a look at the DFT expression above, we see that it is nothing more than a straightforward linear operation: a matrix-vector multiplication of$\vec{x}$



Forward Discrete Fourier Transform (DFT):

$$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N}$$

Inverse Discrete Fourier Transform (IDFT):

$$x_n = \frac{1}{N}\sum_{k=0}^{N-1} X_k e^{i~2\pi~k~n~/~N}$$


```python
import numpy as np
def DFT_slow(x):
    """Compute the discrete Fourier Transform of the 1D array x"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)
#We can double-check the result by comparing to numpy's built-in FFT function:
x = np.random.random(1024)
np.allclose(DFT_slow(x), np.fft.fft(x))
```




    True



Just to confirm the sluggishness of our algorithm, we can compare the execution times of these two approaches:


```python
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
```

    111 ms ± 19.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
    3.89 µs ± 59.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)


We are over 1000 times slower, which is to be expected for such a simplistic implementation. But that's not the worst of it. For an input vector of length N, the FFT algorithm scales as$ O[NlogN]$, while our slow algorithm scales as $O[N2]$. That means that for $N=10^6$ elements, we'd expect the FFT to complete in somewhere around 50 ms, while our slow algorithm would take nearly 20 hours!

So how does the FFT accomplish this speedup? The answer lies in exploiting symmetry.

# Symmetries in the Discrete Fourier Transform

One of the most important tools in the belt of an algorithm-builder is to exploit symmetries of a problem. If you can show analytically that one piece of a problem is simply related to another, you can compute the subresult only once and save that computational cost. Cooley and Tukey used exactly this approach in deriving the FFT.

We'll start by asking what the value of $X_{N+k}$ is. From our above expression:

$$\begin{align*}
X_{N + k} &=  \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~(N + k)~n~/~N}\\
          &= \sum_{n=0}^{N-1} x_n \cdot e^{- i~2\pi~n} \cdot e^{-i~2\pi~k~n~/~N}\\
          &= \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N}
\end{align*}$$

where we've used the identity $exp[2π i n]=1$ which holds for any integer n.

The last line shows a nice symmetry property of the DFT:

$$X_{N+k}=X_k$$

By a simple extension,

$$X_{k + i \cdot N} = X_k$$

for any integer i. As we'll see below, this symmetry can be exploited to compute the DFT much more quickly.

# DFT to FFT: Exploiting Symmetry

Cooley and Tukey showed that it's possible to divide the DFT computation into two smaller parts. From the definition of the DFT we have:

$$\begin{align}
X_k &= \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N} \\
    &= \sum_{m=0}^{N/2 - 1} x_{2m} \cdot e^{-i~2\pi~k~(2m)~/~N} + \sum_{m=0}^{N/2 - 1} x_{2m + 1} \cdot e^{-i~2\pi~k~(2m + 1)~/~N} \\
    &= \sum_{m=0}^{N/2 - 1} x_{2m} \cdot e^{-i~2\pi~k~m~/~(N/2)} + e^{-i~2\pi~k~/~N} \sum_{m=0}^{N/2 - 1} x_{2m + 1} \cdot e^{-i~2\pi~k~m~/~(N/2)}
\end{align}$$

We've split the single Discrete Fourier transform into two terms which themselves look very similar to smaller Discrete Fourier Transforms, one on the odd-numbered values, and one on the even-numbered values. So far, however, we haven't saved any computational cycles. Each term consists of $(N/2)∗N $computations, for a total of $N^2$.

The trick comes in making use of symmetries in each of these terms. Because the range of $k$ is $0≤k<N$, while the range of$ n$ is $0≤n<M≡N/2$, we see from the symmetry properties above that we need only perform half the computations for each sub-problem. Our$ O[N^2]$ computation has become $O[M^2]$, with M half the size of $N$.

But there's no reason to stop there: as long as our smaller Fourier transforms have an even-valued M, we can reapply this divide-and-conquer approach, halving the computational cost each time, until our arrays are small enough that the strategy is no longer beneficial. In the asymptotic limit, this recursive approach scales as $O[NlogN]$.

This recursive algorithm can be implemented very quickly in Python, falling-back on our slow DFT code when the size of the sub-problem becomes suitably small:


```python
import numpy as np
def FFT(x):
    """A recursive implementation of the 1D Cooley-Tukey FFT"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    
    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    elif N <= 32:  # this cutoff should be optimized
        return DFT_slow(x)
    else:
        X_even = FFT(x[::2])
        X_odd = FFT(x[1::2])
        factor = np.exp(-2j * np.pi * np.arange(N) / N)
        return np.concatenate([X_even + factor[:N / 2] * X_odd,
                               X_even + factor[N / 2:] * X_odd])                        
```


```python
x = np.random.random(1024)
np.allclose(FFT(x), np.fft.fft(x))
```


```python

```

6d489446b1c4f127d69c0291491078a5ef663e4e
ipynb
Jupyter Notebook
code/fft/Understanding the FFT.ipynb
xing710/ModSimPy
87f0f481926c40855223e2843bd728edb235c516
MIT
code/fft/Understanding the FFT.ipynb
xing710/ModSimPy
87f0f481926c40855223e2843bd728edb235c516
MIT
code/fft/Understanding the FFT.ipynb
xing710/ModSimPy
87f0f481926c40855223e2843bd728edb235c516
MIT
Qwen/Qwen-72B
1. YES
2. YES