(numpy) Wrong amplitude(?) of FFT'd array?

Question

I'm using numpy and matplotlib to analyze data output form my simulations. There is one (apparent) inconsistency that I can't find the roots of. It's the following:

I have a signal that has a given energy a^2~1. When I use rfft to take the FFT and compute the energy in the Fourier space, it comes out to be significantly larger. To void giving the details of my data etc., here is an example with a simple sin wave:

from pylab import *
xx=np.linspace(0.,2*pi,128)
a=np.zeros(128)
for i in range(0,128):
    a[i]=sin(xx[i])
aft=rfft(a)
print mean(abs(aft)**2),mean(a**2) 

In principle both the numbers should be the same (at least in the numerical sense) but this is what I get out of this code:

62.523081632 0.49609375

I tried to go through numpy.fft documentation but could not find anything. A search here gave the following but I was not able to understand the explanations there:

Big FFT amplitude difference between the existing (synthesized) signal and the filtered signal

What am I missing/ misunderstanding? Any help/ pointer in this regard would be greatly appreciated.

Thanks!

Answer 1

Henry is right on the non-normalization part, but there is a little more to it, because you are using rfft, not fft. The following is consistent with his answer:

>>> x = np.linspace(0, 2 * np.pi, 128)
>>> y = 1 - np.sin(x)
>>> fft = np.fft.fft(y)
>>> np.mean((fft * fft.conj()).real)
191.49999999999991
>>> np.mean(y**2)
1.4960937500000004
>>> fft = fft / np.sqrt(len(fft))
>>> np.mean((fft * fft.conj()).real)
1.4960937499999991

But if you now try the same with rfft, things don't quite work out:

>>> rfft = np.fft.rfft(y)
>>> np.mean((rfft * rfft.conj()).real)
314.58462009358772
>>> rfft /= np.sqrt(len(rfft))
>>> np.mean((rfft * rfft.conj()).real)
4.8397633860551954
65
>>> np.mean((rfft * rfft.conj()).real) / len(rfft)
4.8397633860551954

The following does work properly, though:

>>> (rfft[0] * rfft[0].conj() +
...  2 * np.sum(rfft[1:] * rfft[1:].conj())).real / len(y)
1.4960937873636722

When you use rfft what you are getting is not properly the DFT of your data, but only the positive half of it, since the negative would be symmetric to it. To compute the mean, you need to consider every value other than the DC component twice, which is what the last line of code does.