[Numpy-discussion] rfft different in numpy vs scipy
a.h.jaffe at gmail.com
Thu Sep 7 19:04:52 CDT 2006
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
> 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 == zero frequency, f + 1j* f 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.
This is complicated by the fact that (I think) numpy.fft shows up as
scipy.fft, so functions with the same name in scipy.fft and
scipy.fftpack aren't actually the same!
More information about the Numpy-discussion