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
