strange runtimes of numpy fft

classic Classic list List threaded Threaded
21 messages Options
12
Reply | Threaded
Open this post in threaded view
|

strange runtimes of numpy fft

Max Linke
Hi

I noticed some strange scaling behavior of the fft runtime. For most
array-sizes the fft will be calculated in a couple of seconds, even for
very large ones. But there are some array sizes in between where it will
take about ~20 min (e.g. 400000). This is really odd for me because an
array with 10 million entries is transformed in ~2s. Is this typical for
numpy?

I attached a plot and an ipynb to reproduce and illustrate it.

best Max

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion

test_fft.ipynb (75K) Download Attachment
fft.png (20K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Charles Waldman
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.

It's typical for FFT to perform better at a sample size that is a power of 2, and algorithms like FFTW take advantage of factoring the size, and "sizes that are products of small factors are transformed most efficiently."

  - Charles

On Thu, Nov 14, 2013 at 10:18 AM, Max Linke <[hidden email]> wrote:
Hi

I noticed some strange scaling behavior of the fft runtime. For most
array-sizes the fft will be calculated in a couple of seconds, even for
very large ones. But there are some array sizes in between where it will
take about ~20 min (e.g. 400000). This is really odd for me because an
array with 10 million entries is transformed in ~2s. Is this typical for
numpy?

I attached a plot and an ipynb to reproduce and illustrate it.

best Max

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion



_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

David Cournapeau



On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <[hidden email]> 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).

David

It's typical for FFT to perform better at a sample size that is a power of 2, and algorithms like FFTW take advantage of factoring the size, and "sizes that are products of small factors are transformed most efficiently."

  - Charles

On Thu, Nov 14, 2013 at 10:18 AM, Max Linke <[hidden email]> wrote:
Hi

I noticed some strange scaling behavior of the fft runtime. For most
array-sizes the fft will be calculated in a couple of seconds, even for
very large ones. But there are some array sizes in between where it will
take about ~20 min (e.g. 400000). This is really odd for me because an
array with 10 million entries is transformed in ~2s. Is this typical for
numpy?

I attached a plot and an ipynb to reproduce and illustrate it.

best Max

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion



_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion



_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Robert Kern-2
In reply to this post by Charles Waldman
On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <[hidden email]> 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.
>
> It's typical for FFT to perform better at a sample size that is a power of 2, and algorithms like FFTW take advantage of factoring the size, and "sizes that are products of small factors are transformed most efficiently."

These are the sizes, as given in the notebook he supplied:

array([     100,      161,      261,      421,      681,     1100,
           1778,     2872,     4641,     7498,    12115,    19573,
          31622,    51089,    82540,   133352,   215443,   348070,
         562341,   908517,  1467799,  2371373,  3831186,  6189658, 10000000])

--
Robert Kern

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Max Linke
In reply to this post by Charles Waldman
You can check everything in the notebook, it will generate all the data.
I checked the runtime for sizes in logspace(2, 7, 25). I know that the
fft will work faster for array sizes that are a power of 2 but
differences on 3 orders of magnitude is huge. I also attached a script
generated from the notebook that you can run as well

best Max

On Thu, 2013-11-14 at 10:45 -0600, Charles Waldman 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.
>
> It's typical for FFT to perform better at a sample size that is a power of
> 2, and algorithms like FFTW take advantage of factoring the size, and
> "sizes that are products of small factors are transformed most efficiently."
>
>   - Charles
>
> On Thu, Nov 14, 2013 at 10:18 AM, Max Linke <[hidden email]> wrote:
>
> > Hi
> >
> > I noticed some strange scaling behavior of the fft runtime. For most
> > array-sizes the fft will be calculated in a couple of seconds, even for
> > very large ones. But there are some array sizes in between where it will
> > take about ~20 min (e.g. 400000). This is really odd for me because an
> > array with 10 million entries is transformed in ~2s. Is this typical for
> > numpy?
> >
> > I attached a plot and an ipynb to reproduce and illustrate it.
> >
> > best Max
> >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > [hidden email]
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >
> >
> _______________________________________________
> NumPy-Discussion mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion

test_fft.py (1K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Robert Kern-2
On Thu, Nov 14, 2013 at 5:30 PM, Max Linke <[hidden email]> wrote:
>
> You can check everything in the notebook, it will generate all the data.
> I checked the runtime for sizes in logspace(2, 7, 25). I know that the
> fft will work faster for array sizes that are a power of 2 but
> differences on 3 orders of magnitude is huge.

Primes are especially bad. 215443, for example, is prime. Going from O(N*logN) to O(N**2) where N=215443 accounts for 3 orders of magnitude, at least.

--
Robert Kern

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

David Cournapeau
In reply to this post by Max Linke



On Thu, Nov 14, 2013 at 5:30 PM, Max Linke <[hidden email]> wrote:
You can check everything in the notebook, it will generate all the data.
I checked the runtime for sizes in logspace(2, 7, 25). I know that the
fft will work faster for array sizes that are a power of 2 but
differences on 3 orders of magnitude is huge.

This is expected if you go from N log N to N ** 2 for large N :)

You can for example compare np.fft.fft(a) for 2**16 and 2**16+1 (and 2**16-1 that while bad is not prime, so only 1 order of magnitude slower).

David 

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Jaime Fernández del Río
On Thu, Nov 14, 2013 at 9:37 AM, David Cournapeau <[hidden email]> wrote:

You can for example compare np.fft.fft(a) for 2**16 and 2**16+1 (and 2**16-1 that while bad is not prime, so only 1 order of magnitude slower).

I actually did...

Each step of a FFT basically splits a DFT of size N = P*Q into P DFTs of size Q followed by Q DFTs of size P. If prime length DFTs take time proportional to their length squared, the time complexity can be written as N times the sum of the prime factors of N, which is N log N (actually 2 N log N) for sizes a power of two.

So for 2^16 you actually have 2^16 * 2 * 16 = 2.1e6
For 2^16-1 = 3*5*17*257 you get (2^16-1)*(3+5+17+257) = 18.5e6 (8x slower)
For 2^16+1, which is prime, you get (2^16+1)^2 = 4295.0e6 (2045x slower)

On my system I get:
%timeit np.fft.fft(a) # a = np.ones(2**16)
100 loops, best of 3: 3.18 ms per loop 

%timeit np.fft.fft(b) # b = np.ones(2**16 - 1)
100 loops, best of 3: 13.6 ms per loop (4x slower)

%timeit np.fft.fft(c) # c = np.ones(2**16 + 1)
1 loops, best of 3: 25.1 s per loop (8000x slower)

There clearly are some constant factors missing in this analysis, although it gives reasonable order of magnitude predictions. Doing some more timings it seems that:
 * prime sized inputs perform at least 2x worse than these formulas predict.
 * sizes with factors different than 2, perform about 2x better than these formulas predict.

Jaime

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

David Cournapeau



On Thu, Nov 14, 2013 at 7:05 PM, Jaime Fernández del Río <[hidden email]> wrote:
On Thu, Nov 14, 2013 at 9:37 AM, David Cournapeau <[hidden email]> wrote:

You can for example compare np.fft.fft(a) for 2**16 and 2**16+1 (and 2**16-1 that while bad is not prime, so only 1 order of magnitude slower).

I actually did...

Each step of a FFT basically splits a DFT of size N = P*Q into P DFTs of size Q followed by Q DFTs of size P. If prime length DFTs take time proportional to their length squared, the time complexity can be written as N times the sum of the prime factors of N, which is N log N (actually 2 N log N) for sizes a power of two.

It does not really make sense to talk about 2 N Log N, since the whole point of O notation is to ignore the constant factors.

If you want to care about the # operations, the actual complexity of fft is not n log n, but C n log n with C > 1, and C depends quite a bit on the implementation (the FFTW guys published some details of their implementations). I believe it also depends on how well you can factor the number (i.e. the algos in FFTW are all O(n log n), but the constant factor if you count # operations is actually a function of the input).

David

So for 2^16 you actually have 2^16 * 2 * 16 = 2.1e6
For 2^16-1 = 3*5*17*257 you get (2^16-1)*(3+5+17+257) = 18.5e6 (8x slower)
For 2^16+1, which is prime, you get (2^16+1)^2 = 4295.0e6 (2045x slower)

On my system I get:
%timeit np.fft.fft(a) # a = np.ones(2**16)
100 loops, best of 3: 3.18 ms per loop 

%timeit np.fft.fft(b) # b = np.ones(2**16 - 1)
100 loops, best of 3: 13.6 ms per loop (4x slower)

%timeit np.fft.fft(c) # c = np.ones(2**16 + 1)
1 loops, best of 3: 25.1 s per loop (8000x slower)

There clearly are some constant factors missing in this analysis, although it gives reasonable order of magnitude predictions. Doing some more timings it seems that:
 * prime sized inputs perform at least 2x worse than these formulas predict.
 * sizes with factors different than 2, perform about 2x better than these formulas predict.

Jaime

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion



_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Max Linke
In reply to this post by David Cournapeau
On Thu, 2013-11-14 at 17:19 +0000, David Cournapeau wrote:
> On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <[hidden email]> 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).
Ok I thought the fft is always O(N log N). But it makes sense that this
only works if the input-size factorizes well. Thanks for clearing this
up.

best Max

>
> 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).
>
> David
>
> >
> > It's typical for FFT to perform better at a sample size that is a power of
> > 2, and algorithms like FFTW take advantage of factoring the size, and
> > "sizes that are products of small factors are transformed most efficiently."
> >
> >   - Charles
> >
> > On Thu, Nov 14, 2013 at 10:18 AM, Max Linke <[hidden email]> wrote:
> >
> >> Hi
> >>
> >> I noticed some strange scaling behavior of the fft runtime. For most
> >> array-sizes the fft will be calculated in a couple of seconds, even for
> >> very large ones. But there are some array sizes in between where it will
> >> take about ~20 min (e.g. 400000). This is really odd for me because an
> >> array with 10 million entries is transformed in ~2s. Is this typical for
> >> numpy?
> >>
> >> I attached a plot and an ipynb to reproduce and illustrate it.
> >>
> >> best Max
> >>
> >> _______________________________________________
> >> NumPy-Discussion mailing list
> >> [hidden email]
> >> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >>
> >>
> >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > [hidden email]
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >
> >
> _______________________________________________
> NumPy-Discussion mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/numpy-discussion


_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

David Cournapeau



On Fri, Nov 15, 2013 at 10:47 PM, Max Linke <[hidden email]> wrote:
On Thu, 2013-11-14 at 17:19 +0000, David Cournapeau wrote:
> On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <[hidden email]> 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).
Ok I thought the fft is always O(N log N). But it makes sense that this
only works if the input-size factorizes well. Thanks for clearing this
up.

To be exact, there *are* FFT-like algorithms for prime number size, but that's not implemented in numpy or scipy (see here: http://numpy-discussion.10968.n7.nabble.com/Prime-size-FFT-bluestein-transform-vs-general-chirp-z-transform-td3171.html for an old discussion on that topic).

cheers,

David

 

best Max
>
> 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).
>
> David
>
> >
> > It's typical for FFT to perform better at a sample size that is a power of
> > 2, and algorithms like FFTW take advantage of factoring the size, and
> > "sizes that are products of small factors are transformed most efficiently."
> >
> >   - Charles
> >
> > On Thu, Nov 14, 2013 at 10:18 AM, Max Linke <[hidden email]> wrote:
> >
> >> Hi
> >>
> >> I noticed some strange scaling behavior of the fft runtime. For most
> >> array-sizes the fft will be calculated in a couple of seconds, even for
> >> very large ones. But there are some array sizes in between where it will
> >> take about ~20 min (e.g. 400000). This is really odd for me because an
> >> array with 10 million entries is transformed in ~2s. Is this typical for
> >> numpy?
> >>
> >> I attached a plot and an ipynb to reproduce and illustrate it.
> >>
> >> best Max
> >>
> >> _______________________________________________
> >> NumPy-Discussion mailing list
> >> [hidden email]
> >> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >>
> >>
> >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > [hidden email]
> > http://mail.scipy.org/mailman/listinfo/numpy-discussion
> >
> >
> _______________________________________________
> NumPy-Discussion mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/numpy-discussion


_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion


_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Oscar Benjamin
In reply to this post by David Cournapeau
On 14 November 2013 17:19, David Cournapeau <[hidden email]> wrote:

> On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <[hidden email]> 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
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Charles R Harris



On Mon, Nov 18, 2013 at 3:35 PM, Oscar Benjamin <[hidden email]> wrote:
On 14 November 2013 17:19, David Cournapeau <[hidden email]> wrote:
> On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <[hidden email]> 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
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[hidden email]> wrote:



On Mon, Nov 18, 2013 at 3:35 PM, Oscar Benjamin <[hidden email]> wrote:
On 14 November 2013 17:19, David Cournapeau <[hidden email]> wrote:
> On Thu, Nov 14, 2013 at 4:45 PM, Charles Waldman <[hidden email]> 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
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion



_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Henry Gomersall-2
On 19/11/13 16:00, Charles Waldman wrote:
> How about FFTW?  I think there are wrappers out there for that ...

Yes there are! (complete with the numpy.fft API)
https://github.com/hgomersall/pyFFTW

However, FFTW is dual licensed GPL/commercial and so the wrappers are
also GPL by necessity.

Cheers,

Henry
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Stéfan van der Walt
On Tue, Nov 19, 2013 at 6:03 PM, Henry Gomersall <[hidden email]> wrote:
> However, FFTW is dual licensed GPL/commercial and so the wrappers are
> also GPL by necessity.

I'm not sure if that is true, strictly speaking--you may license your
wrapper code under any license you wish.  It's just that it becomes
confusing when the combined work then has to be released under GPL.

Stéfan
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Henry Gomersall-2
On 19/11/13 16:08, Stéfan van der Walt wrote:
> On Tue, Nov 19, 2013 at 6:03 PM, Henry Gomersall<[hidden email]>  wrote:
>> >However, FFTW is dual licensed GPL/commercial and so the wrappers are
>> >also GPL by necessity.
> I'm not sure if that is true, strictly speaking--you may license your
> wrapper code under any license you wish.  It's just that it becomes
> confusing when the combined work then has to be released under GPL.

This is on shaky GPL ground. I'm inclined to agree with you, but given
any usage necessarily has to link against the FFTW libs, any resultant
redistribution is going to be GPL by necessity. I.e. pyFFTW is useless
without FFTW so it's just simpler to make it GPL for the time being.

Cheers,

Henry
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Nathaniel Smith
On Tue, Nov 19, 2013 at 9:17 AM, Henry Gomersall <[hidden email]> wrote:

> On 19/11/13 16:08, Stéfan van der Walt wrote:
>> On Tue, Nov 19, 2013 at 6:03 PM, Henry Gomersall<[hidden email]>  wrote:
>>> >However, FFTW is dual licensed GPL/commercial and so the wrappers are
>>> >also GPL by necessity.
>> I'm not sure if that is true, strictly speaking--you may license your
>> wrapper code under any license you wish.  It's just that it becomes
>> confusing when the combined work then has to be released under GPL.
>
> This is on shaky GPL ground. I'm inclined to agree with you, but given
> any usage necessarily has to link against the FFTW libs, any resultant
> redistribution is going to be GPL by necessity. I.e. pyFFTW is useless
> without FFTW so it's just simpler to make it GPL for the time being.

The case where it makes a difference is when someone has purchased a
commercial license to FFTW and then wants to use it from their
proprietary Python application.

-n
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Henry Gomersall-2
On 19/11/13 17:52, Nathaniel Smith wrote:

> On Tue, Nov 19, 2013 at 9:17 AM, Henry Gomersall<[hidden email]>  wrote:
>> >On 19/11/13 16:08, Stéfan van der Walt wrote:
>>> >>On Tue, Nov 19, 2013 at 6:03 PM, Henry Gomersall<[hidden email]>   wrote:
>>>>> >>> >However, FFTW is dual licensed GPL/commercial and so the wrappers are
>>>>> >>> >also GPL by necessity.
>>> >>I'm not sure if that is true, strictly speaking--you may license your
>>> >>wrapper code under any license you wish.  It's just that it becomes
>>> >>confusing when the combined work then has to be released under GPL.
>> >
>> >This is on shaky GPL ground. I'm inclined to agree with you, but given
>> >any usage necessarily has to link against the FFTW libs, any resultant
>> >redistribution is going to be GPL by necessity. I.e. pyFFTW is useless
>> >without FFTW so it's just simpler to make it GPL for the time being.
> The case where it makes a difference is when someone has purchased a
> commercial license to FFTW and then wants to use it from their
> proprietary Python application.

Yes, this didn't occur to me as an option, mostly because I'm keen for a
commercial FFTW license myself and it would gall me somewhat if I
couldn't gain the same benefit from my own code as others.

So, given that, if anyone has an FFTW license and is keen for decent
Python wrappers, I'd be more than happy to discuss a sub-license to FFTW
in exchange for a more liberal (say MIT) license for pyFFTW.

Cheers,

Henry
_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
Reply | Threaded
Open this post in threaded view
|

Re: strange runtimes of numpy fft

Chris Barker - NOAA Federal
On Wed, Nov 20, 2013 at 3:06 AM, Henry Gomersall <[hidden email]> wrote:
Yes, this didn't occur to me as an option, mostly because I'm keen for a
commercial FFTW license myself and it would gall me somewhat if I
couldn't gain the same benefit from my own code as others.

So, given that, if anyone has an FFTW license and is keen for decent
Python wrappers, I'd be more than happy to discuss a sub-license to FFTW
in exchange for a more liberal (say MIT) license for pyFFTW.


OT, and IANAL, but  I think what you'd want to do is dual-licence pyFFTW. Heck, you could even charge for a commercially-licenced version, as it would only be useful to folks that were already paying for a commercially-licenced FFTW

-Chris



--

Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

[hidden email]

_______________________________________________
NumPy-Discussion mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
12