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, non-truncated 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 non-truncated 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: numpy-discussion@scipy.org Subject: Re: [Numpy-discussion] 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 non-truncated one had a length that was either prime or could only be factored into large primes.
Phil _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On 2015-08-28 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: 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 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 _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
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 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 _______________________________________________ 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
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 numpy-like 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, 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 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 _______________________________________________ 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
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 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 2015-08-28 16:20:33, Pierre-Andre 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 zero-padding. And while you zero-pad, you can zero-pad to a sequence that is a power of two, thus preventing awkward factorizations.
Zero-padding won't help with the non-periodicity, will it? For that you may want to window instead.
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.
On Aug 28, 2015 5:17 PM, "Pierre-Andre 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 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
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" 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
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
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
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
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
NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
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, 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