(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!

possible duplicate of [Parseval's theorem in Python](http://stackoverflow.com/questions/14011506/parsevals-theorem-in-python) — Jaime, Mar 01 '13 at 00:19

score 8 · Accepted Answer · answered Mar 01 '13 at 00:17

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.

Thanks Jamie! That was super helpful! I don't remember the weird normalization scheme in the documentation. Coming from IDL, I was used to having that done automatically by the fft routine. Now I know better in case I start using some other language at some point. — toylas, Mar 01 '13 at 03:35
It's not weird, it's the DFT. According to [this IDL documentation](http://www.exelisvis.com/docs/FFT.html), the full normalisation there is applied on forward operation, so you'd still get incorrect power measurements. — Henry Gomersall, Mar 01 '13 at 11:25
Apparently, the normalization of `rfft` has changed in the meantime. The given normalization no longer applies to numpy 1.8.2. — fuesika, Sep 01 '16 at 11:23
@fuesika I had a look at my old code and found myself dividing the result of `irfft()` by `n` to get the inverse transform result. So has that changed too? If so the answers should give some hint on version differences. — Jason, Oct 01 '16 at 08:48

Henry Gomersall · Answer 2 · 2013-03-01T11:28:10.477

3

In most FFT libraries, the various DFT flavours are not orthogonal. The numpy.fft library applies the necessary normalizations only during the inverse transform.

Consider the Wikipedia description of the DFT; the inverse DFT has the 1/N term that the DFT does not have (in which N is the length of the transform). To make an orthogonal version of the DFT, you need to scale the result of the un-normalised DFT by 1/sqrt(N). In this case, the transform is orthogonal (that is, if we define the orthogonal DFT as F, then the inverse DFT is the conjugate, or hermitian, transpose of F).

In your case, you can get the correct answer by simply scaling aft by 1.0/sqrt(len(a)) (note that N is found from the length of the transform; the real FFT just throws about half the values away, so it's the length of a that is important).

I suspect that the reason for leaving the normalization until the end is that in most situations, it doesn't matter and you therefore save the computational cost of doing the normalization twice. Indeed, the very quick FFTW library doesn't do any normalization in either direction, and leaves it entirely up to the user to deal with.

Edit: Just to be clear, the explanation above is not quite correct. The correct answer will not be arrived at with that simple scaling, as in your case the DC component will be added in twice, although 1.0/sqrt(len(a)) is still the correct scaling to produce the unitary transform.

edited Mar 01 '13 at 11:28

answered Feb 28 '13 at 23:42

Henry Gomersall

8,434
3
31
54

You are very right about the lack of a normalization factor. But your remarks on how to get the right result from `rfft` are very wrong, as simply trying it out will show. – Jaime Mar 01 '13 at 00:18
1

*very wrong* is disingenuous; the scaling is correct. I was wrong, but only where I said that scaling is sufficient. I suppose that would teach me to write SO answers with jet-lag. – Henry Gomersall Mar 01 '13 at 11:20
Please accept my apologies if my comment came in as too harsh. But unfortunately you are still wrong. `rfft`, apart from repeated terms excluded, and an almost 2x speed-up, returns the exact same you would get from `fft`. And it must therefore be scaled in the same way to get a unitary transform: by the square root of the length of the data, not of the unique values in the transform. [Parseval's theorem](http://en.wikipedia.org/wiki/Discrete_Fourier_transform#The_Plancherel_theorem_and_Parseval.27s_theorem) is the same regardless of the nature of the input. – Jaime Mar 01 '13 at 14:24
That's what I said; the length of the *transform* is what is important, and then I wrote it too to be clear "it's the length of `a` that is important". Since `a` is clearly the input vector in the original question, this is pretty unambiguous. – Henry Gomersall Mar 01 '13 at 14:42
It's worth point out the distinction between the length of the input or output vectors and the length of the transform, because, of course, in the case of the _inverse_ real DFT, it's the length of the *output* vector that is important for the normalisation. – Henry Gomersall Mar 01 '13 at 14:49
1

Please, accept my apologies: you have been right all along. I guess it was your "note that N is found from the length of the transform" that threw me off course. I will edit my answer, and have upvoted yours. Should I erase my comments, or would you rather have them stand as a testimony of my stupidity? – Jaime Mar 01 '13 at 16:13
1

No, of course not. We all learn from one another and I was as negligent as you :) – Henry Gomersall Mar 01 '13 at 17:08

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

2 Answers2