Fast Fourier Transform
31st May 2018
I've recently been looking for good floating-point benchmarks, to test a forthcoming version of uLisp  which will include floating-point support, and one I selected was a Fast Fourier Transform.
My starting point was this routine on Rosetta Code :
(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))
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 .
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.
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