About me

About me


RSS feed

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 
        ((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))))
          (loop for k from 0 below (/ (length x) 2)
                (* (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
     (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))

  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.

blog comments powered by Disqus