PSD amplitudes

Joseph_Park · October 26, 2007, 12:05am

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.

wrote:

···

brett.mcsweeney@…1537…

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:
26/10/2007 06:50 AM

To

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 ______________________________________________________________________`##----------------------------------------------------------------------------

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 * # 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 >>
Matplotlib-users mailing list
`

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.
<jpark@…1765…>matplotlib-users-bounces@lists.sourceforge.netmatplotlib-users@lists.sourceforge.nethttp://www.messagelabs.com/email
www.numpy.org

http://get.splunk.com/_______________________________________________

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

--

system · October 26, 2007, 12:14am

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@…1765…>

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 10:05 AM

To

matplotlib-users@lists.sourceforge.net
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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 06:50 AM

To

matplotlib-users@lists.sourceforge.net
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 * # [www.numpy.org](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/_______________________________________________](http://get.splunk.com/_______________________________________________)

Matplotlib-users mailing list``

Matplotlib-users@lists.sourceforge.net

[https://lists.sourceforge.net/lists/listinfo/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@lists.sourceforge.net

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.

Joseph_Park · October 26, 2007, 12:49am

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:

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.

wrote:

···

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

brett.mcsweeney@…1537…

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:
26/10/2007 10:05 AM

To

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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net

26/10/2007
06:50 AM

To

matplotlib-users@lists.sourceforge.net

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 * # [www.numpy.org](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 >> [Free Trials and Downloads | Splunk](Free Trials and Downloads | Splunk)`

Matplotlib-users mailing list``

Matplotlib-users@lists.sourceforge.net

[https://lists.sourceforge.net/lists/listinfo/matplotlib-users`](matplotlib-users List Signup and Options)

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 ______________________________________________________________________`-------------------------------------------------------------------------

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 >>
Matplotlib-users mailing list
`

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.
<jpark@…1765…>matplotlib-users-bounces@lists.sourceforge.netmatplotlib-users@lists.sourceforge.nethttp://www.messagelabs.com/email
http://get.splunk.com/_______________________________________________

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

--

system · October 26, 2007, 12:59am

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@…1765…>

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 10:49 AM

To

cc

matplotlib-users@lists.sourceforge.net
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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 10:05 AM

To

matplotlib-users@lists.sourceforge.net
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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 06:50 AM

To

matplotlib-users@lists.sourceforge.net
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 * # [www.numpy.org](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/_______________________________________________](http://get.splunk.com/_______________________________________________)

Matplotlib-users mailing list`

Matplotlib-users@lists.sourceforge.net

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/_______________________________________________](http://get.splunk.com/_______________________________________________)

Matplotlib-users mailing list``

Matplotlib-users@lists.sourceforge.net

[https://lists.sourceforge.net/lists/listinfo/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@lists.sourceforge.net

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.

Joseph_Park · October 26, 2007, 1:50am

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?

wrote:

···

brett.mcsweeney@…1537…

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:
26/10/2007 10:49 AM

To

cc

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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net

26/10/2007
10:05 AM

To

matplotlib-users@lists.sourceforge.net

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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net

26/10/2007
06:50 AM

To

matplotlib-users@lists.sourceforge.net

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 * # [www.numpy.org](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 >> [Free Trials and Downloads | Splunk](Free Trials and Downloads | Splunk)`

Matplotlib-users mailing list`

Matplotlib-users@lists.sourceforge.net

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 >> [Free Trials and Downloads | Splunk](Free Trials and Downloads | Splunk)`

Matplotlib-users mailing list``

Matplotlib-users@lists.sourceforge.net

[https://lists.sourceforge.net/lists/listinfo/matplotlib-users`](matplotlib-users List Signup and Options)

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 ______________________________________________________________________`-------------------------------------------------------------------------

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 >>
Matplotlib-users mailing list
`

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.
<jpark@…1765…>matplotlib-users-bounces@lists.sourceforge.netmatplotlib-users@lists.sourceforge.nethttp://www.messagelabs.com/email
http://get.splunk.com/_______________________________________________

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

--

system · October 26, 2007, 2:42am

I didn’t say infinite power, but infinite
power density at the sine wave frequemcy.

Being per Hz doesn’t mean that one computes
the PSD using a 1 Hz band! It means that one divides the power in
the band by the width of the band, which can be anything one chooses.

The formula for S(f) of a sine wave
is a delta function!

Joseph Park <jpark@…1765…>

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 11:50 AM

To

cc

matplotlib-users@lists.sourceforge.net
Subject

Re: [Matplotlib-users] PSD amplitudes

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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 10:49 AM

To

cc

matplotlib-users@lists.sourceforge.net
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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 10:05 AM

To

matplotlib-users@lists.sourceforge.net
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@…1537…
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@…1765…>**

Sent by: matplotlib-users-bounces@lists.sourceforge.net
26/10/2007 06:50 AM

To

matplotlib-users@lists.sourceforge.net
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 * # [www.numpy.org](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/_______________________________________________](http://get.splunk.com/_______________________________________________)

Matplotlib-users mailing list`

Matplotlib-users@lists.sourceforge.net

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/_______________________________________________](http://get.splunk.com/_______________________________________________)

Matplotlib-users mailing list`

Matplotlib-users@lists.sourceforge.net

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/_______________________________________________](http://get.splunk.com/_______________________________________________)

Matplotlib-users mailing list``

Matplotlib-users@lists.sourceforge.net

[https://lists.sourceforge.net/lists/listinfo/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@lists.sourceforge.net

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.

PSD amplitudes

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:

Default Parameters

Main module

Setup Plots

Create some synthetic data with unity RMS amplitude = 0

Show the Plot

Main module

Provide for cmd line invocation

if name == “main”: main()

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:

Default Parameters

Main module

Setup Plots

----------------------------------------------------------------------

Create some synthetic data with unity RMS amplitude = 0 dB

----------------------------------------------------------------------

Show the Plot

----------------------------------------------------------------------

Main module

Provide for cmd line invocation

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:

Default Parameters

Main module

Setup Plots

Create some synthetic data with unity RMS amplitude = 0 dB

Show the Plot

Main module

Provide for cmd line invocation

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:

Default Parameters

Main module

Setup Plots

----------------------------------------------------------------------

Create some synthetic data with unity RMS amplitude = 0 dB

----------------------------------------------------------------------

Show the Plot

----------------------------------------------------------------------

Main module

Provide for cmd line invocation

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:

Default Parameters

Main module

Setup Plots

Create some synthetic data with unity RMS amplitude = 0 dB

Show the Plot

Main module

Provide for cmd line invocation

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:

Default Parameters

if name == “main”:
main()