Possible error in psd, specgram and company

Hello!

I have a doubt about the way matplotlib computes (and plots) a psd. I have
the following example code (I've imported pylab for convenience here):

  from pylab import *

  A = 10. ** (10. / 20.) / sqrt(2.) * 2. # This amplitude should give 10 dB
  frequency = 5.
  sampling_frequency = 1000.
  t_max = 12.
  t = arange(0., t_max, 1. / sampling_frequency)
  x = A * cos(2. * pi * frequency * t)

  Pxx, freqs = mlab.psd(x, Fs=sampling_frequency, NFFT=len(x),
      window=mlab.window_none)

  dB_actual = 20. * log10(rms_flat(x)) # dB_actual = 10., as expected
  dB_psd = 10. * log10(Pxx.max()) # dB_psd = 20.791812460476248

As I understand, this should not happen this way. dB_psd should be equal to
10, as dB_actual. I can get it to compute the correct value with this:

  Pxx, freqs = mlab.psd(x / sqrt(t_max), Fs=sampling_frequency, NFFT=len(x),
      window=mlab.window_none)

(i.e., dividing the variable by the square root of the time). Is this really
the expected behavior or this is a bug?

If it does matter, I'm using using matplotlib 1.3.1.

Thanks for your help!

···

--
Luis Miguel García-Cuevas González

Hey Luis,

I think in general the application of normalizations in DFT / FFT (or even the continuous transforms, for that matter) can vary from implementation to another, in terms of whether they’re applied during the forward transform, inverse transform, or half on each. That sounds like what you’re running into here. MATLAB and other packages (FFTW, and apparently whatever matplotlib uses internally for FFTs) tend to normalize during the inverse transform, so you would need to put the normalization in yourself to see Parseval’s theorem work out, which is essentially what you’re showing.

To my mind, this isn’t a bug, so much as something users should be made aware of so they understand the issues at stake. It would be good to add a note to the PSD computation about these normalization issues, assuming one doesn’t already exist. I have not contributed to matplotlib before, but I’m happy to take a stab at this if people think it would be helpful. My guess is that matplotlib won’t be changed to include the suggested normalization by default (but maybe it could be added as an option?) because this would break backward-compatibility for users that have taken this factor into account in their own code.

Eric

···

On Thu, Nov 21, 2013 at 10:30 AM, Luis Miguel García-Cuevas González <luismiguelgcg@…287…> wrote:

Hello!

I have a doubt about the way matplotlib computes (and plots) a psd. I have

the following example code (I’ve imported pylab for convenience here):

from pylab import *

A = 10. ** (10. / 20.) / sqrt(2.) * 2. # This amplitude should give 10 dB

frequency = 5.

sampling_frequency = 1000.

t_max = 12.

t = arange(0., t_max, 1. / sampling_frequency)

x = A * cos(2. * pi * frequency * t)

Pxx, freqs = mlab.psd(x, Fs=sampling_frequency, NFFT=len(x),

  window=mlab.window_none)

dB_actual = 20. * log10(rms_flat(x)) # dB_actual = 10., as expected

dB_psd = 10. * log10(Pxx.max()) # dB_psd = 20.791812460476248

As I understand, this should not happen this way. dB_psd should be equal to

10, as dB_actual. I can get it to compute the correct value with this:

Pxx, freqs = mlab.psd(x / sqrt(t_max), Fs=sampling_frequency, NFFT=len(x),

  window=mlab.window_none)

(i.e., dividing the variable by the square root of the time). Is this really

the expected behavior or this is a bug?

If it does matter, I’m using using matplotlib 1.3.1.

Thanks for your help!

Luis Miguel García-Cuevas González


Shape the Mobile Experience: Free Subscription

Software experts and developers: Be at the forefront of tech innovation.

Intel® Software Adrenaline delivers strategic insight and game-changing

conversations that shape the rapidly evolving mobile landscape. Sign up now.

http://pubads.g.doubleclick.net/gampad/clk?id=63431311&iu=/4140/ostg.clktrk


Matplotlib-users mailing list

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users

Hey Luis,

I think in general the application of normalizations in DFT / FFT (or
even the continuous transforms, for that matter) can vary from
implementation to another, in terms of whether they're applied during
the forward transform, inverse transform, or half on each. That sounds
like what you're running into here. MATLAB and other packages (FFTW, and
apparently whatever matplotlib uses internally for FFTs) tend to
normalize during the inverse transform, so you would need to put the
normalization in yourself to see Parseval's theorem work out, which is
essentially what you're showing.

It is not a matter of how the FFT is normalized, but rather of whether one is calculating a power spectrum, or power spectral density (psd).

For example, see:

http://pubman.mpdl.mpg.de/pubman/item/escidoc:152164:1/component/escidoc:152163/395068.pdf

To my mind, this isn't a bug, so much as something users should be made
aware of so they understand the issues at stake. It would be good to add
a note to the PSD computation about these normalization issues, assuming
one doesn't already exist. I have not contributed to matplotlib before,
but I'm happy to take a stab at this if people think it would be
helpful. My guess is that matplotlib won't be changed to include the
suggested normalization by default (but maybe it could be added as an
option?) because this would break backward-compatibility for users that
have taken this factor into account in their own code.

I think the mpl psd is consistent with other psd implementations and with the definition of the psd.

Eric

···

On 2013/11/21 8:45 AM, Eric Larson wrote:

Eric

On Thu, Nov 21, 2013 at 10:30 AM, Luis Miguel Garc�a-Cuevas Gonz�lez > <luismiguelgcg@…287… <mailto:luismiguelgcg@…287…>> wrote:

    Hello!

    I have a doubt about the way matplotlib computes (and plots) a psd.
      I have
    the following example code (I've imported pylab for convenience here):

       from pylab import *

       A = 10. ** (10. / 20.) / sqrt(2.) * 2. # This amplitude should
    give 10 dB
       frequency = 5.
       sampling_frequency = 1000.
       t_max = 12.
       t = arange(0., t_max, 1. / sampling_frequency)
       x = A * cos(2. * pi * frequency * t)

       Pxx, freqs = mlab.psd(x, Fs=sampling_frequency, NFFT=len(x),
           window=mlab.window_none)

       dB_actual = 20. * log10(rms_flat(x)) # dB_actual = 10., as expected
       dB_psd = 10. * log10(Pxx.max()) # dB_psd = 20.791812460476248

    As I understand, this should not happen this way. dB_psd should be
    equal to
    10, as dB_actual. I can get it to compute the correct value with this:

       Pxx, freqs = mlab.psd(x / sqrt(t_max), Fs=sampling_frequency,
    NFFT=len(x),
           window=mlab.window_none)

    (i.e., dividing the variable by the square root of the time). Is
    this really
    the expected behavior or this is a bug?

    If it does matter, I'm using using matplotlib 1.3.1.

    Thanks for your help!

    --
    Luis Miguel Garc�a-Cuevas Gonz�lez

    ------------------------------------------------------------------------------
    Shape the Mobile Experience: Free Subscription
    Software experts and developers: Be at the forefront of tech innovation.
    Intel(R) Software Adrenaline delivers strategic insight and
    game-changing
    conversations that shape the rapidly evolving mobile landscape. Sign
    up now.
    http://pubads.g.doubleclick.net/gampad/clk?id=63431311&iu=/4140/ostg.clktrk
    _______________________________________________
    Matplotlib-users mailing list
    Matplotlib-users@lists.sourceforge.net
    <mailto:Matplotlib-users@lists.sourceforge.net>
    https://lists.sourceforge.net/lists/listinfo/matplotlib-users

------------------------------------------------------------------------------
Shape the Mobile Experience: Free Subscription
Software experts and developers: Be at the forefront of tech innovation.
Intel(R) Software Adrenaline delivers strategic insight and game-changing
conversations that shape the rapidly evolving mobile landscape. Sign up now.
http://pubads.g.doubleclick.net/gampad/clk?id=63431311&iu=/4140/ostg.clktrk

_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users