Revision: 6518
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=6518&view=rev
Author: ryanmay
Date: 2008-12-08 21:15:13 +0000 (Mon, 08 Dec 2008)
Log Message:
-----------
Update spectral methods (psd, csd, etc.) to scale one-sided densities by a
factor of 2 and, optionally, scale all densities by the sampling frequency.
This gives better MatLab compatibility. Update corresponding Axes methods, and
make Axes.psd() label the y-axis of the plot as appropriate.
Modified Paths:
--------------
trunk/matplotlib/lib/matplotlib/axes.py
trunk/matplotlib/lib/matplotlib/mlab.py
Modified: trunk/matplotlib/lib/matplotlib/axes.py
===================================================================
--- trunk/matplotlib/lib/matplotlib/axes.py 2008-12-08 21:06:49 UTC (rev
6517)
+++ trunk/matplotlib/lib/matplotlib/axes.py 2008-12-08 21:15:13 UTC (rev
6518)
@@ -6719,13 +6719,13 @@
def psd(self, x, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
window=mlab.window_hanning, noverlap=0, pad_to=None,
- sides='default', **kwargs):
+ sides='default', scale_by_freq=None, **kwargs):
"""
call signature::
psd(x, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
window=mlab.window_hanning, noverlap=0, pad_to=None,
- sides='default', **kwargs)
+ sides='default', scale_by_freq=None, **kwargs)
The power spectral density by Welch's average periodogram
method. The vector *x* is divided into *NFFT* length
@@ -6764,13 +6764,18 @@
"""
if not self._hold: self.cla()
pxx, freqs = mlab.psd(x, NFFT, Fs, detrend, window, noverlap, pad_to,
- sides)
+ sides, scale_by_freq)
pxx.shape = len(freqs),
freqs += Fc
+ if scale_by_freq in (None, True):
+ psd_units = 'dB/Hz'
+ else:
+ psd_units = 'dB'
+
self.plot(freqs, 10*np.log10(pxx), **kwargs)
self.set_xlabel('Frequency')
- self.set_ylabel('Power Spectrum (dB)')
+ self.set_ylabel('Power Spectral Density (%s)' % psd_units)
self.grid(True)
vmin, vmax = self.viewLim.intervaly
intv = vmax-vmin
@@ -6791,13 +6796,13 @@
def csd(self, x, y, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
window=mlab.window_hanning, noverlap=0, pad_to=None,
- sides='default', **kwargs):
+ sides='default', scale_by_freq=None, **kwargs):
"""
call signature::
csd(x, y, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
window=mlab.window_hanning, noverlap=0, pad_to=None,
- sides='default', **kwargs)
+ sides='default', scale_by_freq=None, **kwargs)
The cross spectral density :math:`P_{xy}` by Welch's average
periodogram method. The vectors *x* and *y* are divided into
@@ -6837,7 +6842,7 @@
"""
if not self._hold: self.cla()
pxy, freqs = mlab.csd(x, y, NFFT, Fs, detrend, window, noverlap,
- pad_to, sides)
+ pad_to, sides, scale_by_freq)
pxy.shape = len(freqs),
# pxy is complex
freqs += Fc
@@ -6859,13 +6864,13 @@
def cohere(self, x, y, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
window=mlab.window_hanning, noverlap=0, pad_to=None,
- sides='default', **kwargs):
+ sides='default', scale_by_freq=None, **kwargs):
"""
call signature::
cohere(x, y, NFFT=256, Fs=2, Fc=0, detrend = mlab.detrend_none,
window = mlab.window_hanning, noverlap=0, pad_to=None,
- sides='default', **kwargs)
+ sides='default', scale_by_freq=None, **kwargs)
cohere the coherence between *x* and *y*. Coherence is the normalized
cross spectral density:
@@ -6902,7 +6907,8 @@
.. plot:: mpl_examples/pylab_examples/cohere_demo.py
"""
if not self._hold: self.cla()
- cxy, freqs = mlab.cohere(x, y, NFFT, Fs, detrend, window, noverlap)
+ cxy, freqs = mlab.cohere(x, y, NFFT, Fs, detrend, window, noverlap,
+ scale_by_freq)
freqs += Fc
self.plot(freqs, cxy, **kwargs)
@@ -6915,13 +6921,15 @@
def specgram(self, x, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
window=mlab.window_hanning, noverlap=128,
- cmap=None, xextent=None, pad_to=None, sides='default'):
+ cmap=None, xextent=None, pad_to=None, sides='default',
+ scale_by_freq=None):
"""
call signature::
specgram(x, NFFT=256, Fs=2, Fc=0, detrend=mlab.detrend_none,
window=mlab.window_hanning, noverlap=128,
- cmap=None, xextent=None, pad_to=None, sides='default')
+ cmap=None, xextent=None, pad_to=None, sides='default',
+ scale_by_freq=None)
Compute a spectrogram of data in *x*. Data are split into
*NFFT* length segments and the PSD of each section is
@@ -6965,7 +6973,7 @@
if not self._hold: self.cla()
Pxx, freqs, bins = mlab.specgram(x, NFFT, Fs, detrend,
- window, noverlap, pad_to, sides)
+ window, noverlap, pad_to, sides, scale_by_freq)
Z = 10. * np.log10(Pxx)
Z = np.flipud(Z)
Modified: trunk/matplotlib/lib/matplotlib/mlab.py
===================================================================
--- trunk/matplotlib/lib/matplotlib/mlab.py 2008-12-08 21:06:49 UTC (rev
6517)
+++ trunk/matplotlib/lib/matplotlib/mlab.py 2008-12-08 21:15:13 UTC (rev
6518)
@@ -243,14 +243,15 @@
#This is a helper function that implements the commonality between the
#psd, csd, and spectrogram. It is *NOT* meant to be used outside of mlab
def _spectral_helper(x, y, NFFT=256, Fs=2, detrend=detrend_none,
- window=window_hanning, noverlap=0, pad_to=None, sides='default'):
+ window=window_hanning, noverlap=0, pad_to=None, sides='default',
+ scale_by_freq=None):
#The checks for if y is x are so that we can use the same function to
#implement the core of psd(), csd(), and spectrogram() without doing
#extra calculations. We return the unaveraged Pxy, freqs, and t.
+ same_data = y is x
#Make sure we're dealing with a numpy array. If y and x were the same
#object to start with, keep them that way
- same_data = y is x
x = np.asarray(x)
if not same_data:
@@ -270,15 +271,31 @@
if pad_to is None:
pad_to = NFFT
+ if scale_by_freq is None:
+ warnings.warn("psd, csd, and specgram have changed to scale their "
+ "densities by the sampling frequency for better MatLab "
+ "compatibility. You can pass scale_by_freq=False to disable "
+ "this behavior. Also, one-sided densities are scaled by a "
+ "factor of 2.")
+ scale_by_freq = True
+
# For real x, ignore the negative frequencies unless told otherwise
if (sides == 'default' and np.iscomplexobj(x)) or sides == 'twosided':
numFreqs = pad_to
+ scaling_factor = 1.
elif sides in ('default', 'onesided'):
numFreqs = pad_to//2 + 1
+ scaling_factor = 2.
else:
raise ValueError("sides must be one of: 'default', 'onesided', or "
"'twosided'")
+ # Matlab divides by the sampling frequency so that density function
+ # has units of dB/Hz and can be integrated by the plotted frequency
+ # values. Perform the same scaling here.
+ if scale_by_freq:
+ scaling_factor /= Fs
+
if cbook.iterable(window):
assert(len(window) == NFFT)
windowVals = window
@@ -305,8 +322,10 @@
Pxy[:,i] = np.conjugate(fx[:numFreqs]) * fy[:numFreqs]
# Scale the spectrum by the norm of the window to compensate for
- # windowing loss; see Bendat & Piersol Sec 11.5.2
- Pxy /= (np.abs(windowVals)**2).sum()
+ # windowing loss; see Bendat & Piersol Sec 11.5.2. Also include
+ # scaling factors for one-sided densities and dividing by the sampling
+ # frequency, if desired.
+ Pxy *= scaling_factor / (np.abs(windowVals)**2).sum()
t = 1./Fs * (ind + NFFT / 2.)
freqs = float(Fs) / pad_to * np.arange(numFreqs)
@@ -363,12 +382,18 @@
*sides*: [ 'default' | 'onesided' | 'twosided' ]
Specifies which sides of the PSD to return. Default gives the
default behavior, which returns one-sided for real data and both
- for complex data. 'one' forces the return of a one-sided PSD, while
- 'both' forces two-sided.
+ for complex data. 'onesided' forces the return of a one-sided PSD,
+ while 'twosided' forces two-sided.
+
+ *scale_by_freq*: boolean
+ Specifies whether the resulting density values should be scaled
+ by the scaling frequency, which gives density in units of Hz^-1.
+ This allows for integration over the returned frequency values.
+ The default is True for MatLab compatibility.
"""
-def psd(x, NFFT=256, Fs=2, detrend=detrend_none,
- window=window_hanning, noverlap=0, pad_to=None, sides='default'):
+def psd(x, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning,
+ noverlap=0, pad_to=None, sides='default', scale_by_freq=None):
"""
The power spectral density by Welch's average periodogram method.
The vector *x* is divided into *NFFT* length blocks. Each block
@@ -388,13 +413,14 @@
Bendat & Piersol -- Random Data: Analysis and Measurement
Procedures, John Wiley & Sons (1986)
"""
- Pxx,freqs = csd(x, x, NFFT, Fs, detrend, window, noverlap, pad_to, sides)
+ Pxx,freqs = csd(x, x, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
+ scale_by_freq)
return Pxx.real,freqs
psd.__doc__ = psd.__doc__ % kwdocd
-def csd(x, y, NFFT=256, Fs=2, detrend=detrend_none,
- window=window_hanning, noverlap=0, pad_to=None, sides='default'):
+def csd(x, y, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning,
+ noverlap=0, pad_to=None, sides='default', scale_by_freq=None):
"""
The cross power spectral density by Welch's average periodogram
method. The vectors *x* and *y* are divided into *NFFT* length
@@ -417,7 +443,7 @@
Procedures, John Wiley & Sons (1986)
"""
Pxy, freqs, t = _spectral_helper(x, y, NFFT, Fs, detrend, window,
- noverlap, pad_to, sides)
+ noverlap, pad_to, sides, scale_by_freq)
if len(Pxy.shape) == 2 and Pxy.shape[1]>1:
Pxy = Pxy.mean(axis=1)
@@ -425,9 +451,8 @@
csd.__doc__ = csd.__doc__ % kwdocd
-def specgram(x, NFFT=256, Fs=2, detrend=detrend_none,
- window=window_hanning, noverlap=128, pad_to=None,
- sides='default'):
+def specgram(x, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning,
+ noverlap=128, pad_to=None, sides='default', scale_by_freq=None):
"""
Compute a spectrogram of data in *x*. Data are split into *NFFT*
length segements and the PSD of each section is computed. The
@@ -458,7 +483,7 @@
assert(NFFT > noverlap)
Pxx, freqs, t = _spectral_helper(x, x, NFFT, Fs, detrend, window,
- noverlap, pad_to, sides)
+ noverlap, pad_to, sides, scale_by_freq)
Pxx = Pxx.real #Needed since helper implements generically
if (np.iscomplexobj(x) and sides == 'default') or sides == 'twosided':
@@ -473,8 +498,8 @@
_coh_error = """Coherence is calculated by averaging over *NFFT*
length segments. Your signal is too short for your choice of *NFFT*.
"""
-def cohere(x, y, NFFT=256, Fs=2, detrend=detrend_none,
- window=window_hanning, noverlap=0, pad_to=None, sides='default'):
+def cohere(x, y, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning,
+ noverlap=0, pad_to=None, sides='default', scale_by_freq=None):
"""
The coherence between *x* and *y*. Coherence is the normalized
cross spectral density:
@@ -487,7 +512,9 @@
Array or sequence containing the data
%(PSD)s
The return value is the tuple (*Cxy*, *f*), where *f* are the
- frequencies of the coherence vector.
+ frequencies of the coherence vector. For cohere, scaling the
+ individual densities by the sampling frequency has no effect, since
+ the factors cancel out.
.. seealso::
:func:`psd` and :func:`csd`:
@@ -497,9 +524,12 @@
if len(x)<2*NFFT:
raise ValueError(_coh_error)
- Pxx, f = psd(x, NFFT, Fs, detrend, window, noverlap, pad_to, sides)
- Pyy, f = psd(y, NFFT, Fs, detrend, window, noverlap, pad_to, sides)
- Pxy, f = csd(x, y, NFFT, Fs, detrend, window, noverlap, pad_to, sides)
+ Pxx, f = psd(x, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
+ scale_by_freq)
+ Pyy, f = psd(y, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
+ scale_by_freq)
+ Pxy, f = csd(x, y, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
+ scale_by_freq)
Cxy = np.divide(np.absolute(Pxy)**2, Pxx*Pyy)
Cxy.shape = (len(f),)
@@ -3246,4 +3276,3 @@
c2x, c2y = c1x + 1./3. * (q2x - q0x), c1y + 1./3. * (q2y - q0y)
# c3x, c3y = q2x, q2y
return q0x, q0y, c1x, c1y, c2x, c2y, q2x, q2y
-
This was sent by the SourceForge.net collaborative development platform, the
world's largest Open Source development site.
------------------------------------------------------------------------------
SF.Net email is Sponsored by MIX09, March 18-20, 2009 in Las Vegas, Nevada.
The future of the web can't happen without you. Join us at MIX09 to help
pave the way to the Next Web now. Learn more and register at
http://ad.doubleclick.net/clk;208669438;13503038;i?http://2009.visitmix.com/
_______________________________________________
Matplotlib-checkins mailing list
Matplotlib-checkins@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-checkins
