I forgot something important : I'm using numpy 1.8.2 and python 3.4.1 on linux (fedora 21)

I beleive numpy doesn't support float32 for fft. If you want that, you will need to use Scipy, although, IIRC, there may be some problems there also due to inaccuracy in the transform of data whose length has large prime factors. Probably fixable, but I don't think it is fixed yet.

It's okay with scipy.fftpack rfft, but... while digging the problem, I've also found some similar behaviors in scipy.interpolate, for instance.
I thought that numpy and scipy both transparently supported float32 and complex64. Am I wrong?

Numpy is more basic than Scipy. The ufuncs suppport all the relevant types, as do most of the linalg routines these days, but in general Scipy is more complete. Numpy support for all the types does evolve but it isn't a high priority.

And scipy.fftpack lacks rfft2, whose convenience I was enjoying in numpy. :(

For anyone coming from Google, here's how to replace numpy's rfft2 with scipy's 1D functions that works for me, specifically the rfft2 and irfft2 lambdas.

import numpy.fft as fft
import scipy.fftpack as scifft
import numpy as np

x = np.random.randn(10, 10).astype(np.float32)

# Using two invocations of scipy.fftpack.rfft & irfft. Adjust your call to these
# functions with extra arguments (n, overwrite_x, etc.) as needed.
rfft2 = lambda x: scifft.rfft(scifft.rfft(x, axis=1), axis=0)
irfft2 = lambda X: scifft.irfft(scifft.irfft(X, axis=0))

xScipy = irfft2(rfft2(x))

# Using numpy.fft.rfft2 & irfft2
xNumpy = fft.irfft2(fft.rfft2(x))

# Verify all equal. Note that the scipy version is float32, while the numpy
# version was implictly upcast to float64. For this reason, `allclose` needs a
# looser absolute tolerance, more in line with single-precision.
assert (np.allclose(x, xScipy, atol=2e-7))
assert (np.allclose(x, xNumpy))
print 'Tests pass. Scipy type: %s, Numpy type: %s' % (xScipy.dtype,
                                                      xNumpy.dtype)

(I realize this test just shows that the forward- and backwards-transforms are inverses of each other, and not that the correct frequency domain content is generated—but I've verified that frequency-domain filtering works using both versions here.)

This is also discussed in gh-17801.

Even though this issue is the older one, I'm closing it as the duplicate, because the discussion in gh-17801 is probably more up-to-date.