[Numpy-discussion] strange runtimes of numpy fft

[Charles R Harris]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20131118/60e013df/attachment.html>