8000 Implement bluestein algorithm for prime size FFT (Trac #579) · Issue #1177 · numpy/numpy · GitHub
[go: up one dir, main page]
Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Closed
numpy-gitbot opened this issue Oct 19, 2012 · 13 comments
Closed

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

numpy-gitbot opened this issue Oct 19, 2012 · 13 comments
Labels
17 - Task component: numpy.fft

Comments

@numpy-gitbot
Copy link

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
@endolith endolith mentioned this issue Oct 19, 2012
Open
@numpy-gitbot
Copy link
Author

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

@numpy-gitbot
Copy link
Author

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

@numpy-gitbot
Copy link
Author

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

@numpy-gitbot
Copy link
Author

@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]

This was referenced Oct 23, 2012
Open
Closed
@endolith endolith mentioned this issue Oct 20, 2014
Closed
@endolith
Copy link
Contributor

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 8000 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 endolith mentioned this issue Dec 18, 2014
Closed
@endolith
Copy link
Contributor

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
Copy link
Contributor

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
Copy link
Contributor
mreineck commented Apr 5, 2019

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

@endolith
Copy link
Contributor
endolith commented Apr 5, 2019

@mreineck Which pull request is that from?

@mreineck
Copy link
Contributor
mreineck commented Apr 5, 2019

It's #11888

@seberg
Copy link
Member
seberg commented Apr 5, 2019

Cool, one more closed issue :), thanks @mreineck.

@seberg seberg closed this as completed Apr 5, 2019
@fasiha
Copy link
89CF fasiha commented Nov 29, 2019

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

@mreineck, is Bluestein's algorithm available for all functions in numpy.fft, including rfft* irfft* (1D and up, direct and inverse, real FFT)? Or just some subset? I didn't see any mention of "Bluestein" in the 1.17 docs.

@mreineck
Copy link
Contributor

Yes, Bluestein's algorithm is used in all sorts of transforms if they have a length containing too large prime factors. Unfortunately, if Bluestein is used, a real-valued transform has no speed benefit compared to a complex-valued one any more, since the FFTs carried out within Bluestein's algorithm are intrinsically complex-valued, but at least the O(N**2) scaling is avoided.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
17 - Task component: numpy.fft
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
6 participants
@fasiha @endolith @seberg @mattip @mreineck @numpy-gitbot
0