Menu
▾
▴
SF.net SVN: matplotlib:[6518] trunk/matplotlib/lib/matplotlib
|
From: <ry...@us...> - 2008-12-08 21:15:23
|
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.
|
