Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Parseval's theorem Python NumPy example

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.

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
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))
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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment