next up previous index
Next: A Recursive Implementation of Up: Spectral Methods Previous: Discrete Fourier Transform

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

\begin{displaymath}\omega_n^k = e^{-2\pi i k/n}
\end{displaymath}

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:

 \begin{displaymath}
{}^8\boldsymbol{\Lambda} =
\left(
\begin{array}{llllllll}
...
...8^4& \omega_8^3& \omega_8^2& \omega_8^1\\
\end{array}\right)
\end{displaymath} (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 fk in the following way:

\begin{displaymath}[1,2,3,4,5,6,7,8]\longrightarrow [1,3,5,7,2,4,6,8],
\end{displaymath}

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

 \begin{displaymath}
\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)
\end{displaymath} (2.11)

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

and that

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

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

\begin{displaymath}\left(
\begin{array}{llll}
1 & 0 & 0 & 0 \\
0 & \omega_8 ...
...\Lambda} = {}^8\boldsymbol{D}_4 \cdot {}^4\boldsymbol{\Lambda}
\end{displaymath} (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

\begin{displaymath}\left(
\begin{array}{llll}
\omega_8^4 & \omega_8^4& \omega_...
...3 & \omega_8^1& \omega_8^7& \omega_8^5\\
\end{array} \right)
\end{displaymath} (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:

 
8F(f[1:8]) = $\displaystyle \left(
\begin{array}{cc}
{}^4\boldsymbol{\Lambda}, & {}^8\boldsym...
...y}{ll}
\boldsymbol{f}_{[1:8:2]} \\
\boldsymbol{f}_{[2:8:2]}
\end{array}\right)$  
  = $\displaystyle \left(
\begin{array}{cc}
\boldsymbol{1} & {}^8\boldsymbol{D}_4 \\...
... \\
{}^4\boldsymbol{\Lambda} \cdot \boldsymbol{f}_{[2:8:2]}
\end{array}\right)$  
  = $\displaystyle \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 4F(f[1:8:2]) and to 4F(f[2:8:2]) we can write without much further ado:

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

And we repeat it once more to get:

 
2F(f[1:8:4]) = $\displaystyle \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)$  
2F(f[3:8:4]) = $\displaystyle \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)$  
2F(f[2:8:4]) = $\displaystyle \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)$  
2F(f[4:8:4]) = $\displaystyle \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)$  

where I have made use of the fact that 1F(fk) = fk.

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 2kDk 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 2kDk, are the ones that are really needed. Periodicity of $\omega_n^k$is automatically taken care of.



 
next up previous index
Next: A Recursive Implementation of Up: Spectral Methods Previous: Discrete Fourier Transform
Zdzislaw Meglicki
2001-02-26