Stefan van der Walt wrote:

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.

Results attached. Bottom line is that both new versions are several times faster than the old, and the concat (hstack) version is somewhat faster than the other, but it depends on n. I'm on OSX with the 2.4.3 universal build on PPC and the latest SVN numpy.

...

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:

<< removed >>

Please produce a code snippet to demonstrate the problem. We can then fix the bug and use your code as a unit test.

Aha, the problem is that scipy and numpy define rfft differently! numpy returns n/2+1 complex numbers (so the first and last numbers are actually real) with the frequencies equivalent to the positive part of the fftfreq, whereas scipy returns n real numbers with the frequencies as in rfftfreq. I think the numpy behavior makes more sense, as it doesn't require any unpacking after the fact, at the expense of a tiny amount of wasted space. But would this in fact require scipy doing extra work from whatever the 'native' real_fft (fftw, I assume) produces? Anyone else have an opinion?