[Numpy-discussion] strange runtimes of numpy fft

Oscar Benjamin oscar.j.benjamin at gmail.com
Mon Nov 18 17:35:13 EST 2013


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()


Oscar



More information about the NumPy-Discussion mailing list