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