 # Matplotlib PSD bug?

[I'm not sure if this is best in 'devel' or 'users']

I'm trying to compute PSDs using matplotlib.mlab.psd and came across the "PSD amplitudes" thread from last year:

Using the latest version of psd on svn trunk (rev 6429) that added support for pad_to, I can now compute the Matlab pwelch
example using matplotlib. This example is given in the Signal Processing Toolbox User's Guide:

(look on pages 2-23 and 2-24). Note I do not have easy access to Matlab itself, so I'm just using this published example.

The Matlab code is as follows:

randn('state',1)
fs = 1000; % Sampling frequency
t = (0:0.3*fs)./fs; % 301 samples
A = [2 8]; % Sinusoid amplitudes (row vector)
f = [150;140]; % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 5*randn(size(t));
Hs = spectrum.welch('rectangular',150,50);
psd(Hs,xn,'Fs',fs,'NFFT',512);

This produces a fairly noisy signal from -20 to -10 dB, with a strong peak of ~6 dB at 150 Hz (see the plot on page 2-24).

While my equivalent (?) python code is:

from scipy import *
from mlabsvn import psd # This is a local copy of svn revision 6429 of matplotlib.mlab
from pylab import *
fs=1000
t=linspace(0,0.3,0.3*fs+1)
A=[2,8]
f=[150,140]
xn=A*sin(2*pi*f*t) + A*sin(2*pi*f*t) + 5*randn(len(t))

figure()
plot(freqs, 10*log10(Pxx) )
show()

However, this produces a peak of over 30 dB at 150 Hz. Unless there is a mistake in my code above, there seems to be a
significant difference between the matplotlib and matlab implementations.

I noticed that the values 10*log10(Pxx/len(xn)) produces results that match much better. Without looking more closely at the
code for psd and reviewing Bendat and Piersol I cannot be sure that this is the correct fix.

Does anyone else have any insight? When is the next release planned, and how likely is a fix?

Thanks,
Matt

This e-mail and any files transmitted with it may be proprietary and are intended solely for the use of the individual or entity to whom they are addressed. If you have received this e-mail in error please notify the sender.
Please note that any views or opinions presented in this e-mail are solely those of the author and do not necessarily represent those of ITT Corporation. The recipient should check this e-mail and any attachments for the presence of viruses. ITT accepts no liability for any damage caused by any virus transmitted by this e-mail.

Fago, Matt - AES wrote:

[I'm not sure if this is best in 'devel' or 'users']

I'm trying to compute PSDs using matplotlib.mlab.psd and came across the "PSD amplitudes" thread from last year:

Using the latest version of psd on svn trunk (rev 6429) that added support for pad_to, I can now compute the Matlab pwelch
example using matplotlib. This example is given in the Signal Processing Toolbox User's Guide:

(look on pages 2-23 and 2-24). Note I do not have easy access to Matlab itself, so I'm just using this published example.

The Matlab code is as follows:

randn('state',1)
fs = 1000; % Sampling frequency
t = (0:0.3*fs)./fs; % 301 samples
A = [2 8]; % Sinusoid amplitudes (row vector)
f = [150;140]; % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 5*randn(size(t));
Hs = spectrum.welch('rectangular',150,50);
psd(Hs,xn,'Fs',fs,'NFFT',512);

This produces a fairly noisy signal from -20 to -10 dB, with a strong peak of ~6 dB at 150 Hz (see the plot on page 2-24).

While my equivalent (?) python code is:

from scipy import *
from mlabsvn import psd # This is a local copy of svn revision 6429 of matplotlib.mlab
from pylab import *
fs=1000
t=linspace(0,0.3,0.3*fs+1)
A=[2,8]
f=[150,140]
xn=A*sin(2*pi*f*t) + A*sin(2*pi*f*t) + 5*randn(len(t))

figure()
plot(freqs, 10*log10(Pxx) )
show()

However, this produces a peak of over 30 dB at 150 Hz. Unless there is a mistake in my code above, there seems to be a
significant difference between the matplotlib and matlab implementations.

I noticed that the values 10*log10(Pxx/len(xn)) produces results that match much better. Without looking more closely at the
code for psd and reviewing Bendat and Piersol I cannot be sure that this is the correct fix.

Does anyone else have any insight? When is the next release planned, and how likely is a fix?

I don't have any insight yet, but since I'm the guy who just tweaked it, I guess I'm on the hook to fix it. I'll try to take a look this afternoon.

Ryan

···

--
Ryan May
School of Meteorology
University of Oklahoma

I found bug number 1859027 and have appended the below to the bug report.

When is the next release due and how likely is this to get fixed? I might have time myself
to help in a week or so, but would appreciate some help if someone else has time too (who
has looked at the source before...)

Thanks,
Matt

···

________________________________________
From: Fago, Matt - AES
Sent: Tuesday, November 25, 2008 1:04 PM
To: matplotlib-users@lists.sourceforge.net
Subject: Matplotlib PSD bug?

[I'm not sure if this is best in 'devel' or 'users']

I'm trying to compute PSDs using matplotlib.mlab.psd and came across the "PSD amplitudes" thread from last year:

Using the latest version of psd on svn trunk (rev 6429) that added support for pad_to, I can now compute the Matlab pwelch
example using matplotlib. This example is given in the Signal Processing Toolbox User's Guide:

(look on pages 2-23 and 2-24). Note I do not have easy access to Matlab itself, so I'm just using this published example.

The Matlab code is as follows:

randn('state',1)
fs = 1000; % Sampling frequency
t = (0:0.3*fs)./fs; % 301 samples
A = [2 8]; % Sinusoid amplitudes (row vector)
f = [150;140]; % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 5*randn(size(t));
Hs = spectrum.welch('rectangular',150,50);
psd(Hs,xn,'Fs',fs,'NFFT',512);

This produces a fairly noisy signal from -20 to -10 dB, with a strong peak of ~6 dB at 150 Hz (see the plot on page 2-24).

While my equivalent (?) python code is:

from scipy import *
from mlabsvn import psd # This is a local copy of svn revision 6429 of matplotlib.mlab
from pylab import *
fs=1000
t=linspace(0,0.3,0.3*fs+1)
A=[2,8]
f=[150,140]
xn=A*sin(2*pi*f*t) + A*sin(2*pi*f*t) + 5*randn(len(t))

figure()
plot(freqs, 10*log10(Pxx) )
show()

However, this produces a peak of over 30 dB at 150 Hz. Unless there is a mistake in my code above, there seems to be a
significant difference between the matplotlib and matlab implementations.

I noticed that the values 10*log10(Pxx/len(xn)) produces results that match much better. Without looking more closely at the
code for psd and reviewing Bendat and Piersol I cannot be sure that this is the correct fix.

Does anyone else have any insight? When is the next release planned, and how likely is a fix?

Thanks,
Matt

This e-mail and any files transmitted with it may be proprietary and are intended solely for the use of the individual or entity to whom they are addressed. If you have received this e-mail in error please notify the sender.
Please note that any views or opinions presented in this e-mail are solely those of the author and do not necessarily represent those of ITT Corporation. The recipient should check this e-mail and any attachments for the presence of viruses. ITT accepts no liability for any damage caused by any virus transmitted by this e-mail.

Fago, Matt - AES wrote:

I found bug number 1859027 and have appended the below to the bug report.

When is the next release due and how likely is this to get fixed? I might have time myself
to help in a week or so, but would appreciate some help if someone else has time too (who
has looked at the source before...)

I'm sitting down to look at this right now. Can you tell me if you have this problem with 0.98.3 as well?

Thanks,

Ryan

Ryan May
School of Meteorology
University of Oklahoma

Ryan May wrote:

Fago, Matt - AES wrote:

I found bug number 1859027 and have appended the below to the bug report.

When is the next release due and how likely is this to get fixed? I might have time myself
to help in a week or so, but would appreciate some help if someone else has time too (who
has looked at the source before...)

I'm sitting down to look at this right now. Can you tell me if you have this problem with 0.98.3 as well?

Nevermind, I just realized you can't really do the example without the new support for pad_to. I'm still looking...

Ryan

···

--
Ryan May