[Numpy-discussion] Prime size FFT: bluestein transform vs general chirp/z transform ?

[David Cournapeau]

Hi,

I finally took the time to clean up my code to speed up prime-size FFT
(which use a O(N^2) algo in both numpy and scipy). The code is there:
https://github.com/cournape/numpy/tree/bluestein (most of the code is
tests, because numpy.fft had almost none).

Bottom line: it is used only for prime numbers, and is faster than the
current code for complex transforms > 500. Because of python +
inherent bluestein overhead, this is mostly useful for "long" fft
(where the speed up is significant - already 100x speed up for prime
size ~ 50000).

Several comments:
  - the overhead is pretty significant (on my machine, bluestein
transfrom is slower for prime size < 500)
  - it could be used as such for real transforms, but the overhead
would be even more significant (there is no bluestein transform for
real transforms, so one needs to re-rexpress real transforms in term
of complex ones, multiplying the overhead by 2x). There are several
alternatives to make things faster (Rader-like transform, as used by
fftw), but I think this would be quite hard to do in python without
significant slowdown, because the code cannot be vectorized.
  - one could also decide to provide a chirp-z transform, of which
Bluestein transform is a special case. Maybe this is more adapted to
scipy ?
  - more generic code will require a few simple (but not trivial)
arithmetic-like functions (find prime factors, find generator of Z/nZ
groups with n prime, etc...). Where should I put those ?

cheers,

David