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
import numpy as N
from numpy.testing import *
import timeit
def fftfreq0(n,d=1.0):
""" fftfreq(n, d=1.0) -> f
DFT sample frequencies
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,...,n/2-1,-n/2,...,-1]/(d*n) if n is even
f = [0,1,...,(n-1)/2,-(n-1)/2,...,-1]/(d*n) if n is odd
"""
assert isinstance(n,int) or isinstance(n,integer)
k = range(0,(n-1)/2+1)+range(-(n/2),0)
return N.array(k,'d')/(n*d)
def fftfreq1(n,d=1.0):
""" fftfreq(n, d=1.0) -> f
DFT sample frequencies
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,...,n/2-1,-n/2,...,-1]/(d*n) if n is even
f = [0,1,...,(n-1)/2,-(n-1)/2,...,-1]/(d*n) if n is odd
"""
assert isinstance(n,int) or isinstance(n,integer)
k = N.hstack((N.arange(0,(n-1)/2 + 1), N.arange(-(n/2),0))) / (n*d)
return k
def fftfreq2(n,d=1.0):
""" fftfreq(n, d=1.0) -> f
DFT sample frequencies
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,...,n/2-1,-n/2,...,-1]/(d*n) if n is even
f = [0,1,...,(n-1)/2,-(n-1)/2,...,-1]/(d*n) if n is odd
"""
assert isinstance(n,int) or isinstance(n,integer)
k = N.empty(n)
midpoint = (n-1)/2+1
k[:midpoint] = N.arange(0,(n-1)/2 + 1)
k[midpoint:] = N.arange(-(n/2),0)
k *= 1./(n*d)
return k
for i in [int(x) for x in 1e5,1e5+1,1e6,1e6+1]:
print "Benchmarking for n=%d" % i
def bench(fname,out="x"):
return timeit.Timer("__main__.%s=__main__.%s(%d)" % (out,fname,i),
"import __main__").timeit(number=10)
print "Old: ", bench("fftfreq0",out="a")
print "New_concat: ", bench("fftfreq1",out="b")
print "New_inplace: ", bench("fftfreq2",out="c")
print
assert_array_almost_equal(a,b)
assert_array_almost_equal(b,c)
-------------------------------------------------------------------------
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
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion
