The Fourier analysis in general and the Fourier transform in particular has always been one of my favourite topics in mathematics. When I started learning the J programming language I thought that implementing the discrete Fourier transform (DFT) in J must be really elegant. One rainy weekend I finally found time to do that, and this post is a result of my experiments.

Math

The definition of one-dimensional DFT is pretty simple. It is just a function that maps a given vector f of length N to a vector F of the same length by the following formula

To implement this formula the only thing we need is basic arithmetic for complex numbers. Not only does J provide such an arithmetic, it has special code for complex exponentiation, which makes the implementation much easier than it would be in other languages.

Before we start coding let’s rewrite the formula in terms of matrix multiplication so that we can utilize this feature of J. Let $\omega$ be the primitive Nth root of unity $\omega = e^{2\pi i/N}$, then the DFT formula becomes

where f is an input vector and $\Omega$ is a symmetric matrix $(\omega^{-mn})$ with $m,n=0,\dots,N-1$.

Programming

We are now ready to code. Let’s start with $\omega = e^{2\pi i/N}$. As I said before J has special code for function $e^{\pi x}$, namely the composition ^@o.@. In our case $x$ is equal to $2i/N$, where $N$ is an input to the function $\omega$. In other words, $x$ is a function that divides $2i$ by its input, which in J is written as 0j2&%. That means the complete little omega implementation is

omg=: ^@o.@(0j2&%)

The next step is the matrix $\Omega=(\omega^{-mn})$. We can think of it this way: make a multiplication table of vectors $[0,\dots,N-1]$ and $-[0,\dots,N-1]$, then raise $\omega$ to the power of every element from that table. In J the multiplication table of a vector by its negative is implemented as the hook */ -. For building a range-vector $[0,\dots,N-1]$ we can use the function i.. Therefore, the implementation of a square matrix $(-mn)$ in J is the composition (*/ -)@i., and the function of $\Omega$ is

Omg=: omg ^ (*/ -)@i.

The last piece of the DFT program is the formula $\mathbf{f}\cdot\Omega$. Our Omg function accepts as a parameter the length of the input vector, while the dft function does the vector itself. That means we need to compose Omg function with the length function #. Then the composition Omg@# will work on vectors.

The matrix (and vector) multiplication function in J is +/ .*. Therefore, the product of the input vector and the function Omg@# is just a hook

dft=: +/ .* Omg@#

That’s it. This is a complete implementation of DFT in J. Let’s put all three lines together.

omg=: ^@o.@(0j2&%)

Omg=: omg ^ (*/ -)@i.

dft=: +/ .* Omg@#

The whole program fits in one tweet. If you like one-liners, you can do this as well

dft=: +/ .* ^@o.@(0j2&% * (*/ -)@i.)@#

References

  • Source code with examples on GitHub