## Fast Fourier Transform

31st May 2018

I've recently been looking for good floating-point benchmarks, to test a forthcoming version of uLisp [1] which will include floating-point support, and one I selected was a Fast Fourier Transform.

```(defun fft (x)
(if (<= (length x) 1) x
(let*
((even (fft (loop for i from 0 below (length x) by 2 collect (nth i x))))
(odd (fft (loop for i from 1 below (length x) by 2 collect (nth i x))))
(aux
(loop for k from 0 below (/ (length x) 2)
collect
(* (exp (/ (* (complex 0 -2) pi k ) (length x)))
(nth k odd)))))
(append (mapcar #'+ even aux) (mapcar #'- even aux)))))```

The formatting of the original suggests it was written by a C programmer rather than a Lisp programmer, and I've tidied it up a bit above.

The test case is this simple waveform:

`(fft '(1 1 1 1 0 0 0 0))`

which gives:

```CL-USER > (pprint (fft '(1 1 1 1 0 0 0 0)))

(#C(4.0D0 -0.0D0)
#C(1.0D0 -2.414213562373095D0)
#C(0.0D0 0.0D0)
#C(1.0D0 -0.4142135623730949D0)
#C(0.0D0 0.0D0)
#C(0.9999999999999999D0 0.4142135623730949D0)
#C(0.0D0 0.0D0)
#C(0.9999999999999997D0 2.414213562373095D0))```

### A more elegant solution

However, in searching for solutions to this problem I found an elegant version in Scheme using a purely functional style of programming with no variable assignments, by Prasenjit Saha [3].

Inspired by his version here's a functional Common Lisp version that avoids variable assignments, or the use of the loop macro, and to my mind is much easier to understand:

```(defun evens (f)
(if (null f) nil
(cons (car f) (evens (cddr f)))))

(defun odds (f)
(if (null f) nil
(cons (cadr f) (odds (cddr f)))))

(defun rotate (fun k l lis)
(if (null lis) nil
(cons
(funcall fun k l (car lis))
(rotate fun (1+ k) l (cdr lis)))))

(defun ph (k l j)
(* (exp (/ (* (complex 0 -2) (* pi k)) l)) j))

(defun plusminus (a b)
(append (mapcar #'+ a b) (mapcar #'- a b)))

(defun fft2 (x)
(if (= (length x) 1) x
(plusminus (fft2 (evens x))
(rotate #'ph 0 (length x) (fft2 (odds x))))))```

The routines evens and odds return lists containing every alternate element:

```CL-USER > (evens '(0 1 2 3 4 5 6 7))
(0 2 4 6)

CL-USER > (odds '(0 1 2 3 4 5 6 7))
(1 3 5 7)```

The function rotate is a bit more complicated. It takes four parameters: The first is a function, the second and third are integers, and the fourth is a list. It applies the function to the first two parameters and each element of the list, with the first number incremented in successive calls. For example:

```CL-USER > (rotate #'list 0 'x '(a b c d e))
((0 X A) (1 X B) (2 X C) (3 X D) (4 X E))```

The function ph performs the central calculation of the FFT, and plusminus takes two lists of equal length, and returns a single list of their pairwise sums followed by their pairwise differences:

```CL-USER > (plusminus '(1 3 5) '(0 2 4))
(1 5 9 1 1 1)```

Finally fft2 recursively applies the Cooley–Tukey algorithm to calculate the FFT.

Result:

```CL-USER > (pprint (fft2 '(1 1 1 1 0 0 0 0)))

(#C(4.0D0 -0.0D0)
#C(1.0D0 -2.414213562373095D0)
#C(0.0D0 0.0D0)
#C(1.0D0 -0.4142135623730949D0)
#C(0.0D0 0.0D0)
#C(0.9999999999999999D0 0.4142135623730949D0)
#C(0.0D0 0.0D0)
#C(0.9999999999999997D0 2.414213562373095D0))```

1. ^ uLisp - Lisp for the Arduino, Micro Bit, and MSP430.
2. ^ Fast Fourier transform - Common Lisp on Rosetta Code.
3. ^ A Fast Fourier Transform in Lisp on physik.uzh.ch.