Parseval's theorem - Average Power I want to calculate averaged power of a time-domain signal by means of its spectrum. I guess Parseval is the right tool.

So I sample a sinus of 100 Hz 10000x within one second.

Unfortunately the sum of the squared samples euqals not the sum of the FFT amplitudes (weighted by the number of FFT bins). Where is the mistake?

# Some python code

import numpy
import matplotlib.pyplot as plt

# Create Time Domain Signal for 1 sec
fs = n = 10000    # Samplingfrequency
ti = numpy.linspace(0,1,num=fs)
sx = 1*numpy.sin(2*numpy.pi*100*ti)    

# Calculate spectrum via FFT and account for scaling n/2
# taking the real fft (rfft) only the positive frequencies are calculated
fx = numpy.fft.rfft(sx)/(n/2)
no_of_points = fx.shape[0]

# Calculate RMS for time domains signal + spectrum
parseval_sx = numpy.sum(sx**2)
parseval_fx = numpy.sum(numpy.abs(fx)**2)/no_of_points

print parseval_sx, " equals not ", parseval_fx


Output:

4999.5  equals not  0.000199940012002
It's possible to calculate average power (proportional to RMS amplitude) with the rfft, but it's more complicated because the DC and Nyquist bins are not doubled in the full fft (and the Nyquist bin doesn't exist if n is odd).

NumPy's rfft is preferred over SciPy's, because it outputs complex numbers, which are compatible with N-dimensional arrays, rather than SciPy's "packed" format.

def rms_rfft(spectrum, n=None):
    """
    Use Parseval's theorem to find the RMS value of an even-length signal
    from its rfft, without wasting time doing an inverse real FFT.

    spectrum is produced as spectrum = numpy.fft.rfft(signal)

    For a signal x with an even number of samples, these should produce the
    same result, to within numerical accuracy:

    rms_flat(x) ~= rms_rfft(rfft(x))

    If len(x) is odd, n must be included, or the result will only be
    approximate, due to the ambiguity of rfft for odd lengths.
    """
    if n is None:
        n = (len(spectrum) - 1) * 2
    sq = real(spectrum * conj(spectrum))
    if n % 2:  # odd
        mean = (sq[0] + 2*sum(sq[1:])           )/n
    else:
        mean = (sq[0] + 2*sum(sq[1:-1]) + sq[-1])/n
    root = sqrt(mean)
    return root/sqrt(n)


More details here: https://gist.github.com/endolith/1257010