Fast Fourier Transform

The following uLisp program calculates a Fast Fourier Transform (FFT) for a list of numbers. It is a standard benchmark for floating-point performance, and you can see other versions on Rosetta Code [1].

This program will only run on the 32-bit versions of uLisp with floating point.

Complex numbers

First we need to write routines to perform arithmetic with complex numbers. A complex number is represented as a cons of the real and imaginary parts, and the function cx is provided as a convenient way of constructing a complex number:

(defun cx (x y) (cons x y))

Next we provide complex versions of the operators +, -, *, /. These are designed to work with real numbers too:

(defun c+ (x y)
  (let ((a (if (consp x) (car x) x))
        (b (if (consp x) (cdr x) 0))
        (c (if (consp y) (car y) y))
        (d (if (consp y) (cdr y) 0)))
  (cx (+ a c) (+ b d))))

(defun c- (x y)
    (let ((a (if (consp x) (car x) x))
          (b (if (consp x) (cdr x) 0)) 
          (c (if (consp y) (car y) y))
          (d (if (consp y) (cdr y) 0)))
      (cx (- a c) (- b d))))

(defun c* (x y)
  (let ((a (if (consp x) (car x) x))
        (b (if (consp x) (cdr x) 0)) 
        (c (if (consp y) (car y) y))
        (d (if (consp y) (cdr y) 0)))
    (cx
     (- (* a c) (* b d))
     (+ (* b c) (* a d)))))

(defun c/ (x y)
  (let ((a (if (consp x) (car x) x))
        (b (if (consp x) (cdr x) 0)) 
        (c (if (consp y) (car y) y))
        (d (if (consp y) (cdr y) 0)))
    (cx
     (/ (+ (* a c) (* b d)) (+ (* c c) (* d d)))
     (/ (- (* b c) (* a d)) (+ (* c c) (* d d))))))

For FFT we also need exp, and for convenience we also provide abs:

(defun cexp (x)
  (let ((a (if (consp x) (car x) x))
        (b (if (consp x) (cdr x) 0)))
    (cx
     (* (exp a) (cos b))
     (* (exp a) (sin b)))))

(defun cabs (x)
  (let ((a (if (consp x) (car x) x))
        (b (if (consp x) (cdr x) 0)))
    (sqrt (+ (* a a) (* b b)))))

Fast Fourier Transform

Finally, here is the definition of fft:

(defvar pi 3.14159)

(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 fft (x)
  (if (= (length x) 1) x 
    (let*
        ((even (fft (evens x)))
         (odd (fft (odds x)))
         (k -1)
         (aux (mapcar 
               (lambda (j) 
                 (c* (cexp (c/ (c* (cx 0 -2) (c* pi (incf k))) (length x))) j)) 
               odd)))
      (append (mapcar c+ even aux) (mapcar c- even aux)))))

Testing it out:

> (pprint (fft '(1 1 1 1 0 0 0 0)))

((4.0 . 0)
  (1.0 . -2.41421)
  (0 . 0)
  (1.0 . -0.414212)
  (0 . 0)
  (0.999999 . 0.414215)
  (0 . 0)
  (0.999994 . 2.41421))

Here's the whole FFT program: Fast Fourier Transform program.

Benchmark

Here are the results for running a 64-point FFT on the different uLisp platforms, as a measure of floating-point performance. The test used was:

(for-millis ()
  (fft (let (z) (dotimes (x 32) (push 0 z)) (dotimes (x 32) (push 1 z)) z)))

which returns the time in milliseconds for calculating the FFT of a sequence of 32 '1's followed by 32 '0's.

The following table gives a summary of the different versions:

Platform Processor Clock Cells Benchmark
Arduino Due SAM3X8E 84 MHz 10240 0.6 secs
Arduino Zero SAMD21 48 MHz 3072 1.0 secs
Arduino MKRZero SAMD21 48 MHz 3072 1.3 secs
BBC Micro Bit nRF51822 16 MHz 1024 No room

  1. ^ Fast Fourier transform on Rosetta Code.