1
$\begingroup$

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
Improve this question
$\endgroup$
1

2 Answers 2

Reset to default
1
$\begingroup$

There are several things going on. First of all, the real FFT function (rfft) returns data in a "packed" format, which means that you don't get the actual FFT data, but its shortened convenient representation. Use the regunal fft instead.

Secondly, I don't know where you got the number by which you scale your results on these two lines:

fx = numpy.fft.rfft(sx)/(n/2)
parseval_fx = numpy.sum(numpy.abs(fx)**2)/no_of_points

None of those make sense to me. So let's fix that.

In order to take an orthonormal FFT transform (no scaling of components), we need to scale by the square root of the number of points. So here's what I propose:

# Some python code

import numpy
import math

# 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 1/sqrt(N)
# taking the real fft (rfft) only the positive frequencies are calculated
fx = numpy.fft.fft(sx)/(math.sqrt(len(sx)))

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

print parseval_sx, " equals ", parseval_fx
Improve this answer
$\endgroup$
9
  • $\begingroup$ Your example works fine. However: 1. By means of the scaling factor n/2 (n is number of time-domain samples) the absolute of the value of the 100th FFT bin (numpy.abs(fx[100])) is exactly the amplitude of the time-domain sinus. This is now no longer valid. 2. Scaling factor no_of_points is due to Parselval's relation for DFT, see last formular on Parseval's theorem. Why is it missing in your code but nevertheless well working? $\endgroup$
    – Frank
    May 23, 2014 at 10:05
  • $\begingroup$ According to Numpy rfft Doc rfft actually outputs the positive frequencies as complex numbers rather than in a "packed" representation. However, when using numpy.fft.rfft in your example one needs to applay a factor of 2 on parseval_fx in order to obtain the correct value. However, why don't you need to use the right scaling of the FFT amplitudes in order to fulfuill parceval? $\endgroup$
    – Frank
    May 23, 2014 at 10:47
  • $\begingroup$ SciPy's rfft does the packed representation, while NumPy's rfft does the complex representation, which is easier to apply to N-dimensional arrays $\endgroup$
    – endolith
    May 23, 2014 at 13:41
  • 2
    $\begingroup$ @Frank Applying a value of 2 is not the right solution, since the DC frequency and Nyquist frequency bins are only listed once. You may just be careful and multiply the correct bins by 2 after squaring, but Parseval's theorem is a mathematical statement. And code that would get closest to the statement of the theorem, in my opinion, is the one above. $\endgroup$
    – Phonon
    May 23, 2014 at 18:55
  • 1
    $\begingroup$ @endolith The mean operator will not change the outcome (you're multiplying both sides of the equation by a constant), but Parseval's theorem does not deal with means. It's simply a statement about orthonormality of the Fourier transform. $\endgroup$
    – Phonon
    May 23, 2014 at 18:57
3
$\begingroup$

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.