Instantly share code, notes, and snippets.

# endolith/example.txt

Last active Oct 3, 2022
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.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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))
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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))