uLisp

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.

28th December 2018: I've improved the program; it's now more compact and runs faster.

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 all assume the arguments are complex numbers:

(defun c+ (x y)
  (cx (+ (car x) (car y)) (+ (cdr x) (cdr y))))

(defun c- (x y)
  (cx (- (car x) (car y)) (- (cdr x) (cdr y))))

(defun c* (x y)
  (cx
   (- (* (car x) (car y)) (* (cdr x) (cdr y)))
   (+ (* (cdr x) (car y)) (* (car x) (cdr y)))))

(defun c/ (x y)
  (let ((den (+ (* (car y) (car y)) (* (cdr y) (cdr y)))))
    (cx
     (/ (+ (* (car x) (car y)) (* (cdr x) (cdr y))) den)
     (/ (- (* (cdr x) (car y)) (* (car x) (cdr y))) den))))

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

(defun cexp (x)
  (cx
   (* (exp (car x)) (cos (cdr x)))
   (* (exp (car x)) (sin (cdr x)))))

(defun cabs (x)
  (sqrt (+ (* (car x) (car x)) (* (cdr x) (cdr x)))))

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) (cx (* pi (incf k)) 0)) (cx (length x) 0))) j)) 
               odd)))
      (append (mapcar c+ even aux) (mapcar c- even aux)))))

For convenience, here's a version of fft that takes a list of real numbers and converts them to complex numbers before calling fft:

(defun fftr (r) (fft (mapcar (lambda (x) (cx x 0)) r)))

Testing it out:

> (pprint (fftr '(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 32-point FFT on the different uLisp platforms, as a measure of floating-point performance. The test used was:

(time (fftr (let (z) (dotimes (x 16) (push 0 z)) (dotimes (x 16) (push 1 z)) z)))

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

For a summary of the performance on different 32-bit platforms see Performance.

1. ^ Fast Fourier transform on Rosetta Code.