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.

CL-FFT v1.4.2011.03.24 released March 24th, 2011
Patrick Stein

Elliott Johnson provided me with some patches for my CL-FFT library so that it will work with Allegro modern mode (mlisp). Here is the new tarball: fft_1.4.2011.03.24.tar.gz and its GPG signature: fft_1.4.2011.03.24.tar.gz.asc Thank you!

Common Lisp vs. C/ASM: Battle FFT October 16th, 2009
Patrick Stein

Dr. David McClain of Refined Audiometrics sent me his code for interfacing with the Apple vDSP library. Here it is along with a Makefile that I put together for it: vdsp-wrapper.tar.gz. I have not actually run the code myself. I started to convert his Lispworks FLI code to CFFI code, but bailed on it in […]

Adventures in Package Dependencies October 14th, 2009
Patrick Stein

A variety of people have asked how the speed of my Lisp Fourier Transform Library compares to an FFI call to the FFTW Library. I downloaded the FFTW Library, built it, and installed it. Then, I tried to get someone else’s FFI wrappers for the FFTW Library going. I worked on this for two hours […]

Lisp Fourier Transform Library Faster Yet (v1.3) October 13th, 2009
Patrick Stein

I released version 1.3 of my Common Lisp Fourier Transform library today. It is significantly faster than yesterday’s version. On my MacBook Pro in SBCL 1.0.30 with a 512 by 512 by 16 array, version 1.3 takes 3.74 seconds where version 1.2 took 9.77 seconds and version 1.0 took 25.31 seconds. For a breakdown of […]

Lisp Fourier Transform Library Faster October 12th, 2009
Patrick Stein

I released version 1.2 of my Common Lisp Fourier Transform library today. It is significantly faster than the old version. On my MacBook Pro in SBCL 1.0.30 with a 512 by 512 by 16 array, version 1.2 takes 9.77 seconds where version 1.0 took 25.31 seconds. Version 1.2 consed 2.4M while version 1.0 consed 7.4M. […]

Updates In Email

Email:

l