[Numpy-discussion] strange runtimes of numpy fft
Charles Waldman
charles at crunch.io
Tue Nov 19 11:00:49 EST 2013
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>
More information about the NumPy-Discussion
mailing list