Parseval's theorem in Python

Question

I'm trying to get some grip on Python's fft functionality, and one of the weird things that I've stumbled on is that Parseval's theorem doesn't seem to apply, as it gives a difference of about 50 now, while it should be 0.

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as fftpack

pi = np.pi

tdata = np.arange(5999.)/300
dt = tdata[1]-tdata[0]

datay = np.sin(pi*tdata)+2*np.sin(pi*2*tdata)
N = len(datay)

fouriery = abs(fftpack.rfft(datay))/N

freqs = fftpack.rfftfreq(len(datay), d=(tdata[1]-tdata[0]))

df = freqs[1] - freqs[0]

parceval = sum(datay**2)*dt - sum(fouriery**2)*df
print parceval

plt.plot(freqs, fouriery, 'b-')
plt.xlim(0,3)
plt.show()

I'm pretty sure that it's a normalisation factor, but I don't seem to be able to find it, as all the information I can find about this function is the scipy.fftpack.rfft documentation.

Answer 1

Your normalization factor is coming from trying to apply Parseval's theorem for the Fourier transform of a continuous signal to a discrete sequence. On the side panel of the wikipedia article on the Discrete Fourier transform there is some discussion on the relationship of the Fourier transform, the Fourier series, the Discrete Fourier Transform and sampling with Dirac combs.

To make a long story short, Parseval's theorem, when applied to DFTs, doesn't require integration, but summation: a 2*pi you are creating by multipliying by dt and df your summations.

Note also that, because you are using scipy.fftpack.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. So since you are only adding half the data, plus the 0in the DC term, there's the missing 2 to get to the 4*pi that @unutbu found.

In any case, if datay holds your sequence, you can verify Parseval's theorem as follows:

fouriery = fftpack.rfft(datay)
N = len(datay)
parseval_1 = np.sum(datay**2)
parseval_2 = (fouriery[0]**2 + 2 * np.sum(fouriery[1:]**2)) / N
print parseval_1 - parseval_2

Using scipy.fftpack.fft or numpy.fft.fft the second summation does not need to take on such a weird form:

fouriery_1 = fftpack.fft(datay)
fouriery_2 = np.fft.fft(datay)
N = len(datay)
parseval_1 = np.sum(datay**2)
parseval_2_1 = np.sum(np.abs(fouriery_1)**2) / N
parseval_2_2 = np.sum(np.abs(fouriery_2)**2) / N
print parseval_1 - parseval_2_1
print parseval_1 - parseval_2_2