## 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.

My starting point was this routine on Rosetta Code ^{[2]}:

(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))

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

blog comments powered by Disqus