Menu
▾
▴
SF.net SVN: matplotlib:[6392] trunk/matplotlib
|
From: <ry...@us...> - 2008-11-11 18:42:16
|
Revision: 6392
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=6392&view=rev
Author: ryanmay
Date: 2008-11-11 18:42:11 +0000 (Tue, 11 Nov 2008)
Log Message:
-----------
Update the mlab.csd() to match the new options added to mlab.psd(). Factor out the keyword argument docs so that they can be shared between the two functions.
Modified Paths:
--------------
trunk/matplotlib/CHANGELOG
trunk/matplotlib/lib/matplotlib/mlab.py
Modified: trunk/matplotlib/CHANGELOG
===================================================================
--- trunk/matplotlib/CHANGELOG 2008-11-11 18:39:46 UTC (rev 6391)
+++ trunk/matplotlib/CHANGELOG 2008-11-11 18:42:11 UTC (rev 6392)
@@ -1,3 +1,7 @@
+2008-11-11 Update the mlab.csd() to match the new options added to
+ mlab.psd(). Factor out the keyword argument docs so that
+ they can be shared between the two functions. - RM
+
2008-11-11 Update the axes.psd() method to reflect the new options
available in mlab.psd. Add an example to show how the
options change things. - RM
Modified: trunk/matplotlib/lib/matplotlib/mlab.py
===================================================================
--- trunk/matplotlib/lib/matplotlib/mlab.py 2008-11-11 18:39:46 UTC (rev 6391)
+++ trunk/matplotlib/lib/matplotlib/mlab.py 2008-11-11 18:42:11 UTC (rev 6392)
@@ -252,51 +252,7 @@
*x*
Array or sequence containing the data
-
- *NFFT*
- The number of data points used in each block for the FFT.
- Must be even; a power 2 is most efficient. The default value is 256.
-
- *Fs*
- The sampling frequency (samples per time unit). It is used
- to calculate the Fourier frequencies, freqs, in cycles per time
- unit. The default value is 2.
-
- *detrend*
- Any callable function (unlike in matlab where it is a vector).
- For examples, see :func:`detrend`, :func:`detrend_none`, and
- :func:`detrend_mean`. The default is :func:`detrend_none`.
-
- *window*
- A function or a vector of length *NFFT*. To create window
- vectors see :func:`window_hanning`, :func:`window_none`,
- :func:`numpy.blackman`, :func:`numpy.hamming`,
- :func:`numpy.bartlett`, :func:`scipy.signal`,
- :func:`scipy.signal.get_window`, etc. The default is
- :func:`window_hanning`. If a function is passed as the
- argument, it must take a data segment as an argument and
- return the windowed version of the segment.
-
- *noverlap*
- The number of points of overlap between blocks. The default value
- is 0 (no overlap).
-
- *pad_to*
- The number of points to which the data segment is padd when
- performing the FFT. This can be different from *NFFT*, which
- specifies the number of data points used. While not increasing
- the actual resolution of the psd (the minimum distance between
- resolvable peaks), this can give more points in the plot,
- allowing for more detail. This corresponds to the *n* parameter
- in the call to fft(). The default is None, which sets *pad_to*
- equal to *NFFT*
-
- *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.
-
+ %(PSD)s
Returns the tuple (*Pxx*, *freqs*).
Refs:
@@ -314,14 +270,14 @@
if pad_to is None:
pad_to = NFFT
- # for real x, ignore the negative frequencies
+ # For real x, ignore the negative frequencies unless told otherwise
if (sides == 'default' and np.iscomplexobj(x)) or sides == 'twosided':
numFreqs = pad_to
elif sides in ('default', 'onesided'):
numFreqs = pad_to//2 + 1
else:
raise ValueError("sides must be one of: 'default', 'onesided', or "
- "twosided")
+ "'twosided'")
if cbook.iterable(window):
assert(len(window) == NFFT)
@@ -351,23 +307,9 @@
return Pxx, freqs
-def csd(x, y, NFFT=256, Fs=2, detrend=detrend_none,
- window=window_hanning, noverlap=0):
- """
- The cross power spectral density by Welch's average periodogram
- method. The vectors *x* and *y* are divided into *NFFT* length
- blocks. Each block is detrended by the function *detrend* and
- windowed by the function *window*. *noverlap* gives the length
- of the overlap between blocks. The product of the direct FFTs
- of *x* and *y* are averaged over each segment to compute *Pxy*,
- with a scaling to correct for power loss due to windowing.
-
- If len(*x*) < *NFFT* or len(*y*) < *NFFT*, they will be zero
- padded to *NFFT*.
-
- *x*, *y*
- Array or sequence containing the data
-
+#Split out these keyword docs so that they can be used elsewhere
+kwdocd = dict()
+kwdocd['PSD'] ="""
*NFFT*
The number of data points used in each block for the FFT.
Must be even; a power 2 is most efficient. The default value is 256.
@@ -388,12 +330,49 @@
:func:`numpy.blackman`, :func:`numpy.hamming`,
:func:`numpy.bartlett`, :func:`scipy.signal`,
:func:`scipy.signal.get_window`, etc. The default is
- :func:`window_hanning`.
+ :func:`window_hanning`. If a function is passed as the
+ argument, it must take a data segment as an argument and
+ return the windowed version of the segment.
*noverlap*
The number of points of overlap between blocks. The default value
is 0 (no overlap).
+ *pad_to*
+ The number of points to which the data segment is padd when
+ performing the FFT. This can be different from *NFFT*, which
+ specifies the number of data points used. While not increasing
+ the actual resolution of the psd (the minimum distance between
+ resolvable peaks), this can give more points in the plot,
+ allowing for more detail. This corresponds to the *n* parameter
+ in the call to fft(). The default is None, which sets *pad_to*
+ equal to *NFFT*
+
+ *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.
+"""
+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'):
+ """
+ The cross power spectral density by Welch's average periodogram
+ method. The vectors *x* and *y* are divided into *NFFT* length
+ blocks. Each block is detrended by the function *detrend* and
+ windowed by the function *window*. *noverlap* gives the length
+ of the overlap between blocks. The product of the direct FFTs
+ of *x* and *y* are averaged over each segment to compute *Pxy*,
+ with a scaling to correct for power loss due to windowing.
+
+ If len(*x*) < *NFFT* or len(*y*) < *NFFT*, they will be zero
+ padded to *NFFT*.
+
+ *x*, *y*
+ Array or sequence containing the data
+ %(PSD)s
Returns the tuple (*Pxy*, *freqs*).
Refs:
@@ -401,9 +380,6 @@
Procedures, John Wiley & Sons (1986)
"""
- if NFFT % 2:
- raise ValueError, 'NFFT must be even'
-
x = np.asarray(x) # make sure we're dealing with a numpy array
y = np.asarray(y) # make sure we're dealing with a numpy array
@@ -417,40 +393,50 @@
y = np.resize(y, (NFFT,))
y[n:] = 0
- # for real x, ignore the negative frequencies
- if np.iscomplexobj(x): numFreqs = NFFT
- else: numFreqs = NFFT//2+1
+ if pad_to is None:
+ pad_to = NFFT
+ # For real x, ignore the negative frequencies unless told otherwise
+ if (sides == 'default' and np.iscomplexobj(x)) or sides == 'twosided':
+ numFreqs = pad_to
+ elif sides in ('default', 'onesided'):
+ numFreqs = pad_to//2 + 1
+ else:
+ raise ValueError("sides must be one of: 'default', 'onesided', or "
+ "'twosided'")
+
if cbook.iterable(window):
assert(len(window) == NFFT)
windowVals = window
else:
windowVals = window(np.ones((NFFT,), x.dtype))
- step = NFFT-noverlap
- ind = range(0,len(x)-NFFT+1,step)
+
+ step = NFFT - noverlap
+ ind = range(0, len(x) - NFFT + 1, step)
n = len(ind)
Pxy = np.zeros((numFreqs,n), np.complex_)
# do the ffts of the slices
for i in range(n):
thisX = x[ind[i]:ind[i]+NFFT]
- thisX = windowVals*detrend(thisX)
+ thisX = windowVals * detrend(thisX)
thisY = y[ind[i]:ind[i]+NFFT]
- thisY = windowVals*detrend(thisY)
- fx = np.fft.fft(thisX)
- fy = np.fft.fft(thisY)
- Pxy[:,i] = np.conjugate(fx[:numFreqs])*fy[:numFreqs]
+ thisY = windowVals * detrend(thisY)
+ fx = np.fft.fft(thisX, n=pad_to)
+ fy = np.fft.fft(thisY, n=pad_to)
+ 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
if n>1:
Pxy = Pxy.mean(axis=1)
+
Pxy /= (np.abs(windowVals)**2).sum()
- freqs = Fs/NFFT*np.arange(numFreqs)
+ freqs = float(Fs) / pad_to * np.arange(numFreqs)
return Pxy, freqs
+csd.__doc__ = csd.__doc__ % kwdocd
+
def specgram(x, NFFT=256, Fs=2, detrend=detrend_none,
window=window_hanning, noverlap=128):
"""
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
