Re: [Matplotlib-users] Possible error in psd, specgram and company

[Eric L. <lar...@gm...>]

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
>
>