Fast Fourier Transform

Fast Fourier Transform is a Discrete Fourier Transform calculated more efficiently. In order to evaluate the Discrete Fourier Transform Vandermonde matrix the way we have done it in the previous section we needed to evaluate somewhat less than $n\times n$ terms because the first row and the first column are all ones. But certain economies can be attempted. This derives from the observation that

$$\omega_n^k = e^{-2\pi i k/n}$$

is periodic in k and the period is n. This means that we don't really have to evaluate all powers between $\omega^{1\times 1}$ and $\omega^{(n-1)\times(n-1)}$, because they will keep repeating. For example the Fourier transform matrix for n=8 would have the following terms:

$${}^8\boldsymbol{\Lambda} = \left( \begin{array}{llllllll} ... ...8^4& \omega_8^3& \omega_8^2& \omega_8^1\\ \end{array}\right)$$ (2.10)

It transpires that in order to evaluate a discrete Fourier Transform all that's needed is to evaluate $\omega_n^1 \ldots \omega_n^{n-1}$. Then you just need to insert those terms cleverly in appropriate places in the Vandermonde matrix - also, observe that matrix (2.10) is symmetric. Now the insertion of $\omega_n^1 \ldots \omega_n^{n-1}$ into appropriate spots is not going to be free: those places have to be computed one way or the other and data may have to be moved around too. But the resulting algorithm is not going to be ${\cal O}(n^2)$. It is going to be ${\cal O}(n\log n)$ instead.

Consider matrix given by (2.10). Let us rearrange the columns of the matrix and the entries of the corresponding vector f_k in the following way:

$$[1,2,3,4,5,6,7,8]\longrightarrow [1,3,5,7,2,4,6,8],$$

i.e., we put all the odd columns first and then the even ones. So now the matrix looks as follows:

$$\left( \begin{array}{l} {}^8F_0 \\ {}^8F_1 \\ {}^8F_2 ... ... f(x_1)\\ f(x_3)\\ f(x_5)\\ f(x_7) \end{array}\right)$$ (2.11)

Now, observe that

 

$$\left( \omega_8^1,\omega_8^3,\omega_8^5,\omega_8^7 \right) \omega_8^1 \cdot \left(1,\omega_8^2,\omega_8^4,\omega_8^6\right)$$ =

$$\left( \omega_8^2,\omega_8^6,\omega_8^2,\omega_8^6 \right) \omega_8^2 \cdot \left(1,\omega_8^4,1 ,\omega_8^4\right)$$ =

$$\left( \omega_8^3,\omega_8^1,\omega_8^7,\omega_8^5 \right) \omega_8^3 \cdot \left(1,\omega_8^6,\omega_8^4,\omega_8^2\right)$$ = (2.12)

and that

$${}^4\boldsymbol{\Lambda} = \left( \begin{array}{llll} 1 & 1... ... 1 & \omega_4^3& \omega_4^2& \omega_4^1 \end{array} \right)$$ (2.13)

is a Vandermonde matrix for the 4-point discrete Fourier transform, and that multiplications in equation (2.12) can be expressed as:

$$\left( \begin{array}{llll} 1 & 0 & 0 & 0 \\ 0 & \omega_8 ... ...\Lambda} = {}^8\boldsymbol{D}_4 \cdot {}^4\boldsymbol{\Lambda}$$ (2.14)

Last, observe that the lower right corner of ${}^8\boldsymbol{\Lambda}_{[1,3,5,7,2,4,6,8]}$is simply $-\,{}^8\boldsymbol{D}_4 \cdot {}^4\boldsymbol{\Lambda}$ because

$$\left( \begin{array}{llll} \omega_8^4 & \omega_8^4& \omega_... ...3 & \omega_8^1& \omega_8^7& \omega_8^5\\ \end{array} \right)$$ (2.15)

and $\omega_8^4 = e^{-2\pi i 4/8} = e^{-\pi i} = -1$.

In summary equation (2.11) can be rewritten as follows:

 

$$\left( \begin{array}{cc} {}^4\boldsymbol{\Lambda}, & {}^8\boldsym... ...y}{ll} \boldsymbol{f}_{[1:8:2]} \\ \boldsymbol{f}_{[2:8:2]} \end{array}\right)$$ ⁸F(f_[1:8]) =

$$\left( \begin{array}{cc} \boldsymbol{1} & {}^8\boldsymbol{D}_4 \\... ... \\ {}^4\boldsymbol{\Lambda} \cdot \boldsymbol{f}_{[2:8:2]} \end{array}\right)$$ =

$$\left( \begin{array}{cc} \boldsymbol{1} & {}^8\boldsymbol{D}_4 \\... ..._{[1:8:2]}) \\ {}^4\boldsymbol{F}(\boldsymbol{f}_{[2:8:2]}) \end{array}\right)$$ = (2.16)

Applying the same reasoning to ⁴F(f_[1:8:2]) and to ⁴F(f_[2:8:2]) we can write without much further ado:

 

$$\left( \begin{array}{cc} \boldsymbol{1} & {}^4\boldsymbol{D}_2 \\... ..._{[1:8:4]}) \\ {}^2\boldsymbol{F}(\boldsymbol{f}_{[3:8:4]}) \end{array}\right)$$ ⁴F(f_[1:8:2]) =

$$\left( \begin{array}{cc} \boldsymbol{1} & {}^4\boldsymbol{D}_2 \\... ..._{[2:8:4]}) \\ {}^2\boldsymbol{F}(\boldsymbol{f}_{[4:8:4]}) \end{array}\right)$$ ⁴F(f_[2:8:2]) = (2.17)

And we repeat it once more to get:

 

$$\left( \begin{array}{cc} \boldsymbol{1} & {}^2\boldsymbol{D}_1 \\... ...ray}\right) \cdot \left( \begin{array}{ll} f(x_0) \\ f(x_4) \end{array}\right)$$ ²F(f_[1:8:4]) =

$$\left( \begin{array}{cc} \boldsymbol{1} & {}^2\boldsymbol{D}_1 \\... ...ray}\right) \cdot \left( \begin{array}{ll} f(x_2) \\ f(x_6) \end{array}\right)$$ ²F(f_[3:8:4]) =

$$\left( \begin{array}{cc} \boldsymbol{1} & {}^2\boldsymbol{D}_1 \\... ...ray}\right) \cdot \left( \begin{array}{ll} f(x_1) \\ f(x_5) \end{array}\right)$$ ²F(f_[2:8:4]) =

$$\left( \begin{array}{cc} \boldsymbol{1} & {}^2\boldsymbol{D}_1 \\... ...ray}\right) \cdot \left( \begin{array}{ll} f(x_3) \\ f(x_7) \end{array}\right)$$ ²F(f_[4:8:4]) =

where I have made use of the fact that ¹F(f_k) = f_k.

And this is Fast Fourier Transform.

When it comes to coding this algorithm, you can either start from the top and then write a recursive procedure, or you can start at the bottom and evaluate four 2-point transforms first, then two 4-point transforms and finally the single 8-point transform.

Observe that you never really evaluate and store the full Vandermonde discrete Fourier tranform matrix. Matrices ^2kD_k are diagonal, so they can be stored simply as vectors and the matrix operation that converts two n-point transforms into one 2n-point transform can be implemented as two simple equations. The only $\omega_n^k$ terms that are evaluated while calculating ^2kD_k, are the ones that are really needed. Periodicity of $\omega_n^k$is automatically taken care of.