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. my_1_minute_noise_with_gaps.wavmy_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, nontruncated is 1:00.0833s long
Thanks.
On 08/28/2015 02:02 PM, Joseph Codadeen wrote:
 my_1_minute_noise_with_gaps_truncated took***30.75620985s* to process.
 my_1_minute_noise_with_gaps took *22307.13917s*to process.
You didn't say how long those arrays were, but I can make a good guess that the truncated one had a length that could be factored into small, prime numbers, while the nontruncated one had a length that was either prime or could only be factored into large primes.
Phil
my_1_minute_noise_with_gaps_truncated  Array len is 2646070my_1_minute_noise_with_gaps  Array len is 2649674
Date: Fri, 28 Aug 2015 14:28:49 0400 From: hodge@stsci.edu To: numpydiscussion@scipy.org Subject: Re: [Numpydiscussion] Numpy FFT.FFT slow with certain samples
On 08/28/2015 02:02 PM, Joseph Codadeen wrote:
 my_1_minute_noise_with_gaps_truncated took***30.75620985s* to process.
 my_1_minute_noise_with_gaps took *22307.13917s*to process.
You didn't say how long those arrays were, but I can make a good guess that the truncated one had a length that could be factored into small, prime numbers, while the nontruncated one had a length that was either prime or could only be factored into large primes.
Phil _______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
On 20150828 11:51:47, Joseph Codadeen jdmc80@hotmail.com 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
Great, thanks Stefan and everyone.
From: stefanv@berkeley.edu To: numpydiscussion@scipy.org Date: Fri, 28 Aug 2015 12:03:52 0700 Subject: Re: [Numpydiscussion] Numpy FFT.FFT slow with certain samples
On 20150828 11:51:47, Joseph Codadeen jdmc80@hotmail.com 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 _______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
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 numpylike interface.
 sebastian
On Fr, 20150828 at 19:13 +0000, Joseph Codadeen wrote:
Great, thanks Stefan and everyone.
From: stefanv@berkeley.edu To: numpydiscussion@scipy.org Date: Fri, 28 Aug 2015 12:03:52 0700 Subject: Re: [Numpydiscussion] Numpy FFT.FFT slow with certain
samples
On 20150828 11:51:47, Joseph Codadeen jdmc80@hotmail.com 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 _______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
On 2015/08/28 10:36 AM, Sebastian Berg wrote:
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 numpylike interface.
Sort of; that interface returns a function, not the result.
fftw is still an fft algorithm, so it is still subject to a huge difference in run time depending on how the input array can be factored.
Furthermore, it gets its speed by figuring out how to optimize a calculation for a given size of input array. That initial optimization can be very slow. The overall speed gain is realized only when one saves the result of that optimization, and applies it to many calculations on arrays of the same size.
Eric
 sebastian
On Fr, 20150828 at 19:13 +0000, Joseph Codadeen wrote:
Great, thanks Stefan and everyone.
From: stefanv@berkeley.edu To: numpydiscussion@scipy.org Date: Fri, 28 Aug 2015 12:03:52 0700 Subject: Re: [Numpydiscussion] Numpy FFT.FFT slow with certain
samples
On 20150828 11:51:47, Joseph Codadeen jdmc80@hotmail.com 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 _______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
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 zeropadding. And while you zeropad, you can zeropad 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 zeropadding).
Have a nice day,
PierreAndré
On 08/28/2015 12:03 PM, Stefan van der Walt wrote:
On 20150828 11:51:47, Joseph Codadeen jdmc80@hotmail.com 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
On 20150828 16:20:33, PierreAndre Noel noel.pierre.andre@gmail.com wrote:
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 zeropadding. And while you zeropad, you can zeropad to a sequence that is a power of two, thus preventing awkward factorizations.
Zeropadding won't help with the nonperiodicity, will it? For that you may want to window instead.
Stéfan
Zeropadding won't help with the nonperiodicity, 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 zeropad properly, then the end of the signal may "bleed" on the beginning, and vice versa.
On Aug 28, 2015 5:17 PM, "PierreAndre Noel" noel.pierre.andre@gmail.com wrote:
I had in mind the use of FFT to do convolutions ( https://en.wikipedia.org/wiki/Convolution_theorem ). If you do not zeropad 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
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 nonperiodic.> 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 zeropad, you can zeropad 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: numpydiscussion@scipy.org Subject: Re: [Numpydiscussion] Numpy FFT.FFT slow with certain samples
On Aug 28, 2015 5:17 PM, "PierreAndre Noel" noel.pierre.andre@gmail.com wrote:
I had in mind the use of FFT to do convolutions (
https://en.wikipedia.org/wiki/Convolution_theorem ). If you do not
zeropad 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
_______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
On 1 September 2015 at 11:38, Joseph Codadeen jdmc80@hotmail.com wrote:
And while you zeropad, you can zeropad 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 zeropadding distorting the FFT you can use the Bluestein transform as below. This pads up to a power of two (greater than 2N1) 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
Ahhhhh, looking back I see what PierreAndré did, the penny drops. Sorry I missread.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 zeropadding).
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 32bit 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 rightmost 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: numpydiscussion@scipy.org Subject: Re: [Numpydiscussion] Numpy FFT.FFT slow with certain samples
On 1 September 2015 at 11:38, Joseph Codadeen jdmc80@hotmail.com wrote:
And while you zeropad, you can zeropad 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 zeropadding distorting the FFT you can use the Bluestein transform as below. This pads up to a power of two (greater than 2N1) 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 _______________________________________________ NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
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.
Phil
On Tue, 1 Sep 2015 18:43 Phil Hodge hodge@stsci.edu 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
On Tue, Sep 1, 2015 at 12:06 PM, Oscar Benjamin oscar.j.benjamin@gmail.com wrote:
On Tue, 1 Sep 2015 18:43 Phil Hodge hodge@stsci.edu 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.
I don't see the problem with padding. After all, the spectrum is an estimate and the effect of the padding, as well as the use of windowing is well understood. However, I would recommend subtracting the average before padding, as otherwise the DC frequency is likely to be of large amplitude and its side lobes will mask the lower frequencies.
 Oscar
NumPyDiscussion mailing list NumPyDiscussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpydiscussion
On Fri, Aug 28, 2015 at 11:02 AM, Joseph Codadeen jdmc80@hotmail.com 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.
 my_1_minute_noise_with_gaps.wav
 my_1_minute_noise_with_gaps_truncated.wav
The files are very similar in the following way;
 is white noise with silence gaps on every 15 second interval.
 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, nontruncated 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
participants (10)

Charles R Harris

Eric Firing

Jaime Fernández del Río

Joseph Codadeen

Oscar Benjamin

Phil Hodge

PierreAndre Noel

Sebastian Berg

Stefan van der Walt

Stéfan van der Walt