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: > > > 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?
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 Graduate Research Assistant School of Meteorology University of Oklahoma ------------------------------------------------------------------------- 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
