On 14 November 2013 17:19, David Cournapeau <courn...@gmail.com> wrote:
> On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <char...@crunch.io> wrote:
>>
>> Can you post the raw data? It seems like there are just a couple of "bad"
>> sizes, I'd like to know more precisely what these are.
>
> Indeed. Several of the sizes generated by logspace(2, 7, 25) are prime
> numbers, where numpy.fft is actually O(N^2) and not the usual O(NLogN).
>
> There is unfortunately no freely (aka BSD-like licensed) available fft
> implementation that works for prime (or 'close to prime') numbers, and
> implementing one that is precise enough is not trivial (look at Bernstein
> transform for more details).
I was interested by this comment as I wasn't aware of this aspect of
numpy's fft function (or of fft algorithms in general). Having finally
found a spare minute I've implemented the Bluestein algorithm based
only on the Wikipedia page (feel free to use under any license
including BSD).
Is there anything wrong with the following? It's much faster for e.g.
the prime size 215443 (~1s on this laptop; I didn't wait long enough
to find out what numpy.fft.fft would take).
from numpy import array, exp, pi, arange, concatenate
from numpy.fft import fft, ifft
def ceilpow2(N):
'''
>>> ceilpow2(15)
16
>>> ceilpow2(16)
16
'''
p = 1
while p < N:
p *= 2
return p
def fftbs(x):
'''
>>> data = [1, 2, 5, 2, 5, 2, 3]
>>> from numpy.fft import fft
>>> from numpy import allclose
>>> from numpy.random import randn
>>> for n in range(1, 1000):
... data = randn(n)
... assert allclose(fft(data), fftbs(data))
'''
N = len(x)
x = array(x)
n = arange(N)
b = exp((1j*pi*n**2)/N)
a = x * b.conjugate()
M = ceilpow2(N) * 2
A = concatenate((a, [0] * (M - N)))
B = concatenate((b, [0] * (M - 2*N + 1), b[:0:-1]))
C = ifft(fft(A) * fft(B))
c = C[:N]
return b.conjugate() * c
if __name__ == "__main__":
import doctest
doctest.testmod()
Oscar
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
