[Numpy-discussion] Numpy FFT.FFT slow with certain samples

[Oscar Benjamin]

On 1 September 2015 at 11:38, Joseph Codadeen <jdmc80 at hotmail.com> wrote:
>
>> And while you zero-pad, you can zero-pad to a sequence that is a power of
>> two, thus preventing awkward factorizations.
>
> Does numpy have an easy way to do this, i.e. for a given number, find the
> next highest number (within a range) that could be factored into small,
> prime numbers as Phil explained? It would help if it gave a list,
> prioritised by number of factors.

Just use the next power of 2. Pure powers of 2 are the most efficient
for FFT algorithms so it potentially works out better than finding a
smaller but similarly composite size to pad to. Finding the next power
of 2 is easy to code and never a bad choice.

To avoid the problems mentioned about zero-padding distorting the FFT
you can use the Bluestein transform as below. This pads up to a power
of two (greater than 2N-1) but does it in a special way that ensures
that it is still calculating the exact (up to floating point error)
same DFT:

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