[Numpy-discussion] strange runtimes of numpy fft

[Charles Waldman]

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.harris at gmail.com
> wrote:

>
>
>
> On Mon, Nov 18, 2013 at 3:35 PM, Oscar Benjamin <
> oscar.j.benjamin at gmail.com> wrote:
>
>> On 14 November 2013 17:19, David Cournapeau <cournape at gmail.com> wrote:
>> > On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <charles at 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 at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20131119/603ac322/attachment.html>