`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@gmail.com

> Date: Tue, 1 Sep 2015 16:14:41 +0100

> To: numpy-discussion@scipy.org

> Subject: Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

>

> On 1 September 2015 at 11:38, Joseph Codadeen <jdmc80@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@scipy.org

> http://mail.scipy.org/mailman/listinfo/numpy-discussion

> Date: Tue, 1 Sep 2015 16:14:41 +0100

> To: numpy-discussion@scipy.org

> Subject: Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

>

> On 1 September 2015 at 11:38, Joseph Codadeen <jdmc80@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@scipy.org

> http://mail.scipy.org/mailman/listinfo/numpy-discussion