[SciPy-Dev] Re FFTs in SciPy (and NumPy)

[Sturla Molden]

There has been some discussion on FFTPACK lately. Problems with FFTPACK 
seems to be:

- Written in old Fortran 77.
- Unprecise for single precision.
- Can sometimes be very slow, depending on input size.
- Can only handle a few small prime factors {2,3,4,5} efficiently.
- How to control integer size portably for 64-bit and f2py?
- Only 1D transforms.

There is another FFT library in cwplib from Colorado Mining School. It 
has some very nice properties:

- Written in C, not Fortran, and the FFT-part of the library is just one 
single C file.
- License is BSD it seems.
- Quite fast (see benchmarks in http://www.fftw.org/fftw-paper.pdf and 
http://www.fftw.org/speed/)
- Can handle larger set of prime factors than FFTPACK 
{2,3,4,5,7,8,9,11,13,16}.
- 1D and 2D inplace transforms.

The problem with cwplib is that it does not allow arbitrary sized FFTs. 
But til will be fast whenever FFTPACK is fast (it seems FFTPACK is only 
competitive for N=2**k case). It will also be fast for other 
factorizations where FFTPACK will fallback to O(N**2).

How to handle arbitrary sized arrays with cwplib? One possibility is 
using Bluestein's FFT algorithm, which just requires a few lines og 
Python. This was even suggested to avoid O(N**2) fallback in FFTPACK 
here a while ago. If we are going to do that for FFTPACK, the cwplib's 
FFT will actually be more flexible with respect to array size.

All in all I think it is worth looking at.

The C code is in the "Seismic U*nx" download from CWP, which is huge, 
but contains a lot of nice C library functions for scientific computing 
apart from FFTs as well (all of which are BSD licensed it seems). 
http://www.cwp.mines.edu/cwpcodes/index.html

Regards,
Sturla