Implement bluestein algorithm for prime size FFT (Trac #579)

[numpy-gitbot]

Original ticket http://projects.scipy.org/numpy/ticket/579 on 2007-09-14 by @cournape, assigned to unknown.

See #1104 and http://en.wikipedia.org/wiki/Bluestein's_FFT_algorithm

[numpy-gitbot]

Milestone changed to 1.1 by @cournape on 2007-09-14

[numpy-gitbot]

Milestone changed to 1.3.0 by @cournape on 2008-08-13

[numpy-gitbot]

Milestone changed to Unscheduled by @cournape on 2009-03-02

[numpy-gitbot]

@endolith wrote on 2009-12-01

From #1104:

According to

http://www.fftw.org/fftw-paper-ieee.pdf

FFTW implements Bluestein's algorithm. A python implementation is available at

http://www.mail-archive.com/numpy-discussion@scipy.org/msg01812.html

Also, here is [http://healpix.jpl.nasa.gov/html/libfftpack/bluestein_8c-source.html Bluestein algorithm] written for a C port of fftpack, [http://healpix.jpl.nasa.gov/html/libfftpack/group__fftgroup.html#_details described here]

[endolith]

So as I understand it, the Bluestein algorithm was originally developed to do FFTs, and then the concept was expanded into the Chirp Z-transform, and DFT and zoom-FFT are special cases of CZT (CZT does spirals through Z plane, while circles and arcs are subtypes of spiral), so if we have a CZT then this can do the prime FFTs and zoom FFTs.

Turns out this has already been submitted to SciPy under public domain by Paul Kienzle and Nadav Horesh: [SciPy-user] Chirp Z transform and it works:

x = rand(42073) # prime length

allclose(czt(x), fft(x))
Out[37]: True

timeit czt(x)
10 loops, best of 3: 53.9 ms per loop

timeit fft(x)
1 loops, best of 3: 2.35 s per loop

2.35 / 53.9e-3 # speedup factor
Out[40]: 43.59925788497217

@stefanv said:

Thank you both for the contribution! I've reviewed the code and think
it should be included in SciPy. Would someone else like to have a
look before I commit it?

Where would the best place for this be? scipy.fftpack?

The submission probably does belong in scipy.fftpack because it's specialized and has classes to generate functions, etc. But this bug is for numpy, should some variant be included in numpy too, or should it all be in numpy?

[endolith]

There's also https://github.com/cournape/numpy/compare/bluestein which "is only an implementation optimization to deal with fft size which cannot be factorized in small factors" and doesn't include zoom FFT etc.

So pull request that and put the other solution in scipy?

http://numpy-discussion.10968.n7.nabble.com/Prime-size-FFT-bluestein-transform-vs-general-chirp-z-transform-td3171.html

[endolith]

So I started merging the CZT functions by Paul Kienzle and Nadav Horesh, and modified it to have good accuracy as a prime-length FFT replacement.

N = 100003

is_prime(N)
Out[132]: True

timeit is_prime(N)
1000000 loops, best of 3: 755 ns per loop

a = randn(N)

allclose(czt(a), fft(a))
Out[135]: True

The CZT requires pre-computing large chirp signals, though, which means calling it like this isn't as efficient as it could be. Their code has classes that create a function for a specific-length transform so that it can be called repeatedly on multiple arrays of the same length and the chirps are only generated once. Still, the inefficient function version is faster for prime sizes above ~500:

timeit czt(a)
10 loops, best of 3: 85 ms per loop

timeit fft(a)
1 loops, best of 3: 18.7 s per loop

Should a simplified FFT-only version of this be added to NumPy?

I think in FFTW they always pre-compute which algorithm will be fastest anyway, before applying it to lots of data? So generating the chirps would be part of that.

[mreineck]

Numpy 1.17 contains an FFT implementation supporting Bluestein's algorithm.

[endolith]

@mreineck Which pull request is that from?