Parseval's theorem Python NumPy example

readme.md

parseval.py shows simply how to do the calculation for Parseval's theorem with NumPy's FFT.

The main point is that you have to normalize by the number of samples (depending on your FFT implementation, probably).

Then I made functions.py to calculate RMS faster in the frequency domain and example.txt shows the time savings.

Note that there is no benefit to doing the calculation in the frequency domain if your signal is in the time domain. It's only beneficial when you've already converted to frequency domain and done some frequency-domain processing.

example.txt

In [1]: n = 10000000

In [2]: x = rand(n)

In [3]: X = fft(x)

In [4]: rms_flat(x)
Out[4]: 0.57731639149367697

In [5]: rms_flat(ifft(X))
Out[5]: 0.57731639149367675

In [6]: rms_fft(X)
Out[6]: 0.57731639149370662

In [7]: timeit rms_flat(ifft(X))
1 loops, best of 3: 3.11 s per loop

In [8]: timeit rms_fft(X)
1 loops, best of 3: 762 ms per loop

In [9]: rX = rfft(x)

In [10]: rms_rfft(rX)
Out[10]: 0.57731639149369829

In [11]: timeit rms_flat(irfft(rX))
1 loops, best of 3: 1.61 s per loop

In [12]: timeit rms_rfft(rX)
1 loops, best of 3: 379 ms per loop

parseval.py

from numpy.random import rand
from numpy.fft import fft
from parseval_functions import rms_flat, sqrt

n = 100
x = rand(n)
X = fft(x)

# total energy
print(sum(abs(x)**2))      # 33.803768126149471, for instance
print(sum(abs(X)**2)/n)    # 33.803768126149457, for instance

# root mean square
print(rms_flat(x))              # 0.56180320113171067, for instance
print(rms_flat(X)/sqrt(n))      # 0.56180320113171078, for instance
print(rms_flat(abs(X))/sqrt(n)) # 0.56180320113171078, for instance (faster if you've already computed abs(X))

parseval_functions.py

from numpy import sqrt, mean, absolute, real, conj


def rms_flat(a):
    """
    Return the root mean square of all the elements of *a*, flattened out.
    """
    return sqrt(mean(absolute(a)**2))


def rms_fft(spectrum):
    """
    Use Parseval's theorem to find the RMS value of a signal from its fft,
    without wasting time doing an inverse FFT.

    For a signal x, these should produce the same result, to within numerical
    accuracy:

    rms_flat(x) ~= rms_fft(fft(x))
    """
    return rms_flat(spectrum)/sqrt(len(spectrum))


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-length
        mean = (sq[0] + 2*sum(sq[1:])           )/n
    else:  # even-length
        mean = (sq[0] + 2*sum(sq[1:-1]) + sq[-1])/n
    root = sqrt(mean)
    return root/sqrt(n)


if __name__ == "__main__":
    from numpy.random import randn
    from numpy.fft import fft, ifft, rfft, irfft

    n = 17
    x = randn(n)
    X = fft(x)
    rX = rfft(x)

    print(rms_flat(x))
    print(rms_flat(ifft(X)))
    print(rms_fft(X))
    print()

    # Accurate for odd n:
    print(rms_flat(irfft(rX, n)))
    print(rms_rfft(rX, n))
    print()

    # Only approximate for odd n:
    print(rms_flat(irfft(rX)))
    print(rms_rfft(rX))