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