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:

   
http://sourceforge.net/mailarchive/message.php?msg_id=472101A6.9080206%40isla.hawaii.edu

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:

    http://www.mathworks.com/access/helpdesk/help/pdf_doc/signal/signal_tb.pdf

(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[0]*sin(2*pi*f[0]*t) + A[1]*sin(2*pi*f[1]*t) + 5*randn(len(t))
    Pxx,freqs = psd(xn,Fs=fs,NFFT=150,noverlap=75,pad_to=512)

    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.

-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

Reply via email to