From: Fago, M. - A. <Mat...@it...> - 2008-12-01 16:03:42
|
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: mat...@li... 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. |