[Numpy-discussion] rfft different in numpy vs scipy

Charles R Harris charlesr.harris at gmail.com
Thu Sep 7 19:42:08 CDT 2006


On 9/7/06, Andrew Jaffe <a.h.jaffe at gmail.com> wrote:
>
> Hi Charles,
>
> Charles R Harris wrote:
> > On 9/7/06, *Andrew Jaffe* <a.h.jaffe at gmail.com
> > <mailto:a.h.jaffe at gmail.com>> wrote:
> >
> >     Hi all,
> >
> >     It seems 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.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?
> >
> >     Anyone else have an opinion?
> >
> > Yes, I prefer the scipy version because the result is actually a complex
> > array and can immediately be use as the coefficients of the fft for
> > frequencies <= Nyquist. I suspect, without checking, that what you get
> > in numpy is a real array with f[0] == zero frequency, f[1] + 1j* f[2] as
> > the coefficient of the second frequency, etc. This makes it difficult to
> > multiply by a complex transfer function or phase shift the result to
> > rotate the original points by some fractional amount.
> >
> > As to unpacking, for some algorithms the two real coefficients are
> > packed into the real and complex parts of the zero frequency so all that
> > is needed is an extra complex slot at the end. Other algorithms produce
> > what you describe. I just think it is more convenient to think of the
> > real fft as an efficient complex fft that only computes the coefficients
> > <= Nyquist because Hermitean symmetry determines the rest.
>
> Unless I misunderstand, I think you've got it backwards:
>    - numpy.fft.rfft produces the correct complex array (with no strange
> packing), with frequencies (0, 1, 2, ... n/2)
>
>    - scipy.fftpack.rfft produces a single real array, in the correct
> order, but with frequencies (0, 1, 1, 2, 2, ..., n/2) -- as given by
> scipy.fftpack's rfftfreq function.
>
> So I think you prefer numpy, not scipy.


Ah, well then, yes. I prefer Numpy. IIRC, one way to get the Scipy ordering
is to use the Hartley transform as the front to the real transform. And, now
that you mention it, there was a big discussion about it in the scipy list
way back when with yours truly pushing the complex form. I don't recall that
anything was settled, rather the natural output of the algorithm they were
using prevailed by default. Maybe you could write a front end that did the
right thing?

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/numpy-discussion/attachments/20060907/c5d1abf1/attachment-0001.html 


More information about the Numpy-discussion mailing list