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

[Joseph Codadeen]

Ahhhhh, looking back I see what Pierre-André did, the penny drops. Sorry I miss-read.I'm only interested in this part;>zeropadded_fft_A = fft(seq_A, n=2**(ceil(log(len(seq_A),2))+1))
>zeropadded_fft_B = fft(seq_B, n=2**(ceil(log(len(seq_B),2))+1))>You could remove the "+1" above to get faster results, but then that may >lead to unwanted frequencies (coming from the fact that fft assumes >periodic signals, read online about zero-padding).
I did it in the following way and it seems to work. I've seen an increase in performance for sure.
My Python sucks by the way :)def pad_to_power_of_two(audio):
    '''
    Calculates next power of two and pads the audio sample array with zeros
    '''    # RIFF spec specifies: ckSize    A 32-bit unsigned value identifying the size of ckData.
    max_shift = 32
    shift = 1
    power = len(audio) - 1
    print "{0:b}".format(power)    # Test if a power of two. Turn off the right-most bit (x & (x - 1)) and zero test
    while ((power + 1) & power) and (shift <= max_shift):
        power |= power >> shift
        shift *= 2
        print  "new power: " + "{0:b}".format(power)
    power = power + 1
    print  "{0:b}".format(power)

    ext = np.array([0x00] * (power - len(audio)))
    new_array = np.concatenate((audio, ext), axis=0)

    print  "Next power of 2 greater than actual array length is " + str(power)
    print  "Extending current array length (" + str(len(audio)) + ") by " + str(power - len(audio)) + " null bytes"
    print  "New array length is " + str(len(new_array))

    return new_array
I will try using the suggested methods and see if I can boost the performance further.Thanks.
> From: oscar.j.benjamin at gmail.com
> Date: Tue, 1 Sep 2015 16:14:41 +0100
> To: numpy-discussion at scipy.org
> Subject: Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples
> 
> 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
> _______________________________________________
> 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/20150901/244e003b/attachment.html>