Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

## matplotlib-users

 Re: [Matplotlib-users] PSD amplitudes From: Joseph Park - 2007-10-26 00:07:30 ``` Shouldn't the PSD for a simple sine wave tend to infinity

the spectral resolution will impact the amplitude, if you
are not dealing with a density. by definition a spectral density
has applied the bandwidth resolution correction. the PSD amplitude
should correspond to the RMS amplitude of the sine wave. in the
example a 1VRMS amplitude sine wave (time domain) should have a
PSD power of 20*log(1V) = 0dB. The windowing function will impact
this ideal number a bit, but certainly not by 25dB.

brett.mcsweeney@... wrote:

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)?

Joseph Park <jpark@...>
Sent by: matplotlib-users-bounces@...

26/10/2007 06:50 AM

To
matplotlib-users@...
cc

Subject
[Matplotlib-users] PSD amplitudes

Please try the attached script.
The answer should be ~0 dB for each of the frequencies.
Most likely a simple scaling issue/parameter of which i'm ignorant.

--

______________________________________________________________________
This email has been scanned by the MessageLabs Email Security System.
______________________________________________________________________
##----------------------------------------------------------------------------
## Name:     psd_scale.py
##
## Purpose:  Test Power Spectral Density of 1Vrms data
##           Depends on Python SciPy and NumPy
##
## Author:       J Park
##
## Created:      10/17/07
##
## Modified:
##----------------------------------------------------------------------------

try:
from numpy import *  # www.numpy.org; numpy.scipy.org
except ImportError:
print "Failed to import numpy."

try:
import pylab as mp  # matplotlib.sourceforge.net
from matplotlib.font_manager import fontManager, FontProperties
except ImportError:
print "Failed to import pylab."

# Default Parameters
nFFT          = 1024
overlap       = 512
freqSample    = 100.
PlotAll       = False
WriteOutput   = False

##----------------------------------------------------------------------------
## Main module
def main():

deltaF = freqSample/nFFT # Frequency resolution in Hz
deltaT = 1./freqSample   # Sample interval
print 'Sample interval %e (s)'       % (deltaT)
print 'Frequency resolution %e (Hz)' % (deltaF)

# Setup Plots
# ----------------------------------------------------------------------
mp.figure(1)
mp.title ( "PSD" )
mp.ylabel( "(dB)" )
mp.xlabel( "Frequency (Hz)" )
legendFont = FontProperties(size='small')

ymin = 0
ymax = 30
xmin = 0
xmax = 50
xticks = 5
yticks = 5

if PlotAll:
mp.figure(2)
mp.title ( "Input Timeseries" )
mp.ylabel( "Amplitude" )
mp.xlabel( "time (s)" )

# Create some synthetic data with unity RMS amplitude = 0 dB
# ----------------------------------------------------------------------
t = mp.arange(0., 60., deltaT) # 60 seconds at deltaT interval
A = 1.414

y0 = A * sin( 2. * math.pi * 5  * t )
y1 = A * sin( 2. * math.pi * 10 * t )
y2 = A * sin( 2. * math.pi * 20 * t )
y3 = A * sin( 2. * math.pi * 30 * t )
y4 = A * sin( 2. * math.pi * 40 * t )
y5 = A * sin( 2. * math.pi * 45 * t )

dataList = [ y0, y1, y2, y3, y4, y5 ]

for data in dataList:
inputDataLen = len( data )
numAverages  = math.floor( inputDataLen / (overlap) ) - 1
normalizedRandomError = 1./math.sqrt( numAverages )
print "%d points" % ( inputDataLen ),
print "%d averages" % (numAverages),
print "normalized random error %.3f" % ( normalizedRandomError )

mp.figure(1)
(Pxx, freqs) = mp.psd( data,
NFFT     = nFFT,
Fs       = freqSample,
noverlap = overlap,
lw       = 2,
label    = '' )

Pxx_dB = 10.*log10(Pxx)

if PlotAll:
mp.figure(2)
mp.plot(t, data, label='' )

# Write Output data
# ----------------------------------------------------------------------
if WriteOutput:
PxxLen = len(Pxx)
OutputFile = "PSD.dat"
fdOutFile = open( OutputFile, 'a' )
fdOutFile.write( "Freq\t\tPower(dB)\n" )
for i in range(PxxLen):
fdOutFile.write( "%.4e\t%.3f\n" % ( freqs[i], Pxx_dB[i] ) )
fdOutFile.close()
print "Wrote ", PxxLen, " points to ", OutputFile

# Show the Plot
# ----------------------------------------------------------------------
mp.figure(1)
mp.axis([xmin, xmax, ymin, ymax])
mp.xticks( arange(xmin, xmax+1, xticks) )
mp.yticks( arange(ymin, ymax  , yticks) )
mp.title('')
mp.xlabel('Frequency (Hz)')
mp.ylabel(r'\$\tt{dB re V^2/Hz}\$')
#mp.legend( loc='upper right', prop=legendFont )
if WriteOutput:
plotFileName = "PSD.png"
mp.savefig( plotFileName )
print "Wrote png image to ", plotFileName
if PlotAll:
mp.figure(2)
#mp.legend( loc='lower left', prop=legendFont )
mp.show()

print "Normal Exit"
## Main module
##----------------------------------------------------------------------------

##----------------------------------------------------------------------------
## Provide for cmd line invocation
if __name__ == "__main__":
main()

-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems?  Stop.
Now Search log events and configuration files using AJAX and a browser.
; Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
;

UNITED GROUP
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.
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.

--
```
 Re: [Matplotlib-users] PSD amplitudes From: - 2007-10-26 00:17:43 Attachments: Message as HTML ```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. Joseph Park Sent by: matplotlib-users-bounces@... 26/10/2007 10:05 AM To matplotlib-users@... cc Subject Re: [Matplotlib-users] PSD amplitudes Shouldn't the PSD for a simple sine wave tend to infinity the spectral resolution will impact the amplitude, if you are not dealing with a density. by definition a spectral density has applied the bandwidth resolution correction. the PSD amplitude should correspond to the RMS amplitude of the sine wave. in the example a 1VRMS amplitude sine wave (time domain) should have a PSD power of 20*log(1V) = 0dB. The windowing function will impact this ideal number a bit, but certainly not by 25dB. brett.mcsweeney@... wrote: 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)? Joseph Park Sent by: matplotlib-users-bounces@... 26/10/2007 06:50 AM To matplotlib-users@... cc Subject [Matplotlib-users] PSD amplitudes Please try the attached script. The answer should be ~0 dB for each of the frequencies. Most likely a simple scaling issue/parameter of which i'm ignorant. -- ______________________________________________________________________ This email has been scanned by the MessageLabs Email Security System. For more information please visit http://www.messagelabs.com/email ______________________________________________________________________ ##---------------------------------------------------------------------------- ## Name: psd_scale.py ## ## Purpose: Test Power Spectral Density of 1Vrms data ## Depends on Python SciPy and NumPy ## ## Author: J Park ## ## Created: 10/17/07 ## ## Modified: ##---------------------------------------------------------------------------- try: from numpy import * # http://www.numpy.org numpy.scipy.org except ImportError: print "Failed to import numpy." try: import pylab as mp # matplotlib.sourceforge.net from matplotlib.font_manager import fontManager, FontProperties except ImportError: print "Failed to import pylab." # Default Parameters nFFT = 1024 overlap = 512 freqSample = 100. PlotAll = False WriteOutput = False ##---------------------------------------------------------------------------- ## Main module def main(): deltaF = freqSample/nFFT # Frequency resolution in Hz deltaT = 1./freqSample # Sample interval print 'Sample interval %e (s)' % (deltaT) print 'Frequency resolution %e (Hz)' % (deltaF) # Setup Plots # ---------------------------------------------------------------------- mp.figure(1) mp.title ( "PSD" ) mp.ylabel( "(dB)" ) mp.xlabel( "Frequency (Hz)" ) legendFont = FontProperties(size='small') ymin = 0 ymax = 30 xmin = 0 xmax = 50 xticks = 5 yticks = 5 if PlotAll: mp.figure(2) mp.title ( "Input Timeseries" ) mp.ylabel( "Amplitude" ) mp.xlabel( "time (s)" ) # Create some synthetic data with unity RMS amplitude = 0 dB # ---------------------------------------------------------------------- t = mp.arange(0., 60., deltaT) # 60 seconds at deltaT interval A = 1.414 y0 = A * sin( 2. * math.pi * 5 * t ) y1 = A * sin( 2. * math.pi * 10 * t ) y2 = A * sin( 2. * math.pi * 20 * t ) y3 = A * sin( 2. * math.pi * 30 * t ) y4 = A * sin( 2. * math.pi * 40 * t ) y5 = A * sin( 2. * math.pi * 45 * t ) dataList = [ y0, y1, y2, y3, y4, y5 ] for data in dataList: inputDataLen = len( data ) numAverages = math.floor( inputDataLen / (overlap) ) - 1 normalizedRandomError = 1./math.sqrt( numAverages ) print "%d points" % ( inputDataLen ), print "%d averages" % (numAverages), print "normalized random error %.3f" % ( normalizedRandomError ) mp.figure(1) (Pxx, freqs) = mp.psd( data, NFFT = nFFT, Fs = freqSample, noverlap = overlap, lw = 2, label = '' ) Pxx_dB = 10.*log10(Pxx) if PlotAll: mp.figure(2) mp.plot(t, data, label='' ) # Write Output data # ---------------------------------------------------------------------- if WriteOutput: PxxLen = len(Pxx) OutputFile = "PSD.dat" fdOutFile = open( OutputFile, 'a' ) fdOutFile.write( "Freq\t\tPower(dB)\n" ) for i in range(PxxLen): fdOutFile.write( "%.4e\t%.3f\n" % ( freqs[i], Pxx_dB[i] ) ) fdOutFile.close() print "Wrote ", PxxLen, " points to ", OutputFile # Show the Plot # ---------------------------------------------------------------------- mp.figure(1) mp.axis([xmin, xmax, ymin, ymax]) mp.xticks( arange(xmin, xmax+1, xticks) ) mp.yticks( arange(ymin, ymax , yticks) ) mp.title('') mp.xlabel('Frequency (Hz)') mp.ylabel(r'\$\tt{dB re V^2/Hz}\$') #mp.legend( loc='upper right', prop=legendFont ) if WriteOutput: plotFileName = "PSD.png" mp.savefig( plotFileName ) print "Wrote png image to ", plotFileName if PlotAll: mp.figure(2) #mp.legend( loc='lower left', prop=legendFont ) mp.show() print "Normal Exit" ## Main module ##---------------------------------------------------------------------------- ##---------------------------------------------------------------------------- ## Provide for cmd line invocation if __name__ == "__main__": main() ------------------------------------------------------------------------- This SF.net email is sponsored by: Splunk Inc. Still grepping through log files to find problems? Stop. Now Search log events and configuration files using AJAX and a browser. Download your FREE copy of Splunk now >> http://get.splunk.com/_______________________________________________ Matplotlib-users mailing list Matplotlib-users@... https://lists.sourceforge.net/lists/listinfo/matplotlib-users UNITED GROUP 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. 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. -- ______________________________________________________________________ This email has been scanned by the MessageLabs Email Security System. For more information please visit http://www.messagelabs.com/email ______________________________________________________________________ ------------------------------------------------------------------------- This SF.net email is sponsored by: Splunk Inc. Still grepping through log files to find problems? Stop. Now Search log events and configuration files using AJAX and a browser. Download your FREE copy of Splunk now >> http://get.splunk.com/ _______________________________________________ Matplotlib-users mailing list Matplotlib-users@... https://lists.sourceforge.net/lists/listinfo/matplotlib-users UNITED GROUP 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. 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.```
 Re: [Matplotlib-users] PSD amplitudes From: Joseph Park - 2007-10-26 00:51:24 ``` is the suggestion that the matplotlib algorithm is correct in computing PSD amplitudes?

btw, increasing nFFT increases the number of points used in the FFT, which
increases the spectral frequency resolution (smaller binwidth) but for a limited data set
of N points, as is the case in the example, decreases the number of data averages
thereby decreasing the spectral amplitude resolution (accuracy). keep in mind that
just changing nFFT without making a corresponding change in overlap will oversample
the data, thereby skewing the amplitudes.

in any case, the amplitude change is not approaching infinity, even if you set nFFT to
6000, which is the length of the timeseries, the amplitudes are ~35dB, adjust variable ymax
to see this.

to review issues of spectral/amplitude resolution, windowing/overlap, etc, a good
reference is Random Data by Bendat  &Piersol:
http://www.amazon.com/Random-Data-Analysis-Measurement-Procedures/dp/0471317330
;
i remain unconvinced that the PSD amplitudes are reasonable, which only leaves Matlab
as an alternative... that's a hard pill to swallow... matplotlib is clearly preferable.

brett.mcsweeney@... wrote:

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.

Joseph Park <jpark@...>
Sent by: matplotlib-users-bounces@...

26/10/2007 10:05 AM

To
matplotlib-users@...
cc

Subject
Re: [Matplotlib-users] PSD amplitudes

Shouldn't the PSD for a simple sine wave tend to infinity

the spectral resolution will impact the amplitude, if you
are not dealing with a density. by definition a spectral density
has applied the bandwidth resolution correction. the PSD amplitude
should correspond to the RMS amplitude of the sine wave. in the
example a 1VRMS amplitude sine wave (time domain) should have a
PSD power of 20*log(1V) = 0dB. The windowing function will impact
this ideal number a bit, but certainly not by 25dB.

brett.mcsweeney@... wrote:

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)?

Joseph Park <jpark@...>
Sent by:
matplotlib-users-bounces@...

26/10/2007 06:50 AM

To
matplotlib-users@...
cc

Subject
[Matplotlib-users] PSD amplitudes

Please try the attached script.
The answer should be ~0 dB for each of the frequencies.
Most likely a simple scaling issue/parameter of which i'm ignorant.

--

______________________________________________________________________
This email has been scanned by the MessageLabs Email Security System.
http://www.messagelabs.com/email
______________________________________________________________________
##----------------------------------------------------------------------------
## Name:     psd_scale.py
##
## Purpose:  Test Power Spectral Density of 1Vrms data
##           Depends on Python SciPy and NumPy
##
## Author:       J Park
##
## Created:      10/17/07
##
## Modified:
##----------------------------------------------------------------------------

try:
from numpy import *  #
www.numpy.org numpy.scipy.org
except ImportError:
print "Failed to import numpy."

try:
import pylab as mp  # matplotlib.sourceforge.net
from matplotlib.font_manager import fontManager, FontProperties
except ImportError:
print "Failed to import pylab."

# Default Parameters
nFFT          = 1024
overlap       = 512
freqSample    = 100.
PlotAll       = False
WriteOutput   = False

##----------------------------------------------------------------------------
## Main module
def main():

deltaF = freqSample/nFFT # Frequency resolution in Hz
deltaT = 1./freqSample   # Sample interval
print 'Sample interval %e (s)'       % (deltaT)
print 'Frequency resolution %e (Hz)' % (deltaF)

# Setup Plots
# ----------------------------------------------------------------------
mp.figure(1)
mp.title ( "PSD" )
mp.ylabel( "(dB)" )
mp.xlabel( "Frequency (Hz)" )
legendFont = FontProperties(size='small')

ymin = 0
ymax = 30
xmin = 0
xmax = 50
xticks = 5
yticks = 5

if PlotAll:
mp.figure(2)
mp.title ( "Input Timeseries" )
mp.ylabel( "Amplitude" )
mp.xlabel( "time (s)" )

# Create some synthetic data with unity RMS amplitude = 0 dB
# ----------------------------------------------------------------------
t = mp.arange(0., 60., deltaT) # 60 seconds at deltaT interval
A = 1.414

y0 = A * sin( 2. * math.pi * 5  * t )
y1 = A * sin( 2. * math.pi * 10 * t )
y2 = A * sin( 2. * math.pi * 20 * t )
y3 = A * sin( 2. * math.pi * 30 * t )
y4 = A * sin( 2. * math.pi * 40 * t )
y5 = A * sin( 2. * math.pi * 45 * t )

dataList = [ y0, y1, y2, y3, y4, y5 ]

for data in dataList:
inputDataLen = len( data )
numAverages  = math.floor( inputDataLen / (overlap) ) - 1
normalizedRandomError = 1./math.sqrt( numAverages )
print "%d points" % ( inputDataLen ),
print "%d averages" % (numAverages),
print "normalized random error %.3f" % ( normalizedRandomError )

mp.figure(1)
(Pxx, freqs) = mp.psd( data,
NFFT     = nFFT,
Fs       = freqSample,
noverlap = overlap,
lw       = 2,
label    = '' )

Pxx_dB = 10.*log10(Pxx)

if PlotAll:
mp.figure(2)
mp.plot(t, data, label='' )

# Write Output data
# ----------------------------------------------------------------------
if WriteOutput:
PxxLen = len(Pxx)
OutputFile = "PSD.dat"
fdOutFile = open( OutputFile, 'a' )
fdOutFile.write( "Freq\t\tPower(dB)\n" )
for i in range(PxxLen):
fdOutFile.write( "%.4e\t%.3f\n" % ( freqs[i], Pxx_dB[i] ) )
fdOutFile.close()
print "Wrote ", PxxLen, " points to ", OutputFile

# Show the Plot
# ----------------------------------------------------------------------
mp.figure(1)
mp.axis([xmin, xmax, ymin, ymax])
mp.xticks( arange(xmin, xmax+1, xticks) )
mp.yticks( arange(ymin, ymax  , yticks) )
mp.title('')
mp.xlabel('Frequency (Hz)')
mp.ylabel(r'\$\tt{dB re V^2/Hz}\$')
#mp.legend( loc='upper right', prop=legendFont )
if WriteOutput:
plotFileName = "PSD.png"
mp.savefig( plotFileName )
print "Wrote png image to ", plotFileName
if PlotAll:
mp.figure(2)
#mp.legend( loc='lower left', prop=legendFont )
mp.show()

print "Normal Exit"
## Main module
##----------------------------------------------------------------------------

##----------------------------------------------------------------------------
## Provide for cmd line invocation
if __name__ == "__main__":
main()

-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems?  Stop.
Now Search log events and configuration files using AJAX and a browser.
http://get.splunk.com/_______________________________________________
Matplotlib-users mailing list

Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

UNITED GROUP
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.
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.

--

______________________________________________________________________
This email has been scanned by the MessageLabs Email Security System.
______________________________________________________________________
-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems?  Stop.
Now Search log events and configuration files using AJAX and a browser.
; Matplotlib-users mailing list
Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
;

UNITED GROUP
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.
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.

--
```
 Re: [Matplotlib-users] PSD amplitudes From: - 2007-10-26 01:02:23 Attachments: Message as HTML ```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? 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. Joseph Park Sent by: matplotlib-users-bounces@... 26/10/2007 10:49 AM To cc matplotlib-users@... Subject Re: [Matplotlib-users] PSD amplitudes is the suggestion that the matplotlib algorithm is correct in computing PSD amplitudes? btw, increasing nFFT increases the number of points used in the FFT, which increases the spectral frequency resolution (smaller binwidth) but for a limited data set of N points, as is the case in the example, decreases the number of data averages thereby decreasing the spectral amplitude resolution (accuracy). keep in mind that just changing nFFT without making a corresponding change in overlap will oversample the data, thereby skewing the amplitudes. in any case, the amplitude change is not approaching infinity, even if you set nFFT to 6000, which is the length of the timeseries, the amplitudes are ~35dB, adjust variable ymax to see this. to review issues of spectral/amplitude resolution, windowing/overlap, etc, a good reference is Random Data by Bendat &Piersol: http://www.amazon.com/Random-Data-Analysis-Measurement-Procedures/dp/0471317330 i remain unconvinced that the PSD amplitudes are reasonable, which only leaves Matlab as an alternative... that's a hard pill to swallow... matplotlib is clearly preferable. brett.mcsweeney@... wrote: 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. Joseph Park Sent by: matplotlib-users-bounces@... 26/10/2007 10:05 AM To matplotlib-users@... cc Subject Re: [Matplotlib-users] PSD amplitudes Shouldn't the PSD for a simple sine wave tend to infinity the spectral resolution will impact the amplitude, if you are not dealing with a density. by definition a spectral density has applied the bandwidth resolution correction. the PSD amplitude should correspond to the RMS amplitude of the sine wave. in the example a 1VRMS amplitude sine wave (time domain) should have a PSD power of 20*log(1V) = 0dB. The windowing function will impact this ideal number a bit, but certainly not by 25dB. brett.mcsweeney@... wrote: 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)? Joseph Park Sent by: matplotlib-users-bounces@... 26/10/2007 06:50 AM To matplotlib-users@... cc Subject [Matplotlib-users] PSD amplitudes Please try the attached script. The answer should be ~0 dB for each of the frequencies. Most likely a simple scaling issue/parameter of which i'm ignorant. -- ______________________________________________________________________ This email has been scanned by the MessageLabs Email Security System. For more information please visit http://www.messagelabs.com/email ______________________________________________________________________ ##---------------------------------------------------------------------------- ## Name: psd_scale.py ## ## Purpose: Test Power Spectral Density of 1Vrms data ## Depends on Python SciPy and NumPy ## ## Author: J Park ## ## Created: 10/17/07 ## ## Modified: ##---------------------------------------------------------------------------- try: from numpy import * # http://www.numpy.org numpy.scipy.org except ImportError: print "Failed to import numpy." try: import pylab as mp # matplotlib.sourceforge.net from matplotlib.font_manager import fontManager, FontProperties except ImportError: print "Failed to import pylab." # Default Parameters nFFT = 1024 overlap = 512 freqSample = 100. PlotAll = False WriteOutput = False ##---------------------------------------------------------------------------- ## Main module def main(): deltaF = freqSample/nFFT # Frequency resolution in Hz deltaT = 1./freqSample # Sample interval print 'Sample interval %e (s)' % (deltaT) print 'Frequency resolution %e (Hz)' % (deltaF) # Setup Plots # ---------------------------------------------------------------------- mp.figure(1) mp.title ( "PSD" ) mp.ylabel( "(dB)" ) mp.xlabel( "Frequency (Hz)" ) legendFont = FontProperties(size='small') ymin = 0 ymax = 30 xmin = 0 xmax = 50 xticks = 5 yticks = 5 if PlotAll: mp.figure(2) mp.title ( "Input Timeseries" ) mp.ylabel( "Amplitude" ) mp.xlabel( "time (s)" ) # Create some synthetic data with unity RMS amplitude = 0 dB # ---------------------------------------------------------------------- t = mp.arange(0., 60., deltaT) # 60 seconds at deltaT interval A = 1.414 y0 = A * sin( 2. * math.pi * 5 * t ) y1 = A * sin( 2. * math.pi * 10 * t ) y2 = A * sin( 2. * math.pi * 20 * t ) y3 = A * sin( 2. * math.pi * 30 * t ) y4 = A * sin( 2. * math.pi * 40 * t ) y5 = A * sin( 2. * math.pi * 45 * t ) dataList = [ y0, y1, y2, y3, y4, y5 ] for data in dataList: inputDataLen = len( data ) numAverages = math.floor( inputDataLen / (overlap) ) - 1 normalizedRandomError = 1./math.sqrt( numAverages ) print "%d points" % ( inputDataLen ), print "%d averages" % (numAverages), print "normalized random error %.3f" % ( normalizedRandomError ) mp.figure(1) (Pxx, freqs) = mp.psd( data, NFFT = nFFT, Fs = freqSample, noverlap = overlap, lw = 2, label = '' ) Pxx_dB = 10.*log10(Pxx) if PlotAll: mp.figure(2) mp.plot(t, data, label='' ) # Write Output data # ---------------------------------------------------------------------- if WriteOutput: PxxLen = len(Pxx) OutputFile = "PSD.dat" fdOutFile = open( OutputFile, 'a' ) fdOutFile.write( "Freq\t\tPower(dB)\n" ) for i in range(PxxLen): fdOutFile.write( "%.4e\t%.3f\n" % ( freqs[i], Pxx_dB[i] ) ) fdOutFile.close() print "Wrote ", PxxLen, " points to ", OutputFile # Show the Plot # ---------------------------------------------------------------------- mp.figure(1) mp.axis([xmin, xmax, ymin, ymax]) mp.xticks( arange(xmin, xmax+1, xticks) ) mp.yticks( arange(ymin, ymax , yticks) ) mp.title('') mp.xlabel('Frequency (Hz)') mp.ylabel(r'\$\tt{dB re V^2/Hz}\$') #mp.legend( loc='upper right', prop=legendFont ) if WriteOutput: plotFileName = "PSD.png" mp.savefig( plotFileName ) print "Wrote png image to ", plotFileName if PlotAll: mp.figure(2) #mp.legend( loc='lower left', prop=legendFont ) mp.show() print "Normal Exit" ## Main module ##---------------------------------------------------------------------------- ##---------------------------------------------------------------------------- ## Provide for cmd line invocation if __name__ == "__main__": main() ------------------------------------------------------------------------- This SF.net email is sponsored by: Splunk Inc. Still grepping through log files to find problems? Stop. Now Search log events and configuration files using AJAX and a browser. Download your FREE copy of Splunk now >> http://get.splunk.com/_______________________________________________ Matplotlib-users mailing list Matplotlib-users@... https://lists.sourceforge.net/lists/listinfo/matplotlib-users UNITED GROUP 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. 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. -- ______________________________________________________________________ This email has been scanned by the MessageLabs Email Security System. For more information please visit http://www.messagelabs.com/email ______________________________________________________________________ ------------------------------------------------------------------------- This SF.net email is sponsored by: Splunk Inc. Still grepping through log files to find problems? Stop. Now Search log events and configuration files using AJAX and a browser. Download your FREE copy of Splunk now >> http://get.splunk.com/_______________________________________________ Matplotlib-users mailing list Matplotlib-users@... https://lists.sourceforge.net/lists/listinfo/matplotlib-users UNITED GROUP 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. 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. -- ______________________________________________________________________ This email has been scanned by the MessageLabs Email Security System. For more information please visit http://www.messagelabs.com/email ______________________________________________________________________ ------------------------------------------------------------------------- This SF.net email is sponsored by: Splunk Inc. Still grepping through log files to find problems? Stop. Now Search log events and configuration files using AJAX and a browser. Download your FREE copy of Splunk now >> http://get.splunk.com/ _______________________________________________ Matplotlib-users mailing list Matplotlib-users@... https://lists.sourceforge.net/lists/listinfo/matplotlib-users UNITED GROUP 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. 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.```
 Re: [Matplotlib-users] PSD amplitudes From: Joseph Park - 2007-10-26 01:52:10 Attachments: Message as HTML     moz-screenshot-5.jpg ``` spectral density is by convention a 1Hz binwidth, not an arbitrary one, units of A^2/Hz.

perhaps if you manually compute the spectral density of a sine wave, you will easily see
that they don't have infinite power, R is the autocorrelation of the Asin(wt):

Back to the original question:

Is there evidence that the matplotlib PSD spectral amplitudes are accurate?
say by comparison with Matlab results, or a synthetic signal as in the example, or
from considerations of basic DSP as in the references?

brett.mcsweeney@... wrote:

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?

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.

Joseph Park <jpark@...>
Sent by: matplotlib-users-bounces@...

26/10/2007 10:49 AM

To

cc
matplotlib-users@...
Subject
Re: [Matplotlib-users] PSD amplitudes

is the suggestion that the matplotlib algorithm is correct in computing PSD amplitudes?

btw, increasing nFFT increases the number of points used in the FFT, which
increases the spectral frequency resolution (smaller binwidth) but for a limited data set
of N points, as is the case in the example, decreases the number of data averages
thereby decreasing the spectral amplitude resolution (accuracy). keep in mind that
just changing nFFT without making a corresponding change in overlap will oversample
the data, thereby skewing the amplitudes.

in any case, the amplitude change is not approaching infinity, even if you set nFFT to
6000, which is the length of the timeseries, the amplitudes are ~35dB, adjust variable ymax
to see this.

to review issues of spectral/amplitude resolution, windowing/overlap, etc, a good
reference is Random Data by Bendat  &Piersol:

http://www.amazon.com/Random-Data-Analysis-Measurement-Procedures/dp/0471317330

i remain unconvinced that the PSD amplitudes are reasonable, which only leaves Matlab
as an alternative... that's a hard pill to swallow... matplotlib is clearly preferable.

brett.mcsweeney@... wrote:

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.

Joseph Park <jpark@...>
Sent by:
matplotlib-users-bounces@...

26/10/2007 10:05 AM

To
matplotlib-users@...
cc

Subject
Re: [Matplotlib-users] PSD amplitudes

Shouldn't the PSD for a simple sine wave tend to infinity

the spectral resolution will impact the amplitude, if you
are not dealing with a density. by definition a spectral density
has applied the bandwidth resolution correction. the PSD amplitude
should correspond to the RMS amplitude of the sine wave. in the
example a 1VRMS amplitude sine wave (time domain) should have a
PSD power of 20*log(1V) = 0dB. The windowing function will impact
this ideal number a bit, but certainly not by 25dB.

brett.mcsweeney@... wrote:

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)?

Joseph Park <jpark@...>
Sent by:
matplotlib-users-bounces@...

26/10/2007 06:50 AM

To
matplotlib-users@...
cc

Subject
[Matplotlib-users] PSD amplitudes

Please try the attached script.
The answer should be ~0 dB for each of the frequencies.
Most likely a simple scaling issue/parameter of which i'm ignorant.

--

______________________________________________________________________
This email has been scanned by the MessageLabs Email Security System.
http://www.messagelabs.com/email
______________________________________________________________________
##----------------------------------------------------------------------------
## Name:     psd_scale.py
##
## Purpose:  Test Power Spectral Density of 1Vrms data
##           Depends on Python SciPy and NumPy
##
## Author:       J Park
##
## Created:      10/17/07
##
## Modified:
##----------------------------------------------------------------------------

try:
from numpy import *  #
www.numpy.org numpy.scipy.org
except ImportError:
print "Failed to import numpy."

try:
import pylab as mp  # matplotlib.sourceforge.net
from matplotlib.font_manager import fontManager, FontProperties
except ImportError:
print "Failed to import pylab."

# Default Parameters
nFFT          = 1024
overlap       = 512
freqSample    = 100.
PlotAll       = False
WriteOutput   = False

##----------------------------------------------------------------------------
## Main module
def main():

deltaF = freqSample/nFFT # Frequency resolution in Hz
deltaT = 1./freqSample   # Sample interval
print 'Sample interval %e (s)'       % (deltaT)
print 'Frequency resolution %e (Hz)' % (deltaF)

# Setup Plots
# ----------------------------------------------------------------------
mp.figure(1)
mp.title ( "PSD" )
mp.ylabel( "(dB)" )
mp.xlabel( "Frequency (Hz)" )
legendFont = FontProperties(size='small')

ymin = 0
ymax = 30
xmin = 0
xmax = 50
xticks = 5
yticks = 5

if PlotAll:
mp.figure(2)
mp.title ( "Input Timeseries" )
mp.ylabel( "Amplitude" )
mp.xlabel( "time (s)" )

# Create some synthetic data with unity RMS amplitude = 0 dB
# ----------------------------------------------------------------------
t = mp.arange(0., 60., deltaT) # 60 seconds at deltaT interval
A = 1.414

y0 = A * sin( 2. * math.pi * 5  * t )
y1 = A * sin( 2. * math.pi * 10 * t )
y2 = A * sin( 2. * math.pi * 20 * t )
y3 = A * sin( 2. * math.pi * 30 * t )
y4 = A * sin( 2. * math.pi * 40 * t )
y5 = A * sin( 2. * math.pi * 45 * t )

dataList = [ y0, y1, y2, y3, y4, y5 ]

for data in dataList:
inputDataLen = len( data )
numAverages  = math.floor( inputDataLen / (overlap) ) - 1
normalizedRandomError = 1./math.sqrt( numAverages )
print "%d points" % ( inputDataLen ),
print "%d averages" % (numAverages),
print "normalized random error %.3f" % ( normalizedRandomError )

mp.figure(1)
(Pxx, freqs) = mp.psd( data,
NFFT     = nFFT,
Fs       = freqSample,
noverlap = overlap,
lw       = 2,
label    = '' )

Pxx_dB = 10.*log10(Pxx)

if PlotAll:
mp.figure(2)
mp.plot(t, data, label='' )

# Write Output data
# ----------------------------------------------------------------------
if WriteOutput:
PxxLen = len(Pxx)
OutputFile = "PSD.dat"
fdOutFile = open( OutputFile, 'a' )
fdOutFile.write( "Freq\t\tPower(dB)\n" )
for i in range(PxxLen):
fdOutFile.write( "%.4e\t%.3f\n" % ( freqs[i], Pxx_dB[i] ) )
fdOutFile.close()
print "Wrote ", PxxLen, " points to ", OutputFile

# Show the Plot
# ----------------------------------------------------------------------
mp.figure(1)
mp.axis([xmin, xmax, ymin, ymax])
mp.xticks( arange(xmin, xmax+1, xticks) )
mp.yticks( arange(ymin, ymax  , yticks) )
mp.title('')
mp.xlabel('Frequency (Hz)')
mp.ylabel(r'\$\tt{dB re V^2/Hz}\$')
#mp.legend( loc='upper right', prop=legendFont )
if WriteOutput:
plotFileName = "PSD.png"
mp.savefig( plotFileName )
print "Wrote png image to ", plotFileName
if PlotAll:
mp.figure(2)
#mp.legend( loc='lower left', prop=legendFont )
mp.show()

print "Normal Exit"
## Main module
##----------------------------------------------------------------------------

##----------------------------------------------------------------------------
## Provide for cmd line invocation
if __name__ == "__main__":
main()

-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems?  Stop.
Now Search log events and configuration files using AJAX and a browser.
http://get.splunk.com/_______________________________________________
Matplotlib-users mailing list

Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

UNITED GROUP
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.
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.

--

______________________________________________________________________
This email has been scanned by the MessageLabs Email Security System.
http://www.messagelabs.com/email
______________________________________________________________________
-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems?  Stop.
Now Search log events and configuration files using AJAX and a browser.
http://get.splunk.com/_______________________________________________
Matplotlib-users mailing list

Matplotlib-users@...
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

UNITED GROUP
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.
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.

--

______________________________________________________________________
This email has been scanned by the MessageLabs Email Security System.
______________________________________________________________________
-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems?  Stop.
Now Search log events and configuration files using AJAX and a browser.