In the previous post we discussed the discrete Fourier transform (DFT). Its J implementation was pretty straightforward, and the program looked almost identical to its mathematical definition. In this post, I want to explore how to implement fast Fourier transform (FFT), a recursively defined version of DFT that reduces algorithmic complexity from $O(n^2)$ to $O(N\log{N})$.
Math
Recall the DFT formula from the previous post
I use a slightly different notation here to make recursive definition easier to understand.
We assume that N is a power of 2. Using basic arithmetic and properties of complex exponentials, we can derive the formula for the DFT of order 2N in terms of two DFTs of order N
where $m=0,\dots,N1$. The base case, when $N=1$, is
i.e. the identity function.
Programming
Let’s express the FFT algorithm in words. If $N=1$ then the fft
function is just the identity function (]
in J). Overwise, it’s a function involving a recursive call on an input vector of half the length. Here is the same phrase written in J
1


If $N>1$, we split the input vector into two vectors by even and odd indices. We apply FFT recursively to each of these vectors. We multiply the first output vector by 1 and the second one by $e^{\pi im/N}$. Then, we find the sum and difference of the result vectors, and concatenate them. That should be it.
Let’s start with the split function. To split a vector by even and odd indices in J, we can split it by groups of two elements, stack the groups on top of each other making an $N\times2$ matrix, then transpose it to get $2\times N$ matrix, which gives us the two vectors we need. Here is the implementation of the split function
1


and here is an example of applying this function to a range vector
spl 0 1 2 3 4 5 6 7
0 2 4 6
1 3 5 7
The next step is the function $e^{\pi im/N}$. We implement it in two steps: first, scalar $\omega=e^{\pi i/N}$, then vector $\boldsymbol\omega=\omega^{m}$. The implementation is analogous to the one from the previous post, with the difference of $\pi i/N$ instead of $2\pi i/N$, and vector $\boldsymbol\omega$ instead of matrix $\Omega$
1 2 

The new matrix $\Omega$ is build by stacking a unit vector $[1,\dots,1]$ of length N on top of the vector $\boldsymbol\omega$
1


Here is an example of its application to $N=4$
Omg 4
1 1 1 1
1 0.707107j_0.707107 2.22045e_16j_1 _0.707107j_0.707107
Now we can go back to the fft
function and add some code to the something_recursive_here
placeholder
1


$:"1
means the fft
function is recursively applied to every row of its input matrix (which is the result of the spl
function above). Omg@#@{.
means: take the first row of the input matrix, get the length of it, and pass it to the Omg
function.
final_piece_here
is then what we described above as “find the sum and difference of the result vectors, and concatenate them”. In J it’s implemented as the fork +/ , /
. Now we have all the pieces in place, and the final program is
1 2 3 4 5 
