CL-FFT October 7th, 2009
Patrick Stein

This is a Common Lisp library to do a Fast Fourier Transform on a multi-dimensional array of numbers. The array can be any number of dimensions, but each dimension must be a power-of-two in size.

Acquisition

Use

Here is a simple example of its use on a two-dimensional input array. This example does a Fourier transform on the data, then it does an inverse Fourier transform on the results.

(require :asdf)
(asdf:operate 'asdf:load-op 'fft)

(defparameter *buf* #2A((1 2 3 4)(5 6 7 8)))

(let ((transformed (fft:fft *buf*)))
    (fft:ifft transformed))

The (fft:fft …) and (fft:ifft …) functions take an optional second argument which specifies the destination buffer. The destination buffer must be the same dimensions as the source buffer. Further, the destination buffer must have an :element-type of (complex double-float).

If you are using a Lisp implementation that supports threads, you can use a parallelized version of the FFT with the #:pfft package. Note: The parallelization happens over rows and columns in multi-dimensional Fourier transforms. There is no advantage here to using the parallel version on a one-dimensional array since the whole row will be placed in the same job. For multi-dimensional arrays, however, the speed-up can be significant: processing the (geometrically) parallel rows and columns in parallel (temporally).

(require :asdf)
(asdf:operate 'asdf:load-op 'pfft)

(defparameter *buf* #2A((1 2 3 4)(5 6 7 8)))

(let ((transformed (pfft:pfft *buf*)))
    (pfft:pifft transformed))

Like their single-threaded kin, (pfft:pfft …) and (pfft:pifft …) functions take an optional second argument which specifies the destination buffer. The destination buffer must be the same dimensions as the source buffer. Further, the destination buffer must have an :element-type of (complex double-float).

The parallel version uses PCall which in turn uses Bordeaux threads. If Bordeaux threads is unable to provide threads for your Lisp implementation, the (pfft:pfft …) and (pfft:pifft …) functions revert to their non-parallelized kin. As such, you don’t need to worry about when to use which. You can use the parallel all of the time, if you like.

Performance

Here is a performance breakdown at the moment for my Macbook Pro 2.4GHz Intel Core 2 Duo with 4GB of RAM. All times are in seconds. Note: there is a penalty in the parallel version for having too many small jobs.

  ‘(1048576) ‘(512 512) ‘(1024 1024) ‘(256 256 64)

SBCL 1.0.37 (64-bit)
(single-thread)
0.44 0.07 0.37 2.45
SBCL 1.0.37 (64-bit)
(multi-thread)
0.45 0.07 0.25 4.19
SBCL 1.0.30

0.73 0.13 0.66 3.61
CMUCL 20A Unicode

0.67 0.11 0.58 3.57
Clozure64 v1.3
(single-thread)

28.07 2.39 13.01 129.08
Clozure64 v1.3
(multi-thread)

30.00 1.80 10.39 165.58
CLISP 2.47

36.12 8.45 39.32 215.29

I still have a great deal of work to do to keep things from consing so terribly on platforms other than CMUCL and SBCL.

Comparison to Bordeaux-FFT

After I wrote this, several people pointed out Bordeaux-FFT. I missed it because it wasn’t on the Cliki at the time and doesn’t show up very well on Google when you search for Fourier transform and lisp. It does, however, show up well when you search for fft and lisp. Oops.

Anyhow, Bordeaux-FFT only operates on one-dimensional arrays. This implementations operates on n-dimensional arrays. However, Bordeaux-FFT is somewhat faster than this implementation. Bordeaux-FFT runs about 25% faster than this implementation on the 1,048,576-sample array under SBCL and CLISP (0.55 seconds instead of 0.71 under SBCL). This code is slightly faster than Bordeaux-FFT under Clozure.

Bordeaux-FFT also provides additional functions for windowing a signal that are not provided in this library.

Caveats

Some implementations do not scale the results of the forward transform at all and scale the results of the inverse transform down by the number of samples in the array. I have chosen, instead, to scale both the forward and inverse transforms down by the square root of the number of samples in the array. This keeps them in closer proximity to each other.

Also, I have centered the results so that the lowest frequencies are in the center of the array instead of at the origin. This is the typical way that Fourier data is presented. But, not all libraries do this.