Algorithm - Fast Fourier Transform

Subscribe Send me a message

#algorithm #fast fourier transform #convolution

Introduction

In this post, we present a brief introduction to the fast Fourier transform algorithm and how it can be used to calculate convolutions. We will remove the mathematical details and try to only keep the core algorithm part.

Context

Suppose we have $n$ data points $ (y_0, y_1, ..., y_{n-1}) $, the problem we want to solve is to calculate the Discrete Fourier Transform:

$$ c_k = \sum_{j=0}^{n-1} y_i \; e^{-i2\pi{jk/n}} \;\;\; \textrm{where} \; k = 0, 1, ..., n-1 $$

For each coefficient $c_k$, it takes $O(n)$ to compute. As there are $n$ coefficients to calculate, a naive solution will have $O(n^2)$ complexity.

Reformulate the Problem

Let $z = e^{-i2\pi{}/n} $ be the n-th root of unity and $P(x) = y_0 + y_1 x + ... + y_{n-1}x^{n-1}$. It can be shown that

$$ c_k = P(z^k) $$

Now the problem becomes how to evaluate the polynomial $P$ of degree $n-1$ at the points of the n-th roots of unity. More specifically, how can we construct the following sequence efficiently?

$$ P(1), P(z), ..., P(z^{n-1}) $$

We will assume $n = 2s$ is an even number. We can easily extend the reasoning to cover odder number case.

Divide and Conquer

The Fast Fourier Transform is an elegant algorithm. It illustrates the power of the divide-and-conquer design pattern. To apply the divide-and-conquer technique, we need to find recursive structures in the problem and most importantly, we need to reduce the original problem to a subproblem with a smaller size.

If we look at a polynomial $Q(x) = a_0 + a_1x + ... + a_m x^{m}$ with $m$ being an even number. It can be written as

$$ Q(x) = (a_0 + a_2 x^2 + a_4 x^4 + ... + a_{m} x^{m}) + x(a_1 + a_3 x^2 + ... + a_{m-1} x^{m-2}) $$

As we can see, the two parts on the right side of the equation are two polynomials of degree m/2 evaluated at $x^2$. It indicates that there exists a recursive structure in the polynomial representation and calculation.

The Fast Fourier Transform leverages the same idea. Back to our polynomial $P(x)$, it can be shown that there exist two polynomials of degree $s-1$, such that

$$ \begin{eqnarray} P_{\textrm{even}}(z^{2k}) & = & P(z^{2k}) \\ P_{\textrm{odd}}(z^{2k}) & = & P(z^{2k+1}) \end{eqnarray} $$

And we have the following recursive structure:

(Note that the "box" represents a collection and the actual order of the items may be different.)

Let $T(n)$ be the time complexity of the original problem (i.e. evaluating $P(x)$ at $n$-th root of unity ), then we have

$$ T(n) = 2T(n/2) + O(n) $$

The term $O(n)$ comes from the work of results reordering.

We conclude that the complexity of the Fast Fourier Transform is $O(n\log{}n)$.

Technical Details

In the lecture note Fast Fourier Transform, the author explains how we can re-order the results. To re-order the results, we need to know the mapping between the index and the exponent:

The author uses the example of index = 6. The binary representation of 6 is 110 and the exponent at index 6 is $(011)_2 = 3$. Note that there is a symmetric structure here: the exponent at index 3 is 6 as highlighted in the figure above.

So how can we establish this relationship?

Let $\alpha$ denote the odd/even subscript sequence and $z_{\alpha}$ be the root of unity for that level. Moreover, we assume

$$ P_{\alpha}(z_{\alpha}^k) = A_{\alpha}k + B_{\alpha} $$

The objective is to figure out the relationship between $B_{\alpha}$ and $B_{\alpha{},o/e}$. We know that

$$ \begin{eqnarray} P_{\alpha{},e}(z_{\alpha{}}^{2k}) & = & P_{\alpha{}}(z_{\alpha{}}^{2k}) \\ P_{\alpha{},o}(z_{\alpha{}}^{2k}) & = & P_{\alpha{}}(z_{\alpha{}}^{2k+1}) \end{eqnarray} $$

We also know that

$$ \begin{eqnarray} z_{\alpha{},e} & = & z_{\alpha}^2 \\ z_{\alpha{},o} & = & z_{\alpha}^2 \end{eqnarray} $$

It follows that

$$ \begin{eqnarray} P_{\alpha, e}(z_{\alpha{},e} ^ k) = P_{\alpha}(z_{\alpha} ^ {2k}) = P(z^{A_{\alpha}(2k) + B_{\alpha}}) \\ P_{\alpha, o}(z_{\alpha{},e} ^ k) = P_{\alpha}(z_{\alpha} ^ {2k+1}) = P(z^{A_{\alpha}(2k+1) + B_{\alpha}}) \\ \end{eqnarray} $$

Hence,

$$ \begin{eqnarray} B_{\alpha{}, e} & = & B_{\alpha{}} \\ B_{\alpha{}, o} & = & B_{\alpha{}} + A_{\alpha} \end{eqnarray} $$

The above recursive formula tells us that the late $o$ in the sequence has more contribution to the exponent value. The figure below shows how we construct the exponent value from the index.