Skip to content

Instantly share code, notes, and snippets.

@SchofieChen
Forked from endolith/example.txt
Created December 16, 2020 16:20
Show Gist options
  • Save SchofieChen/c311feb6ff93010e722cd93aec8655b6 to your computer and use it in GitHub Desktop.
Save SchofieChen/c311feb6ff93010e722cd93aec8655b6 to your computer and use it in GitHub Desktop.
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 __future__ import division
n = 100
x = rand(n)
X = fft(x)
# total energy
sum(abs(x)**2) # 33.803768126149471, for instance
sum(abs(X)**2)/n # 33.803768126149457, for instance
# root mean square
rms_flat(x) # 0.56180320113171067, for instance
rms_flat(X)/sqrt(n) # 0.56180320113171078, for instance
rms_flat(abs(X))/sqrt(n) # 0.56180320113171078, for instance (faster if you've already computed abs(X))
# -*- coding: utf-8 -*-
from __future__ import division
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