How about FFTW? I think there are wrappers out there for that ...
On Mon, Nov 18, 2013 at 9:26 PM, Charles R Harris charlesr.har...@gmail.com
wrote:
On Mon, Nov 18, 2013 at 3:35 PM, Oscar Benjamin
oscar.j.benja...@gmail.com wrote:
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()
Where this starts to get tricky is when N is a product of primes not
natively supported in fftpack. The fftpack supports primes 2, 3, 5, 7(?) at
the moment, one would need to do initial transforms to break it down into a
number of smaller transforms whose size would have prime factors supported
by fftpack. Then use fftpack on each of those. Or the other way round,
depending on taste.
Chuck
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
___
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
