Menu
▾
▴
Re: [Matplotlib-users] PSD amplitudes
|
From: Joseph P. <jp...@is...> - 2007-10-26 01:52:10
|
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta content="text/html;charset=ISO-8859-1" http-equiv="Content-Type">
</head>
<body bgcolor="#ffffff" text="#000066">
<font face="Times New Roman">spectral density is by convention a 1Hz
binwidth, not an arbitrary one, units of A^2/Hz.<br>
<br>
perhaps if you manually compute the spectral density of a sine wave,
you will easily see<br>
that they don't have infinite power, R is the autocorrelation of the
Asin(wt):<br>
</font><img src="cid:par...@is..." alt=""><br>
<font face="Times New Roman"><br>
Back to the original question:<br>
<br>
Is there evidence that the matplotlib PSD spectral amplitudes are
accurate?<br>
say by comparison with Matlab results, or a synthetic signal as in the
example, or <br>
from considerations of basic DSP as in the references?<br>
<br>
</font><br>
<a class="moz-txt-link-abbreviated" href="mailto:bre...@un...">bre...@un...</a> wrote:
<blockquote
cite="mid:OFE...@un..."
type="cite"><br>
<font face="sans-serif" size="2">There is certainly differences
(usually
of a factor of PI) in the various definitions used for PSDs, but a
simple
sign wave has an infinite power density at the sine wave frequency.
Are
we agreed on that?</font>
<br>
<br>
<font face="sans-serif" size="2">Use of windowing will modify this
comment
somewhat (so it probably won't really go to infinity) but the basic
fact
remains. The units of a PSD are amp^2/Hz. The MS of a signal
between two frequencies should equal the area under the PSD between
those
frequencies (with allowance for different definitions/factors of PI).
As
I said, for a sign wave the frequency band can be made arbitrarily
small
about the sine wave frequency, but the power between these bands
remains
constant. Therefore the PSD goes to infinity. Otherwise it
isn't a density.</font>
<br>
<br>
<br>
<br>
<br>
<table width="100%">
<tbody>
<tr valign="top">
<td width="40%"><font face="sans-serif" size="1"><b>Joseph Park
<a class="moz-txt-link-rfc2396E" href="mailto:jp...@is..."><jp...@is...></a></b>
</font><br>
<font face="sans-serif" size="1">Sent by:
<a class="moz-txt-link-abbreviated" href="mailto:mat...@li...">mat...@li...</a></font>
<p><font face="sans-serif" size="1">26/10/2007 10:49 AM</font>
</p>
</td>
<td width="59%">
<table width="100%">
<tbody>
<tr valign="top">
<td>
<div align="right"><font face="sans-serif" size="1">To</font></div>
</td>
<td><br>
</td>
</tr>
<tr valign="top">
<td>
<div align="right"><font face="sans-serif" size="1">cc</font></div>
</td>
<td><font face="sans-serif" size="1"><a class="moz-txt-link-abbreviated" href="mailto:mat...@li...">mat...@li...</a></font>
</td>
</tr>
<tr valign="top">
<td>
<div align="right"><font face="sans-serif" size="1">Subject</font></div>
</td>
<td><font face="sans-serif" size="1">Re:
[Matplotlib-users] PSD amplitudes</font></td>
</tr>
</tbody>
</table>
<br>
<table>
<tbody>
<tr valign="top">
<td>
<br>
</td>
<td><br>
</td>
</tr>
</tbody>
</table>
<br>
</td>
</tr>
</tbody>
</table>
<br>
<br>
<br>
<font color="#000066" face="Times New Roman" size="3">is the
suggestion
that the matplotlib algorithm is correct in computing PSD amplitudes?<br>
<br>
btw, increasing nFFT increases the number of points used in the FFT,
which
<br>
increases the spectral <b>frequency</b> resolution (smaller binwidth)
but
for a limited data set<br>
of N points, as is the case in the example, decreases the number of
data
averages <br>
thereby decreasing the spectral <b>amplitude</b> resolution
(accuracy).
keep in mind that<br>
just changing nFFT without making a corresponding change in overlap
will
oversample<br>
the data, thereby skewing the amplitudes. <br>
<br>
in any case, the amplitude change is not approaching infinity, even if
you set nFFT to<br>
6000, which is the length of the timeseries, the amplitudes are ~35dB,
adjust variable ymax<br>
to see this. <br>
<br>
to review issues of spectral/amplitude resolution, windowing/overlap,
etc,
a good<br>
reference is Random Data by Bendat &Piersol:</font><font
color="blue" face="Times New Roman" size="3"><u><br>
</u></font><a moz-do-not-send="true"
href="http://www.amazon.com/Random-Data-Analysis-Measurement-Procedures/dp/0471317330"><font
color="blue" face="Times New Roman" size="3"><u>http://www.amazon.com/Random-Data-Analysis-Measurement-Procedures/dp/0471317330</u></font></a><font
color="#000066" face="Times New Roman" size="3"><br>
<br>
i remain unconvinced that the PSD amplitudes are reasonable, which only
leaves Matlab<br>
as an alternative... that's a hard pill to swallow... matplotlib is
clearly
preferable.<br>
</font><font color="#000066" size="3"><br>
</font><font color="blue" size="3"><u><br>
</u></font><a moz-do-not-send="true"
href="mailto:bre...@un..."><font color="blue"
size="3"><u>bre...@un...</u></font></a><font
color="#000066" size="3">
wrote: </font>
<br>
<font color="#000066" face="sans-serif" size="2"><br>
If you lower the resolution (ie increase nFFT) in your program you will
see that the PSD does indeed increase. I think it may be on the way
to infinity.</font><font color="#000066" size="3"> <br>
<br>
<br>
<br>
</font>
<table width="100%">
<tbody>
<tr valign="top">
<td width="54%"><font color="#000066" face="sans-serif" size="1"><b>Joseph
Park
</b></font><a moz-do-not-send="true"
href="mailto:jp...@is..."><font color="blue"
face="sans-serif" size="1"><b><u><jp...@is...></u></b></font></a><font
color="#000066" face="sans-serif" size="1">
<br>
Sent by: </font><a moz-do-not-send="true"
href="mailto:mat...@li..."><font
color="blue" face="sans-serif" size="1"><u>mat...@li...</u></font></a><font
color="#000066" size="3">
</font>
<p><font color="#000066" face="sans-serif" size="1">26/10/2007
10:05 AM</font><font color="#000066" size="3">
</font></p>
</td>
<td width="45%"><br>
<table width="100%">
<tbody>
<tr valign="top">
<td width="17%">
<div align="right"><font color="#000066" face="sans-serif"
size="1">To</font></div>
</td>
<td width="82%"><a moz-do-not-send="true"
href="mailto:mat...@li..."><font color="blue"
face="sans-serif" size="1"><u>mat...@li...</u></font></a><font
color="#000066" size="3">
</font></td>
</tr>
<tr valign="top">
<td>
<div align="right"><font color="#000066" face="sans-serif"
size="1">cc</font></div>
</td>
<td><br>
</td>
</tr>
<tr valign="top">
<td>
<div align="right"><font color="#000066" face="sans-serif"
size="1">Subject</font></div>
</td>
<td><font color="#000066" face="sans-serif" size="1">Re:
[Matplotlib-users]
PSD amplitudes</font></td>
</tr>
</tbody>
</table>
<br>
<br>
<table width="100%">
<tbody>
<tr valign="top">
<td width="50%">
<br>
</td>
<td width="50%"><br>
</td>
</tr>
</tbody>
</table>
<br>
</td>
</tr>
</tbody>
</table>
<br>
<font color="#000066" size="3"><br>
<br>
</font><font color="#000066" face="sans-serif" size="2"><br>
Shouldn't the PSD for a simple sine wave tend to infinity</font><font
color="#000066" face="Courier New" size="2"><br>
<br>
the spectral resolution will impact the amplitude, if you<br>
are <b>not</b> dealing with a density. by definition a spectral density<br>
has applied the bandwidth resolution correction. the PSD amplitude<br>
should correspond to the RMS amplitude of the sine wave. in the<br>
example a 1VRMS amplitude sine wave (time domain) should have a<br>
PSD power of 20*log(1V) = 0dB. The windowing function will impact<br>
this ideal number a bit, but certainly not by 25dB.</font><font
color="blue" size="3"><u><br>
<br>
</u></font><a moz-do-not-send="true"
href="mailto:bre...@un..."><font color="blue"
size="3"><u>bre...@un...</u></font></a><font
color="#000066" size="3">
wrote: </font><font color="#000066" face="sans-serif" size="2"><br>
<br>
Are you sure that the answer should be zero? Shouldn't the PSD for
a simple sine wave tend to infinity (depending on the resolution)?</font><font
color="#000066" size="3">
<br>
<br>
<br>
</font>
<table width="100%">
<tbody>
<tr valign="top">
<td width="54%"><font color="#000066" face="sans-serif" size="1"><b>Joseph
Park
</b></font><a moz-do-not-send="true"
href="mailto:jp...@is..."><font color="blue"
face="sans-serif" size="1"><b><u><jp...@is...></u></b></font></a><font
color="#000066" face="sans-serif" size="1">
<br>
Sent by: </font><a moz-do-not-send="true"
href="mailto:mat...@li..."><font
color="blue" face="sans-serif" size="1"><u>mat...@li...</u></font></a><font
color="#000066" size="3">
</font>
<p><font color="#000066" face="sans-serif" size="1">26/10/2007
06:50 AM</font><font color="#000066" size="3">
</font></p>
</td>
<td width="45%"><br>
<table width="100%">
<tbody>
<tr valign="top">
<td width="17%">
<div align="right"><font color="#000066" face="sans-serif"
size="1">To</font></div>
</td>
<td width="82%"><a moz-do-not-send="true"
href="mailto:mat...@li..."><font color="blue"
face="sans-serif" size="1"><u>mat...@li...</u></font></a><font
color="#000066" size="3">
</font></td>
</tr>
<tr valign="top">
<td>
<div align="right"><font color="#000066" face="sans-serif"
size="1">cc</font></div>
</td>
<td><br>
</td>
</tr>
<tr valign="top">
<td>
<div align="right"><font color="#000066" face="sans-serif"
size="1">Subject</font></div>
</td>
<td><font color="#000066" face="sans-serif" size="1">[Matplotlib-users]
PSD
amplitudes</font></td>
</tr>
</tbody>
</table>
<br>
<font color="#000066" size="3"><br>
</font>
<br>
<table width="100%">
<tbody>
<tr valign="top">
<td width="50%">
<br>
</td>
<td width="50%"><br>
</td>
</tr>
</tbody>
</table>
<br>
</td>
</tr>
</tbody>
</table>
<br>
<font color="#000066" size="3"><br>
<br>
<br>
<br>
Please try the attached script.<br>
The answer should be ~0 dB for each of the frequencies.<br>
Most likely a simple scaling issue/parameter of which i'm ignorant.</font><font
color="#000066" size="3"><tt><br>
<br>
-- </tt></font><font color="#000066" size="3"><br>
<br>
<br>
______________________________________________________________________<br>
This email has been scanned by the MessageLabs Email Security System.<br>
For more information please visit </font><a moz-do-not-send="true"
href="http://www.messagelabs.com/email"><font color="blue" size="3"><u>http://www.messagelabs.com/email</u></font></a><font
color="#000066" size="3">
<br>
______________________________________________________________________</font><font
color="#000066" size="2"><tt>##----------------------------------------------------------------------------<br>
## Name: psd_scale.py<br>
## <br>
## Purpose: Test Power Spectral Density of 1Vrms data<br>
## Depends on Python SciPy and NumPy<br>
## <br>
## Author: J Park<br>
##<br>
## Created: 10/17/07<br>
##<br>
## Modified: <br>
##----------------------------------------------------------------------------<br>
<br>
try:<br>
from numpy import * # </tt></font><a moz-do-not-send="true"
href="http://www.numpy.org/"><font color="blue" size="2"><tt><u>www.numpy.org</u></tt></font></a><font
color="#000066" size="2"><tt>
numpy.scipy.org<br>
except ImportError: <br>
print "Failed to import numpy."<br>
<br>
try:<br>
import pylab as mp # matplotlib.sourceforge.net<br>
from matplotlib.font_manager import fontManager, FontProperties
<br>
except ImportError: <br>
print "Failed to import pylab."<br>
<br>
<br>
# Default Parameters<br>
nFFT = 1024 <br>
overlap = 512 <br>
freqSample = 100. <br>
PlotAll = False<br>
WriteOutput = False<br>
<br>
##----------------------------------------------------------------------------<br>
## Main module<br>
def main():<br>
<br>
deltaF = freqSample/nFFT # Frequency resolution in Hz<br>
deltaT = 1./freqSample # Sample interval<br>
print 'Sample interval %e (s)' % (deltaT)<br>
print 'Frequency resolution %e (Hz)' % (deltaF)<br>
<br>
# Setup Plots<br>
#
----------------------------------------------------------------------<br>
mp.figure(1)<br>
mp.title ( "PSD" )<br>
mp.ylabel( "(dB)" )<br>
mp.xlabel( "Frequency (Hz)" )<br>
legendFont = FontProperties(size='small')<br>
<br>
ymin = 0<br>
ymax = 30<br>
xmin = 0<br>
xmax = 50<br>
xticks = 5<br>
yticks = 5<br>
<br>
if PlotAll:<br>
mp.figure(2)<br>
mp.title ( "Input Timeseries" )<br>
mp.ylabel( "Amplitude" )<br>
mp.xlabel( "time (s)" )<br>
<br>
# Create some synthetic data with unity RMS amplitude = 0 dB<br>
#
----------------------------------------------------------------------<br>
t = mp.arange(0., 60., deltaT) # 60 seconds at deltaT interval<br>
A = 1.414<br>
<br>
y0 = A * sin( 2. * math.pi * 5 * t )<br>
y1 = A * sin( 2. * math.pi * 10 * t )<br>
y2 = A * sin( 2. * math.pi * 20 * t )<br>
y3 = A * sin( 2. * math.pi * 30 * t )<br>
y4 = A * sin( 2. * math.pi * 40 * t )<br>
y5 = A * sin( 2. * math.pi * 45 * t )<br>
<br>
dataList = [ y0, y1, y2, y3, y4, y5 ]<br>
<br>
for data in dataList:<br>
inputDataLen = len( data )<br>
numAverages = math.floor( inputDataLen / (overlap)
) - 1<br>
normalizedRandomError = 1./math.sqrt( numAverages
)<br>
print "%d points" % ( inputDataLen ),<br>
print "%d averages" % (numAverages),<br>
print "normalized random error %.3f" % (
normalizedRandomError )<br>
<br>
mp.figure(1)<br>
(Pxx, freqs) = mp.psd( data,<br>
NFFT = nFFT,<br>
Fs = freqSample,<br>
noverlap = overlap,<br>
lw = 2,<br>
label = '' )<br>
<br>
Pxx_dB = 10.*log10(Pxx)<br>
<br>
if PlotAll:<br>
mp.figure(2)<br>
mp.plot(t, data, label='' )<br>
<br>
# Write Output data<br>
#
----------------------------------------------------------------------<br>
if WriteOutput:<br>
PxxLen = len(Pxx)<br>
OutputFile = "PSD.dat"<br>
fdOutFile = open( OutputFile, 'a' )<br>
fdOutFile.write( "Freq\t\tPower(dB)\n"
)<br>
for i in range(PxxLen):<br>
fdOutFile.write( "%.4e\t%.3f\n"
% ( freqs[i], Pxx_dB[i] ) )<br>
fdOutFile.close()<br>
print "Wrote ", PxxLen, "
points to ", OutputFile<br>
<br>
<br>
# Show the Plot<br>
#
----------------------------------------------------------------------<br>
mp.figure(1)<br>
mp.axis([xmin, xmax, ymin, ymax])<br>
mp.xticks( arange(xmin, xmax+1, xticks) )<br>
mp.yticks( arange(ymin, ymax , yticks) )<br>
mp.title('')<br>
mp.xlabel('Frequency (Hz)')<br>
mp.ylabel(r'$\tt{dB re V^2/Hz}$')<br>
#mp.legend( loc='upper right', prop=legendFont )<br>
if WriteOutput:<br>
plotFileName = "PSD.png"<br>
mp.savefig( plotFileName )<br>
print "Wrote png image to ", plotFileName<br>
if PlotAll:<br>
mp.figure(2)<br>
#mp.legend( loc='lower left', prop=legendFont )<br>
mp.show()<br>
<br>
print "Normal Exit"<br>
## Main module<br>
##----------------------------------------------------------------------------<br>
<br>
##----------------------------------------------------------------------------<br>
## Provide for cmd line invocation<br>
if __name__ == "__main__":<br>
main()<br>
<br>
-------------------------------------------------------------------------<br>
This SF.net email is sponsored by: Splunk Inc.<br>
Still grepping through log files to find problems? Stop.<br>
Now Search log events and configuration files using AJAX and a browser.<br>
Download your FREE copy of Splunk now >> </tt></font><a
moz-do-not-send="true"
href="http://get.splunk.com/_______________________________________________"><font
color="blue" size="2"><tt><u>http://get.splunk.com/_______________________________________________</u></tt></font></a><font
color="#000066" size="2"><tt><br>
Matplotlib-users mailing list</tt></font><font color="blue" size="3"><u><br>
</u></font><a moz-do-not-send="true"
href="mailto:Mat...@li..."><font color="blue"
size="2"><tt><u>Mat...@li...</u></tt></font></a><font
color="blue" size="3"><u><br>
</u></font><a moz-do-not-send="true"
href="https://lists.sourceforge.net/lists/listinfo/matplotlib-users"><font
color="blue" size="2"><tt><u>https://lists.sourceforge.net/lists/listinfo/matplotlib-users</u></tt></font></a><font
color="#000066" size="3"><br>
<br>
<br>
UNITED GROUP<br>
This email message is the property of United Group. The information in
this email is confidential and may be legally privileged. It is
intended
solely for the addressee. Access to this email by anyone else is
unauthorised.
If you are not the intended recipient, you may not disclose, copy or
distribute
this email, nor take or omit to take any action in reliance on it.
United
Group accepts no liability for any damage caused by this email or any
attachments
due to viruses, interference, interception, corruption or unauthorised
access.<br>
If you have received this email in error, please notify United Group
immediately
by email to the sender's email address and delete this document. <br>
</font><font color="#000066" size="3"><tt><br>
-- </tt></font><font color="#000066" size="3"><br>
<br>
<br>
______________________________________________________________________<br>
This email has been scanned by the MessageLabs Email Security System.<br>
For more information please visit </font><a moz-do-not-send="true"
href="http://www.messagelabs.com/email"><font color="blue" size="3"><u>http://www.messagelabs.com/email</u></font></a><font
color="#000066" size="3">
<br>
______________________________________________________________________</font><font
color="#000066" size="2"><tt>-------------------------------------------------------------------------<br>
This SF.net email is sponsored by: Splunk Inc.<br>
Still grepping through log files to find problems? Stop.<br>
Now Search log events and configuration files using AJAX and a browser.<br>
Download your FREE copy of Splunk now >> </tt></font><a
moz-do-not-send="true"
href="http://get.splunk.com/_______________________________________________"><font
color="blue" size="2"><tt><u>http://get.splunk.com/_______________________________________________</u></tt></font></a><font
color="#000066" size="2"><tt><br>
Matplotlib-users mailing list</tt></font><font color="blue" size="2"><tt><u><br>
</u></tt></font><a moz-do-not-send="true"
href="mailto:Mat...@li..."><font color="blue"
size="2"><tt><u>Mat...@li...</u></tt></font></a><font
color="blue" size="2"><tt><u><br>
</u></tt></font><a moz-do-not-send="true"
href="https://lists.sourceforge.net/lists/listinfo/matplotlib-users"><font
color="blue" size="2"><tt><u>https://lists.sourceforge.net/lists/listinfo/matplotlib-users</u></tt></font></a><font
color="#000066" size="3"><br>
<br>
<br>
UNITED GROUP<br>
This email message is the property of United Group. The information in
this email is confidential and may be legally privileged. It is
intended
solely for the addressee. Access to this email by anyone else is
unauthorised.
If you are not the intended recipient, you may not disclose, copy or
distribute
this email, nor take or omit to take any action in reliance on it.
United
Group accepts no liability for any damage caused by this email or any
attachments
due to viruses, interference, interception, corruption or unauthorised
access.<br>
If you have received this email in error, please notify United Group
immediately
by email to the sender's email address and delete this document.</font>
<br>
<br>
<font color="#000066" size="3"><tt>-- <br>
</tt></font>
<br>
<font color="#000066" size="3"><br>
______________________________________________________________________<br>
This email has been scanned by the MessageLabs Email Security System.<br>
For more information please visit <a class="moz-txt-link-freetext" href="http://www.messagelabs.com/email">http://www.messagelabs.com/email</a> <br>
______________________________________________________________________</font><font
size="2"><tt>-------------------------------------------------------------------------<br>
This SF.net email is sponsored by: Splunk Inc.<br>
Still grepping through log files to find problems? Stop.<br>
Now Search log events and configuration files using AJAX and a browser.<br>
Download your FREE copy of Splunk now >>
<a class="moz-txt-link-freetext" href="http://get.splunk.com/_______________________________________________">http://get.splunk.com/_______________________________________________</a><br>
Matplotlib-users mailing list<br>
<a class="moz-txt-link-abbreviated" href="mailto:Mat...@li...">Mat...@li...</a><br>
<a class="moz-txt-link-freetext" href="https://lists.sourceforge.net/lists/listinfo/matplotlib-users">https://lists.sourceforge.net/lists/listinfo/matplotlib-users</a><br>
</tt></font>
<br>
<br>
UNITED GROUP<br>
This email message is the property of United Group. The information in
this email is confidential and may be legally privileged. It is
intended solely for the addressee. Access to this email by anyone else
is unauthorised. If you are not the intended recipient, you may not
disclose, copy or distribute this email, nor take or omit to take any
action in reliance on it. United Group accepts no liability for any
damage caused by this email or any attachments due to viruses,
interference, interception, corruption or unauthorised access.<br>
If you have received this email in error, please notify United Group
immediately by email to the sender's email address and delete this
document.<br>
</blockquote>
<br>
<pre class="moz-signature" cols="80">--
</pre>
</body>
</html>
|