On Fri, Aug 28, 2015 at 11:02 AM, Joseph Codadeen wrote:

Hi,

I am a numpy newbie.

I have two wav files, one that numpy takes a long time to process the FFT. They was created within audacity using white noise and silence for gaps.

1. my_1_minute_noise_with_gaps.wav 2. my_1_minute_noise_with_gaps_truncated.wav

The files are very similar in the following way;

- 1. is white noise with silence gaps on every 15 second interval. - 2. is 1. but slightly shorter, i.e. I trimmed some ms off the end but it still has the last gap at 60s.

The code I am using processes the file like this;

framerate, data = scipy.io.wavfile.read(filepath) right = data[:, 0] # Align it to be efficient. if len(right) % 2 != 0: right = right[range(len(right) - 1)] noframes = len(right) fftout = np.fft.fft(right) / noframes # <<< I am timing this cmd

Using timeit...

- my_1_minute_noise_with_gaps_truncated took *30.75620985s* to process. - my_1_minute_noise_with_gaps took *22307.13917s* to process.

Could someone tell me why this behaviour is happening please?

Sorry I can't attach the files as this email gets bounced but you could easily create the files yourself. E.g my last gap width is 59.9995 - 1:00.0005, I repeat this every 15 seconds. My truncated file is 1:00.0015s long, non-truncated is 1:00.0833s long

It is almost certainly caused by the number of samples in your signals, i.e. look at what noframes is in one case and the other. You will get best performance when noframes is a power of two, or has a factorization that includes many small integers (2, 3, 5, perhaps also 7), and the worst if the size is a prime number. Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

On 2015-08-28 11:51:47, Joseph Codadeen wrote:

my_1_minute_noise_with_gaps_truncated - Array len is 2646070my_1_minute_noise_with_gaps - Array len is 2649674

In [6]: from sympy import factorint In [7]: max(factorint(2646070)) Out[7]: 367 In [8]: max(factorint(2649674)) Out[8]: 1324837 Those numbers give you some indication of how long the FFT will take to compute. Stéfan

If you don't mind the extra dependency or licensing and this is an issue for you, you can try pyfftw (there are likely other similar projects) which wraps fftw and does not have this problem as far as I know. It exposes a numpy-like interface. - sebastian On Fr, 2015-08-28 at 19:13 +0000, Joseph Codadeen wrote:

Great, thanks Stefan and everyone.

...

From: stefanv@berkeley.edu To: numpy-discussion@scipy.org Date: Fri, 28 Aug 2015 12:03:52 -0700 Subject: Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples

On 2015-08-28 11:51:47, Joseph Codadeen wrote:

...

my_1_minute_noise_with_gaps_truncated - Array len is 2646070my_1_minute_noise_with_gaps - Array len is 2649674

In [6]: from sympy import factorint In [7]: max(factorint(2646070)) Out[7]: 367 In [8]: max(factorint(2649674)) Out[8]: 1324837

Those numbers give you some indication of how long the FFT will take to compute.

Stéfan _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

If your sequence is not meant to be periodic (i.e., if after one minute there is no reason why the signal should start back at the beginning right away), then you should do zero-padding. And while you zero-pad, you can zero-pad to a sequence that is a power of two, thus preventing awkward factorizations. from numpy.fft import fft from numpy.random import rand from math import log, ceil seq_A = rand(2649674) seq_B = rand(2646070) fft_A = fft(seq_A) #Long fft_B = fft(seq_B) 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). Have a nice day, Pierre-André On 08/28/2015 12:03 PM, Stefan van der Walt wrote:

On 2015-08-28 11:51:47, Joseph Codadeen wrote:

...

my_1_minute_noise_with_gaps_truncated - Array len is 2646070my_1_minute_noise_with_gaps - Array len is 2649674

In [6]: from sympy import factorint In [7]: max(factorint(2646070)) Out[7]: 367 In [8]: max(factorint(2649674)) Out[8]: 1324837 Those numbers give you some indication of how long the FFT will take to compute.

Stéfan

Zero-padding won't help with the non-periodicity, will it? For that you may want to window instead.

Umh, it depends what you use the FFT for. You are right Stéfan when saying that Joseph should probably also use a window to get rid of the high frequencies that will come from the sharp steps at the beginning and end of his signal. I had in mind the use of FFT to do convolutions ( https://en.wikipedia.org/wiki/Convolution_theorem ). If you do not zero-pad properly, then the end of the signal may "bleed" on the beginning, and vice versa.

Hi, I cannot see how the following would work when it is np.fft.fft() that takes a long time based on the length of data. In my case my data is non-periodic.> from numpy.fft import fft> from numpy.random import rand> from math import log, ceil> seq_A = rand(2649674)> seq_B = rand(2646070)> fft_A = fft(seq_A) #Long> fft_B = fft(seq_B)>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))Ideally I need to calculate a favourable length to pad, prior to calling np.fft.fft() as Stéfan pointed out by calculating the factors. In [6]: from sympy import factorint In [7]: max(factorint(2646070)) Out[7]: 367 In [8]: max(factorint(2649674)) Out[8]: 1324837 I will try adding some code to calculate the next power of two above my array length as you suggest;> 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.Thanks,Joseph Date: Fri, 28 Aug 2015 17:26:32 -0700 From: stefanv@berkeley.edu To: numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] Numpy FFT.FFT slow with certain samples On Aug 28, 2015 5:17 PM, "Pierre-Andre Noel" wrote:

I had in mind the use of FFT to do convolutions (

https://en.wikipedia.org/wiki/Convolution_theorem ). If you do not

zero-pad properly, then the end of the signal may "bleed" on the

beginning, and vice versa. Ah, gotcha! All these things should also be handled nicely in scipy.signal.fftconvolve. Stéfan

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

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 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

On Tue, 1 Sep 2015 18:43 Phil Hodge wrote: On 09/01/2015 11:14 AM, Oscar Benjamin wrote:

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.

It would be better to pad to a power of three or five, or some other product of small primes, if the pad length would be significantly less than padding to a power of two. The obvious problem with padding is that it's different from the real data, so its presence will affect the result. At the least, it would be good to pad with the mean value of the original data, or to pad with zero after subtracting the mean from the original data. I meant performance wise it's not a bad choice. If you're concerned about distortion then use the Bluestein algorithm as I showed. -- Oscar