Fast Fourier Transform

Definition

FFT rapidly computes DFT by factorizing the DFT matrix As a result, It manages to reduce the complexity of computing the DCT from $O(n^2)$ to $O(n\log n)$, where $n$ is a data size.

Algorithm

Consider a problem evaluating $n$ points in an $n-1$ degree polynomial $P(x)=p_0+p_{1}x+p_2{x}^2+\cdots +p_{n-1}{x}^{n-1}$ It costs $n^2$ multiplication operations.

We can split the terms of the polynomial into odd and even parts $P(x)=(p_0+p_2{x}^2+\cdots +p_{n-1}{x}^{n-1})+(p_{1}x+p_3{x}^3+\cdots +p_{n-2}{x}^{n-2})=\underset{P_e}{\underbrace{(p_0+p_2{x}^2+\cdots +p_{n-1}{x}^{n-1})}}+x\underset{P_o}{\underbrace{(p_1+p_3{x}^2+\cdots +p_{n-2}{x}^{n-1})}}$ $P(x)=P_e(x)+x{P}_o(x)$

Both even $P_e$ and odd $P_o$ parts are $\frac{n}{2}-1$ degree polynomials for $x^2$, so these are even functions. Therefore, the following is satisfied $P(-x)=P_e(-x)-x{P}_o(-x)=P_e(x)-x{P}_o(x)$

So, if we take $n$ points as negative-positive pairs $\pm x_1,\pm x_2,\dots ,\pm x_{n/2}$, we can divide the original problem into two small problems.

Now the problem is to evaluate $\frac{n}{2}$ points $\pm x_1^2,\pm x_2^2,\dots ,\pm x_{n/2}^2$ on the $\frac{n}{2}-1$ degree polynomials $P_e,P_o$ So, the cost become $n^2$ to $2({\frac{n}{2}}^2)$ multiplication operations.

It seems to be the subject of a recursive algorithm. (can apply same logic to $P_e,P_o$ each) However, since the evaluating points are squared form don’t have positive-negative pairs, so recursion breaks.

But if we extend the domain of the points to complex number, this become possible.

The roots of an equation $z^n=1$, which lie on the unit circle on a complex plane with equal distance, satisfy the above condition.

So, If we choose the roots of the equation $1,{\omega }^1,{\omega }_2,\dots ,{\omega }^{n-1}$ as evaluating points, the $2n$-Squares of roots also consist of positive-negative pairs, and can exploit recursion.

Finally, the cost become $n$ to $n{\log }_{2}n$

Let’s express the original problem that evaluating the $n$ points ${\omega }^0,{\omega }^1,{\omega }_2,\dots ,{\omega }^{n-1}$ in the polynomial $P(x)=p_0+p_{1}x+p_2{x}^2+\cdots +p_{n-1}{x}^{n-1}$

$\left[\begin{array}{c}P({\omega }^0)\\ P({\omega }^1)\\ P({\omega }^2)\\ ⋮\\ P({\omega }^{n-1})\end{array}\right]=\underset{\text{DFT matrix}}{\underbrace{\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& \omega & {\omega }^2& \cdots & {\omega }^{n-1}\\ 1& {\omega }^2& {\omega }^4& \cdots & {\omega }^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{n-1}& {\omega }^{2(n-1)}& \cdots & {\omega }^{(n-1)(n-1)}\end{array}\right]}}=\left[\begin{array}{c}{p}_0\\ p_1\\ p_2\\ ⋮\\ p_{n-1}\end{array}\right]$ where $\omega =\exp (\frac{2\pi i}{n})$

The matrix is exactly the DCT matrix. In other words, we can accelerate DFT using the FFT algorithm. The inverse of the matrix also has similar structure, so we can apply the same logic to the calculation $\left[\begin{array}{c}{p}_0\\ p_1\\ p_2\\ ⋮\\ p_{n-1}\end{array}\right]=\frac{1}{n}\underset{\text{IDFT matrix}}{\underbrace{\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }^{-1}& {\omega }^{-2}& \cdots & {\omega }^{-(n-1)}\\ 1& {\omega }^{-2}& {\omega }^{-4}& \cdots & {\omega }^{-2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{-(n-1)}& {\omega }^{-2(n-1)}& \cdots & {\omega }^{-(n-1)(n-1)}\end{array}\right]}}=\left[\begin{array}{c}P({\omega }^0)\\ P({\omega }^1)\\ P({\omega }^2)\\ ⋮\\ P({\omega }^{n-1})\end{array}\right]$ where $\omega =\exp (\frac{2\pi i}{n})$