Thread: [Matplotlib-users] 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 <
lui...@gm...> 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
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>
>

On 2013/11/21 8:45 AM, Eric Larson wrote:
> 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

>
> Eric
>
>
>
> On Thu, Nov 21, 2013 at 10:30 AM, Luis Miguel García-Cuevas González
> <lui...@gm... <mailto:lui...@gm...>> 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
>     Mat...@li...
>     <mailto:Mat...@li...>
>     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
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>