[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.
I'll try to take a look this afternoon.