fftfreq very slow; rfftfreq incorrect?
Hi all, the current implementation of fftfreq (which is meant to return the appropriate frequencies for an FFT) does the following: k = range(0,(n-1)/2+1)+range(-(n/2),0) return array(k,'d')/(n*d) I have tried this with very long (2**24) arrays, and it is ridiculously slow. Should this instead use arange (or linspace?) and concatenate rather than converting the above list? This seems to result in acceptable performance, but we could also perhaps even pre-allocate the space. The numpy.fft.rfftfreq seems just plain incorrect to me. It seems to produce lots of duplicated frequencies, contrary to the actual output of rfft: def rfftfreq(n,d=1.0): """ rfftfreq(n, d=1.0) -> f DFT sample frequencies (for usage with rfft,irfft). The returned float array contains the frequency bins in cycles/unit (with zero at the start) given a window length n and a sample spacing d: f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2]/(d*n) if n is even f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2,n/2]/(d*n) if n is odd **** None of these should be doubled, right? """ assert isinstance(n,int) return array(range(1,n+1),dtype=int)/2/float(n*d) Thanks, Andrew
[copied to the scipy list since rfftfreq is only in scipy] Andrew Jaffe wrote:
Hi all,
the current implementation of fftfreq (which is meant to return the appropriate frequencies for an FFT) does the following:
k = range(0,(n-1)/2+1)+range(-(n/2),0) return array(k,'d')/(n*d)
I have tried this with very long (2**24) arrays, and it is ridiculously slow. Should this instead use arange (or linspace?) and concatenate rather than converting the above list? This seems to result in acceptable performance, but we could also perhaps even pre-allocate the space.
The numpy.fft.rfftfreq seems just plain incorrect to me. It seems to produce lots of duplicated frequencies, contrary to the actual output of rfft:
def rfftfreq(n,d=1.0): """ rfftfreq(n, d=1.0) -> f
DFT sample frequencies (for usage with rfft,irfft).
The returned float array contains the frequency bins in cycles/unit (with zero at the start) given a window length n and a sample spacing d:
f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2]/(d*n) if n is even f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2,n/2]/(d*n) if n is odd
**** None of these should be doubled, right?
""" assert isinstance(n,int) return array(range(1,n+1),dtype=int)/2/float(n*d)
Thanks,
Andrew
------------------------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642
Hello, I would like to discuss the following code: #***start*** import numpy as N a = N.array((200), dtype = N.uint8) print (a * 100) / 100 b = N.array((200, 200), dtype = N.uint8) print (b * 100) / 100 #***stop*** The first print statement will print "200" because the uint8-value is cast "upwards", I suppose. The second statement prints "[0 0]". I suppose this is due to overflows during the calculation. How can I tell numpy to do the upcast also in the second case, returning "[200 200]"? I am interested in the fastest solution regarding execution time. In my application I would like to store the result in an Numeric.UInt8-array. Thanks for every comment Lars
On 8/30/06, Lars Friedrich <lfriedri@imtek.de> wrote:
Hello,
I would like to discuss the following code:
#***start*** import numpy as N
a = N.array((200), dtype = N.uint8) print (a * 100) / 100
This is actually a scalar, i.e., a zero dimensional array. N.uint8(200) would give you the same thing, because (200) is a number, not a tuple like (200,). In any case In [44]:a = array([200], dtype=uint8) In [45]:a*100 Out[45]:array([32], dtype=uint8) In [46]:uint8(100)*100 Out[46]:10000 i.e. , the array arithmetic is carried out in mod 256 because Numpy keeps the array type when multiplying by scalars. On the other hand, when multiplying a *scalar* by a number, the lower precision scalars are upconverted in the conventional way. Numpy makes the choices it does for space efficiency. If you want to work in uint8 you don't have to recast every time you multiply by a small integer. I suppose one could demand using uint8(1) instead of 1, but the latter is more convenient. Integers can be tricky once the ordinary precision is exceeded and modular arithmetic takes over, it just happens more easily for uint8 than for uint32. Chuck
On 8/30/06, Lars Friedrich <lfriedri@imtek.de> wrote:
Hello,
I would like to discuss the following code:
#***start*** import numpy as N
a = N.array((200), dtype = N.uint8) print (a * 100) / 100
b = N.array((200, 200), dtype = N.uint8) print (b * 100) / 100 #***stop***
The first print statement will print "200" because the uint8-value is cast "upwards", I suppose. The second statement prints "[0 0]". I suppose this is due to overflows during the calculation.
How can I tell numpy to do the upcast also in the second case, returning "[200 200]"? I am interested in the fastest solution regarding execution time. In my application I would like to store the result in an Numeric.UInt8-array.
Thanks for every comment
To answer the original question, you need to use a higher precision array or explicitly cast it to higher precision. In [49]:(a.astype(int)*100)/100 Out[49]:array([200]) Chuck
On Wed, Aug 30, 2006 at 12:04:22PM +0100, Andrew Jaffe wrote:
the current implementation of fftfreq (which is meant to return the appropriate frequencies for an FFT) does the following:
k = range(0,(n-1)/2+1)+range(-(n/2),0) return array(k,'d')/(n*d)
I have tried this with very long (2**24) arrays, and it is ridiculously slow. Should this instead use arange (or linspace?) and concatenate rather than converting the above list? This seems to result in acceptable performance, but we could also perhaps even pre-allocate the space.
Please try the attached benchmark.
The numpy.fft.rfftfreq seems just plain incorrect to me. It seems to produce lots of duplicated frequencies, contrary to the actual output of rfft:
def rfftfreq(n,d=1.0): """ rfftfreq(n, d=1.0) -> f
DFT sample frequencies (for usage with rfft,irfft).
The returned float array contains the frequency bins in cycles/unit (with zero at the start) given a window length n and a sample spacing d:
f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2]/(d*n) if n is even f = [0,1,1,2,2,...,n/2-1,n/2-1,n/2,n/2]/(d*n) if n is odd
**** None of these should be doubled, right?
""" assert isinstance(n,int) return array(range(1,n+1),dtype=int)/2/float(n*d)
Please produce a code snippet to demonstrate the problem. We can then fix the bug and use your code as a unit test. Regards Stéfan
participants (4)
-
Andrew Jaffe -
Charles R Harris -
Lars Friedrich -
Stefan van der Walt