Screenshot instructions:
Windows
Mac
Red Hat Linux
Ubuntu
Click URL instructions:
Rightclick on ad, choose "Copy Link", then paste here →
(This may not be possible with some types of ads)
From: Ryan May <rmay31@gm...>  20081125 20:15:22

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 223 and 224). 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 224). > > 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 