Home ~ About

---

%matplotlib inline

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import time

Radix-2 Fast Fourier Transform

This is an introduction to the famous Fast Fourier Transform algorithm devised by Cooley and Tuckey in the sixties. The goal of the FFT algorithm is to solve the Discrete Fourier Transform (DFT) in $O(nlog(n))$ time complexity, significantly improving on the naive $O(n^2)$ implementation.

DFT recap

The definition of the DFT is

\begin{equation} X_k = \sum^{N-1}_{n=0}x_n \cdot e^{-j\frac{2\pi}{N}kn} \end{equation}

, which can be expanded with the Euler's formula to \begin{equation} X_k = \sum^{N-1}_{n=0}x_n \cdot (\cos{\frac{2\pi}{N}kn} - j \sin{\frac{2\pi}{N}kn}) \end{equation} We can also write the exponential term as

Matrix form of the DFT

\begin{equation} e^{-j\frac{2\pi}{N}kn} = W_N^{nk} \end{equation}

Quite simply, W is a N by N matrix, each row of which is multiplied by the x vector and summed up to get the transformed data point.

Now the matrix operation version of the DFT looks like

\begin{equation} X_k = \sum_{n=0}^{N-1} x_n W_N^{nk} = W_N \cdot x \end{equation}\begin{equation} \begin{bmatrix} X_0 \\ X_1 \\ X_2 \\ X_3 \end{bmatrix}= \begin{bmatrix} W^{00} & W^{01} & W^{02} & W^{03} \\ W^{10} & W^{11} & W^{12} & W^{13} \\ W^{20} & W^{21} & W^{22} & W^{23} \\ W^{30} & W^{31} & W^{32} & W^{33} \\ \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \end{bmatrix} \end{equation}

The key thing to understand here is that the W matrix does not depend on the value vector x, but only the integer indexes n and k. Here we can instantly see an improvement, the nk and kn values must be identical for the W matrix.

We can write the W matrix open for a simple data set of N = 4. Use python for that

N = 4
W = np.zeros((N,N))* (0+0j)

## Naive (n^2) filling of W matrix
for k in range(N):
    for n in range(N):
        W[k,n] = np.exp(np.complex(0, -2*np.pi/N*k*n))

print(W.round(0))

[[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j]
 [ 1.+0.j  0.-1.j -1.-0.j -0.+1.j]
 [ 1.+0.j -1.-0.j  1.+0.j -1.-0.j]
 [ 1.+0.j -0.+1.j -1.-0.j  0.-1.j]]

\begin{equation} W_N^{nk} = e^{-j\frac{2\pi}{N}kn} = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & -j & -1 & j \\ 1 & -1 & 1 & -1 \\ 1 & j & -1 & -j \\ \end{bmatrix} \end{equation}

Here we can confirm that the matrix is symmetric about the diagonal, only depending on the values of k and n. The matrix W can only have N different values, so calculating the required values can easily be optimized to $O(n)$. Note that here we could simply get the DFT by multiplying this W matrix with the input vector x.

Now we know that we can optimize the algorithm by halving the computation w.r.t the W matrix, but it still gives us $O(n^2)$ computation time because we still have to calculate the product of the W matrix and the x vector. So this doesn't really solve anything. But it showed us the redundancy that exists, now we just need to eliminate it.

Divide and conguer

For simplicity, we start by assuming that the data contains $n^2$ values. A divide and conguer algorithm can be used to drastically reduce the amount of required computation. The approach is based on splitting the input data to two groups, even and odd values, both groups having N/2 values.

\begin{equation} X_k = E_k + W^k_N O_k \end{equation}

Where E corresponds to the equation for the even values and O to the odd values.

For even values:

\begin{equation} E_k = \sum_{n=0}^{N/2-1}x_{2n}W_{N}^{2kn} \end{equation}

And for odd values:

\begin{equation} O_{k} = \sum_{n=0}^{N/2-1}x_{2n+1}W_{N}^{(2n+1)k} \end{equation}

The most important detail is in the indexing of the W matrix, as we previously noticed, the indexing of the W matrix only depends on the values of k and n. Let's see what this means in the exponential terms.

The even indexing: \begin{equation} W_{N}^{2kn} = e^{-j\frac{2\pi}{N}2kn} = e^{-j\frac{2}{\pi}2kn}= W_{N/2}^{kn} \end{equation}

For odd indexing: \begin{equation} W_{N}^{(2n+1)k} = e^{-j\frac{2\pi}{N}(2n+1)k} = e^{-j\frac{2\pi}{N/2}kn}e^{-j\frac{2\pi}{N}k} = W_{N/2}^{kn}W_{N}^k \end{equation}

This means that we just effectively reduced the size of the W matrix from N x N to N/2 x N/2! It's clearly seen now that repeating this yields a logarithmic descent of the size.

Now we do this until N/2 = 2, and the divide and conquer is ready. The W matrices are now so small that precomputing them is trivial. We can also further exploit the symmetric nature of the W matrix as because of the complex conjugation the last half of the $W^k_{N}$ vector is equal to negative of the first half. For our splitted DFT this means

First half: \begin{equation} X_k = E_n + W^k_N O_n, k<N/2 \end{equation}

Second half: \begin{equation} X_k = E_n - W^k_N O_n, k>=N/2 \end{equation}

Then we can just concatenate the two $X_n$ vectors.

Note that we assume here that the size N is a power of two (Radix-2 FFT). We are ready to implement the algorithm using recursion.

def naiveFFT(x):
    ''' The naive implementation for comparison '''
    N = x.size
    X = np.ones(N)*(0+0j)

    for k in range(N):
        A = np.ones(N)*(0+0j)
        for n in range(N):
            A[n] = x[n]*np.exp(-np.complex(0, 2*np.pi*k*n/N))
        X[k] = sum(A)

    return X

def FFT(x):
    ''' Recursive radix-2 FFT'''
    x = np.array(x, dtype=float)
    N = int(x.size)
    # Use the naive version when the size is small enough
    if N <= 8:
        return naiveFFT(x)
    else:
        # Calculate first half of the W veco
        k = np.arange(N//2)
        W = np.exp(-2j*np.pi*k/N)
        evens = FFT(x[::2])
        odds = FFT(x[1::2])
        return np.concatenate([evens + (W * odds), evens - (W * odds)])

    return 0

# Test with the data from last notebook
def f_sin(t):
    return np.sin(4*t*2*np.pi) + np.sin(1.5*t*2*np.pi)
# Create a time-series signal
N = 2**10
t = np.linspace(0, 10, N)
T = t[1]-t[0]
signal = f_sin(t)
fft = FFT(signal)

# Calculate the frequency scale for the plot
freq_scale = np.linspace(0,1/T, N)

print("sampling_frequency: " + str(1/T))

# Plot the results
f, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(t, signal)
ax1.set_title('Signal')
ax2.plot(freq_scale, np.absolute(fft))
ax2.set_title('DFT')
plt.plot

sampling_frequency: 102.3

<function matplotlib.pyplot.plot(*args, scalex=True, scaley=True, data=None, **kwargs)>

For fun we can now test our FFT against the built in FFT in the numpy library

# Test it with random digits. Again, this has to be a power of two because we made a radix-2 implementation
x = np.random.rand(256)

t0 = time.time()
naive = naiveFFT(x)
t1 = time.time()
print("Naive: {}".format(t1-t0))

better = FFT(x)
t2 = time.time()
print("Better: {}".format(t2-t1))

t3 = time.time()
numpyfft = np.fft.fft(x)
print("Numpy: {}".format(t3-t2))

# Check that all had same result
if np.allclose(naive, numpyfft) and np.allclose(better, numpyfft):
    print("\n Correct results")
else:
    print("\n Wrong results")

Naive: 0.30092906951904297
Better: 0.016707897186279297
Numpy: 0.0001480579376220703

 Correct results

Further optimization by removing recursion

Success, we managed to create a O(n log(n)) algorithm. Hovever, we are still quite far away from the numpy's version. From my knowledge with computer science I know that recursion sucks. Function calls are slow and with too long inputs we run to memory issues on the callstack. So we should improve the algorithm by removing recursion.

In this version of the algorithm, we start from bottom up. We reshape the input to segments with the length where we would stop splitting the arrays in the recursive version, and calculate the DFT of all of the small segments simultaneously, then calculating is just a matter of combining the subarrays using same method as we did in the recursive version.

def nonrecursiveFFT(x):
    t0 = time.time()

    x = np.array(x)
    N = x.shape[0]
    stop = 16
    x_split = x.reshape((stop, -1)) # Split the input to n columns of size 16

    n = np.arange(stop)
    k = n.reshape((stop, 1))
    W = np.exp(-2j*np.pi*k*n/stop) # W matrix is of shape (stop, x)

    X = np.dot(W, x_split) # This calculates the FFT for each of the columns

    # Each loop should halve the number of columns in the splitted X-matrix
    # We continue unitll there is only one column
    while X.shape[0] < N:
        evens = X[:, :X.shape[1] // 2] # // is integer division in python
        odds = X[:, X.shape[1] // 2:]
        W = np.exp(-1j * np.pi * np.arange(X.shape[0])/ X.shape[0])[:, None]
        X = np.vstack([evens + W * odds, evens - W * odds])


    return np.concatenate(X)

x = np.random.rand(16)
nonrecursive = nonrecursiveFFT(x)
numpyfft = np.fft.fft(x)

if (np.allclose(nonrecursive, numpyfft)):
    print("Nonrecursive implementation gives the correct result")
else:
    print("Incorrect result")

Nonrecursive implementation gives the correct result

x = np.random.rand(2048*32)

t0 = time.time()
nonrecursive = nonrecursiveFFT(x)
t1 = time.time()
print("Nonrecursive: {}".format(t1-t0))

t2 = time.time()
numpyfft = np.fft.fft(x)
print("Numpy: {}".format(t2-t1))

if np.allclose(nonrecursive, numpyfft):
    print("\n Correct results")
else:
    print("\n Wrong results")

Nonrecursive: 0.058702945709228516
Numpy: 0.00020313262939453125

 Correct results

Closing remarks

This is good enough for understanding the algorithm, but it could be optimized a lot by counting the W vectors or "Twiddle factors" is still repeated a lot in the nonrecursive version, by doing that more effectively we could speed up the algorithm. Also with larger input sizes the memory accessing becomes an issue and that can't easily be optimized using python. (Numpy relies on backends implemented in more low-level languages).