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(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
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>
matplotlib-users List Signup and Options
------------------------------------------------------------------------------
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
matplotlib-users List Signup and Options