Numpy's real fft (rfft) - losing power

Question

Related to another problem I'm having, I was looking into the workings of numpy's rfft2 and irfft2. These are special versions of the FFT routine, in so far that it needs less input; because you require the real-space image to be real you only need to 'fill' half of Fourier space - due to symmetry, that's all the information you need.

Now, in the docs it is mentioned that performing the inverse transform right after the forward transform should return the original array - to within numerical precision. Time to put this to the test.

If I do the following:

import numpy as np

fourier_image = np.random.uniform(-1,1,(100,51))
fourier_recovered = np.fft.rfft2(np.fft.irfft2(fourier_image))

print 'Original power: ', sum(sum(abs(fourier_image)**2.))
print 'Recovered power: ', sum(sum(abs(fourier_recovered)**2.))

... I find the original power and the recovered power differ at the percent level. That seems a rather big difference to me. Am I doing something wrong, or is this just the best we can do?

Answer 1

I naively assumed any (N,N/2+1) Fourier image had the right symmetry, but this goes wrong at the 0- and Nyquist frequency because they don't have a symmetric counterpart. Here, you have to impose the symmetry explicitly.

So, if we generate an image:

N = 100
fourier_image = np.random.uniform(-1,1,(N,1+N/2))

we impose the required symmetry as follows:

fourier_image[1:N/2,0] = np.conj(fourier_image[N:N/2:-1,0])
fourier_image[1:N/2,-1] = np.conj(fourier_image[N:N/2:-1,-1])

And doing the same as in the question above, we find the input - and recovered power:

print 'Original power: ', sum(sum(abs(fourier_image)**2.))
print 'Recovered power: ', sum(sum(abs(fourier_recovered)**2.))

which are indeed the same:

Original power: 1664.3177945
Recovered power: 1664.3177945

For an odd-shaped image this will be slightly different but the idea is the same.