[Numpy-discussion] rfft different in numpy vs scipy

[Steven G. Johnson]

Andrew Jaffe wrote:
> 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.e., two real numbers at the same frequency, except for
> the highest and lowest) [All of the above for even n; but the difference
> between numpy and scipy remains for odd n.]
>
> 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?

As an author of FFTW, let me interject a couple of points into this
discussion.

First, if you are using FFTW, then its real-input r2c routines
"natively" produce output in the "unpacked" numpy format as described
above: an array of n/2+1 complex numbers.  Any "packed" format would
require some data permutations.  Other FFT implementations use a
variety of formats.

Second, the *reason* why FFTW's r2c routines produce unpacked output is
largely because "packed" formats do not generalize well to
multi-dimensional FFTs, while the "unpacked" format does.  (Packed
formats are *possible* for multidimensional transforms, but become
increasingly intricate as you add more dimensions.)  Additionally, I
personally find the unpacked format more convenient in most
applications.

(FFTW does actually support a different, "packed" format in its r2r
interface, but only for 1d FFTs.  The reason for this has do to with
advantages of that format for odd sizes.  We recommend that most users
employ the r2c interface, however; our r2c interface is generally
faster for even sizes.)

I hope this is helpful.

Cordially,
Steven G. Johnson