Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/matplotlib/mlab.py: 15%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Numerical Python functions written for compatibility with MATLAB
3commands with the same names. Most numerical Python functions can be found in
4the `NumPy`_ and `SciPy`_ libraries. What remains here is code for performing
5spectral computations and kernel density estimations.
7.. _NumPy: https://numpy.org
8.. _SciPy: https://www.scipy.org
10Spectral functions
11------------------
13`cohere`
14 Coherence (normalized cross spectral density)
16`csd`
17 Cross spectral density using Welch's average periodogram
19`detrend`
20 Remove the mean or best fit line from an array
22`psd`
23 Power spectral density using Welch's average periodogram
25`specgram`
26 Spectrogram (spectrum over segments of time)
28`complex_spectrum`
29 Return the complex-valued frequency spectrum of a signal
31`magnitude_spectrum`
32 Return the magnitude of the frequency spectrum of a signal
34`angle_spectrum`
35 Return the angle (wrapped phase) of the frequency spectrum of a signal
37`phase_spectrum`
38 Return the phase (unwrapped angle) of the frequency spectrum of a signal
40`detrend_mean`
41 Remove the mean from a line.
43`detrend_linear`
44 Remove the best fit line from a line.
46`detrend_none`
47 Return the original line.
48"""
50import functools
51from numbers import Number
53import numpy as np
55from matplotlib import _api, _docstring, cbook
58def window_hanning(x):
59 """
60 Return *x* times the Hanning (or Hann) window of len(*x*).
62 See Also
63 --------
64 window_none : Another window algorithm.
65 """
66 return np.hanning(len(x))*x
69def window_none(x):
70 """
71 No window function; simply return *x*.
73 See Also
74 --------
75 window_hanning : Another window algorithm.
76 """
77 return x
80def detrend(x, key=None, axis=None):
81 """
82 Return *x* with its trend removed.
84 Parameters
85 ----------
86 x : array or sequence
87 Array or sequence containing the data.
89 key : {'default', 'constant', 'mean', 'linear', 'none'} or function
90 The detrending algorithm to use. 'default', 'mean', and 'constant' are
91 the same as `detrend_mean`. 'linear' is the same as `detrend_linear`.
92 'none' is the same as `detrend_none`. The default is 'mean'. See the
93 corresponding functions for more details regarding the algorithms. Can
94 also be a function that carries out the detrend operation.
96 axis : int
97 The axis along which to do the detrending.
99 See Also
100 --------
101 detrend_mean : Implementation of the 'mean' algorithm.
102 detrend_linear : Implementation of the 'linear' algorithm.
103 detrend_none : Implementation of the 'none' algorithm.
104 """
105 if key is None or key in ['constant', 'mean', 'default']:
106 return detrend(x, key=detrend_mean, axis=axis)
107 elif key == 'linear':
108 return detrend(x, key=detrend_linear, axis=axis)
109 elif key == 'none':
110 return detrend(x, key=detrend_none, axis=axis)
111 elif callable(key):
112 x = np.asarray(x)
113 if axis is not None and axis + 1 > x.ndim:
114 raise ValueError(f'axis(={axis}) out of bounds')
115 if (axis is None and x.ndim == 0) or (not axis and x.ndim == 1):
116 return key(x)
117 # try to use the 'axis' argument if the function supports it,
118 # otherwise use apply_along_axis to do it
119 try:
120 return key(x, axis=axis)
121 except TypeError:
122 return np.apply_along_axis(key, axis=axis, arr=x)
123 else:
124 raise ValueError(
125 f"Unknown value for key: {key!r}, must be one of: 'default', "
126 f"'constant', 'mean', 'linear', or a function")
129def detrend_mean(x, axis=None):
130 """
131 Return *x* minus the mean(*x*).
133 Parameters
134 ----------
135 x : array or sequence
136 Array or sequence containing the data
137 Can have any dimensionality
139 axis : int
140 The axis along which to take the mean. See `numpy.mean` for a
141 description of this argument.
143 See Also
144 --------
145 detrend_linear : Another detrend algorithm.
146 detrend_none : Another detrend algorithm.
147 detrend : A wrapper around all the detrend algorithms.
148 """
149 x = np.asarray(x)
151 if axis is not None and axis+1 > x.ndim:
152 raise ValueError('axis(=%s) out of bounds' % axis)
154 return x - x.mean(axis, keepdims=True)
157def detrend_none(x, axis=None):
158 """
159 Return *x*: no detrending.
161 Parameters
162 ----------
163 x : any object
164 An object containing the data
166 axis : int
167 This parameter is ignored.
168 It is included for compatibility with detrend_mean
170 See Also
171 --------
172 detrend_mean : Another detrend algorithm.
173 detrend_linear : Another detrend algorithm.
174 detrend : A wrapper around all the detrend algorithms.
175 """
176 return x
179def detrend_linear(y):
180 """
181 Return *x* minus best fit line; 'linear' detrending.
183 Parameters
184 ----------
185 y : 0-D or 1-D array or sequence
186 Array or sequence containing the data
188 See Also
189 --------
190 detrend_mean : Another detrend algorithm.
191 detrend_none : Another detrend algorithm.
192 detrend : A wrapper around all the detrend algorithms.
193 """
194 # This is faster than an algorithm based on linalg.lstsq.
195 y = np.asarray(y)
197 if y.ndim > 1:
198 raise ValueError('y cannot have ndim > 1')
200 # short-circuit 0-D array.
201 if not y.ndim:
202 return np.array(0., dtype=y.dtype)
204 x = np.arange(y.size, dtype=float)
206 C = np.cov(x, y, bias=1)
207 b = C[0, 1]/C[0, 0]
209 a = y.mean() - b*x.mean()
210 return y - (b*x + a)
213def _spectral_helper(x, y=None, NFFT=None, Fs=None, detrend_func=None,
214 window=None, noverlap=None, pad_to=None,
215 sides=None, scale_by_freq=None, mode=None):
216 """
217 Private helper implementing the common parts between the psd, csd,
218 spectrogram and complex, magnitude, angle, and phase spectrums.
219 """
220 if y is None:
221 # if y is None use x for y
222 same_data = True
223 else:
224 # The checks for if y is x are so that we can use the same function to
225 # implement the core of psd(), csd(), and spectrogram() without doing
226 # extra calculations. We return the unaveraged Pxy, freqs, and t.
227 same_data = y is x
229 if Fs is None:
230 Fs = 2
231 if noverlap is None:
232 noverlap = 0
233 if detrend_func is None:
234 detrend_func = detrend_none
235 if window is None:
236 window = window_hanning
238 # if NFFT is set to None use the whole signal
239 if NFFT is None:
240 NFFT = 256
242 if noverlap >= NFFT:
243 raise ValueError('noverlap must be less than NFFT')
245 if mode is None or mode == 'default':
246 mode = 'psd'
247 _api.check_in_list(
248 ['default', 'psd', 'complex', 'magnitude', 'angle', 'phase'],
249 mode=mode)
251 if not same_data and mode != 'psd':
252 raise ValueError("x and y must be equal if mode is not 'psd'")
254 # Make sure we're dealing with a numpy array. If y and x were the same
255 # object to start with, keep them that way
256 x = np.asarray(x)
257 if not same_data:
258 y = np.asarray(y)
260 if sides is None or sides == 'default':
261 if np.iscomplexobj(x):
262 sides = 'twosided'
263 else:
264 sides = 'onesided'
265 _api.check_in_list(['default', 'onesided', 'twosided'], sides=sides)
267 # zero pad x and y up to NFFT if they are shorter than NFFT
268 if len(x) < NFFT:
269 n = len(x)
270 x = np.resize(x, NFFT)
271 x[n:] = 0
273 if not same_data and len(y) < NFFT:
274 n = len(y)
275 y = np.resize(y, NFFT)
276 y[n:] = 0
278 if pad_to is None:
279 pad_to = NFFT
281 if mode != 'psd':
282 scale_by_freq = False
283 elif scale_by_freq is None:
284 scale_by_freq = True
286 # For real x, ignore the negative frequencies unless told otherwise
287 if sides == 'twosided':
288 numFreqs = pad_to
289 if pad_to % 2:
290 freqcenter = (pad_to - 1)//2 + 1
291 else:
292 freqcenter = pad_to//2
293 scaling_factor = 1.
294 elif sides == 'onesided':
295 if pad_to % 2:
296 numFreqs = (pad_to + 1)//2
297 else:
298 numFreqs = pad_to//2 + 1
299 scaling_factor = 2.
301 if not np.iterable(window):
302 window = window(np.ones(NFFT, x.dtype))
303 if len(window) != NFFT:
304 raise ValueError(
305 "The window length must match the data's first dimension")
307 result = np.lib.stride_tricks.sliding_window_view(
308 x, NFFT, axis=0)[::NFFT - noverlap].T
309 result = detrend(result, detrend_func, axis=0)
310 result = result * window.reshape((-1, 1))
311 result = np.fft.fft(result, n=pad_to, axis=0)[:numFreqs, :]
312 freqs = np.fft.fftfreq(pad_to, 1/Fs)[:numFreqs]
314 if not same_data:
315 # if same_data is False, mode must be 'psd'
316 resultY = np.lib.stride_tricks.sliding_window_view(
317 y, NFFT, axis=0)[::NFFT - noverlap].T
318 resultY = detrend(resultY, detrend_func, axis=0)
319 resultY = resultY * window.reshape((-1, 1))
320 resultY = np.fft.fft(resultY, n=pad_to, axis=0)[:numFreqs, :]
321 result = np.conj(result) * resultY
322 elif mode == 'psd':
323 result = np.conj(result) * result
324 elif mode == 'magnitude':
325 result = np.abs(result) / window.sum()
326 elif mode == 'angle' or mode == 'phase':
327 # we unwrap the phase later to handle the onesided vs. twosided case
328 result = np.angle(result)
329 elif mode == 'complex':
330 result /= window.sum()
332 if mode == 'psd':
334 # Also include scaling factors for one-sided densities and dividing by
335 # the sampling frequency, if desired. Scale everything, except the DC
336 # component and the NFFT/2 component:
338 # if we have a even number of frequencies, don't scale NFFT/2
339 if not NFFT % 2:
340 slc = slice(1, -1, None)
341 # if we have an odd number, just don't scale DC
342 else:
343 slc = slice(1, None, None)
345 result[slc] *= scaling_factor
347 # MATLAB divides by the sampling frequency so that density function
348 # has units of dB/Hz and can be integrated by the plotted frequency
349 # values. Perform the same scaling here.
350 if scale_by_freq:
351 result /= Fs
352 # Scale the spectrum by the norm of the window to compensate for
353 # windowing loss; see Bendat & Piersol Sec 11.5.2.
354 result /= (window**2).sum()
355 else:
356 # In this case, preserve power in the segment, not amplitude
357 result /= window.sum()**2
359 t = np.arange(NFFT/2, len(x) - NFFT/2 + 1, NFFT - noverlap)/Fs
361 if sides == 'twosided':
362 # center the frequency range at zero
363 freqs = np.roll(freqs, -freqcenter, axis=0)
364 result = np.roll(result, -freqcenter, axis=0)
365 elif not pad_to % 2:
366 # get the last value correctly, it is negative otherwise
367 freqs[-1] *= -1
369 # we unwrap the phase here to handle the onesided vs. twosided case
370 if mode == 'phase':
371 result = np.unwrap(result, axis=0)
373 return result, freqs, t
376def _single_spectrum_helper(
377 mode, x, Fs=None, window=None, pad_to=None, sides=None):
378 """
379 Private helper implementing the commonality between the complex, magnitude,
380 angle, and phase spectrums.
381 """
382 _api.check_in_list(['complex', 'magnitude', 'angle', 'phase'], mode=mode)
384 if pad_to is None:
385 pad_to = len(x)
387 spec, freqs, _ = _spectral_helper(x=x, y=None, NFFT=len(x), Fs=Fs,
388 detrend_func=detrend_none, window=window,
389 noverlap=0, pad_to=pad_to,
390 sides=sides,
391 scale_by_freq=False,
392 mode=mode)
393 if mode != 'complex':
394 spec = spec.real
396 if spec.ndim == 2 and spec.shape[1] == 1:
397 spec = spec[:, 0]
399 return spec, freqs
402# Split out these keyword docs so that they can be used elsewhere
403_docstring.interpd.update(
404 Spectral="""\
405Fs : float, default: 2
406 The sampling frequency (samples per time unit). It is used to calculate
407 the Fourier frequencies, *freqs*, in cycles per time unit.
409window : callable or ndarray, default: `.window_hanning`
410 A function or a vector of length *NFFT*. To create window vectors see
411 `.window_hanning`, `.window_none`, `numpy.blackman`, `numpy.hamming`,
412 `numpy.bartlett`, `scipy.signal`, `scipy.signal.get_window`, etc. If a
413 function is passed as the argument, it must take a data segment as an
414 argument and return the windowed version of the segment.
416sides : {'default', 'onesided', 'twosided'}, optional
417 Which sides of the spectrum to return. 'default' is one-sided for real
418 data and two-sided for complex data. 'onesided' forces the return of a
419 one-sided spectrum, while 'twosided' forces two-sided.""",
421 Single_Spectrum="""\
422pad_to : int, optional
423 The number of points to which the data segment is padded when performing
424 the FFT. While not increasing the actual resolution of the spectrum (the
425 minimum distance between resolvable peaks), this can give more points in
426 the plot, allowing for more detail. This corresponds to the *n* parameter
427 in the call to `~numpy.fft.fft`. The default is None, which sets *pad_to*
428 equal to the length of the input signal (i.e. no padding).""",
430 PSD="""\
431pad_to : int, optional
432 The number of points to which the data segment is padded when performing
433 the FFT. This can be different from *NFFT*, which specifies the number
434 of data points used. While not increasing the actual resolution of the
435 spectrum (the minimum distance between resolvable peaks), this can give
436 more points in the plot, allowing for more detail. This corresponds to
437 the *n* parameter in the call to `~numpy.fft.fft`. The default is None,
438 which sets *pad_to* equal to *NFFT*
440NFFT : int, default: 256
441 The number of data points used in each block for the FFT. A power 2 is
442 most efficient. This should *NOT* be used to get zero padding, or the
443 scaling of the result will be incorrect; use *pad_to* for this instead.
445detrend : {'none', 'mean', 'linear'} or callable, default: 'none'
446 The function applied to each segment before fft-ing, designed to remove
447 the mean or linear trend. Unlike in MATLAB, where the *detrend* parameter
448 is a vector, in Matplotlib it is a function. The :mod:`~matplotlib.mlab`
449 module defines `.detrend_none`, `.detrend_mean`, and `.detrend_linear`,
450 but you can use a custom function as well. You can also use a string to
451 choose one of the functions: 'none' calls `.detrend_none`. 'mean' calls
452 `.detrend_mean`. 'linear' calls `.detrend_linear`.
454scale_by_freq : bool, default: True
455 Whether the resulting density values should be scaled by the scaling
456 frequency, which gives density in units of 1/Hz. This allows for
457 integration over the returned frequency values. The default is True for
458 MATLAB compatibility.""")
461@_docstring.dedent_interpd
462def psd(x, NFFT=None, Fs=None, detrend=None, window=None,
463 noverlap=None, pad_to=None, sides=None, scale_by_freq=None):
464 r"""
465 Compute the power spectral density.
467 The power spectral density :math:`P_{xx}` by Welch's average
468 periodogram method. The vector *x* is divided into *NFFT* length
469 segments. Each segment is detrended by function *detrend* and
470 windowed by function *window*. *noverlap* gives the length of
471 the overlap between segments. The :math:`|\mathrm{fft}(i)|^2`
472 of each segment :math:`i` are averaged to compute :math:`P_{xx}`.
474 If len(*x*) < *NFFT*, it will be zero padded to *NFFT*.
476 Parameters
477 ----------
478 x : 1-D array or sequence
479 Array or sequence containing the data
481 %(Spectral)s
483 %(PSD)s
485 noverlap : int, default: 0 (no overlap)
486 The number of points of overlap between segments.
488 Returns
489 -------
490 Pxx : 1-D array
491 The values for the power spectrum :math:`P_{xx}` (real valued)
493 freqs : 1-D array
494 The frequencies corresponding to the elements in *Pxx*
496 References
497 ----------
498 Bendat & Piersol -- Random Data: Analysis and Measurement Procedures, John
499 Wiley & Sons (1986)
501 See Also
502 --------
503 specgram
504 `specgram` differs in the default overlap; in not returning the mean of
505 the segment periodograms; and in returning the times of the segments.
507 magnitude_spectrum : returns the magnitude spectrum.
509 csd : returns the spectral density between two signals.
510 """
511 Pxx, freqs = csd(x=x, y=None, NFFT=NFFT, Fs=Fs, detrend=detrend,
512 window=window, noverlap=noverlap, pad_to=pad_to,
513 sides=sides, scale_by_freq=scale_by_freq)
514 return Pxx.real, freqs
517@_docstring.dedent_interpd
518def csd(x, y, NFFT=None, Fs=None, detrend=None, window=None,
519 noverlap=None, pad_to=None, sides=None, scale_by_freq=None):
520 """
521 Compute the cross-spectral density.
523 The cross spectral density :math:`P_{xy}` by Welch's average
524 periodogram method. The vectors *x* and *y* are divided into
525 *NFFT* length segments. Each segment is detrended by function
526 *detrend* and windowed by function *window*. *noverlap* gives
527 the length of the overlap between segments. The product of
528 the direct FFTs of *x* and *y* are averaged over each segment
529 to compute :math:`P_{xy}`, with a scaling to correct for power
530 loss due to windowing.
532 If len(*x*) < *NFFT* or len(*y*) < *NFFT*, they will be zero
533 padded to *NFFT*.
535 Parameters
536 ----------
537 x, y : 1-D arrays or sequences
538 Arrays or sequences containing the data
540 %(Spectral)s
542 %(PSD)s
544 noverlap : int, default: 0 (no overlap)
545 The number of points of overlap between segments.
547 Returns
548 -------
549 Pxy : 1-D array
550 The values for the cross spectrum :math:`P_{xy}` before scaling (real
551 valued)
553 freqs : 1-D array
554 The frequencies corresponding to the elements in *Pxy*
556 References
557 ----------
558 Bendat & Piersol -- Random Data: Analysis and Measurement Procedures, John
559 Wiley & Sons (1986)
561 See Also
562 --------
563 psd : equivalent to setting ``y = x``.
564 """
565 if NFFT is None:
566 NFFT = 256
567 Pxy, freqs, _ = _spectral_helper(x=x, y=y, NFFT=NFFT, Fs=Fs,
568 detrend_func=detrend, window=window,
569 noverlap=noverlap, pad_to=pad_to,
570 sides=sides, scale_by_freq=scale_by_freq,
571 mode='psd')
573 if Pxy.ndim == 2:
574 if Pxy.shape[1] > 1:
575 Pxy = Pxy.mean(axis=1)
576 else:
577 Pxy = Pxy[:, 0]
578 return Pxy, freqs
581_single_spectrum_docs = """\
582Compute the {quantity} of *x*.
583Data is padded to a length of *pad_to* and the windowing function *window* is
584applied to the signal.
586Parameters
587----------
588x : 1-D array or sequence
589 Array or sequence containing the data
591{Spectral}
593{Single_Spectrum}
595Returns
596-------
597spectrum : 1-D array
598 The {quantity}.
599freqs : 1-D array
600 The frequencies corresponding to the elements in *spectrum*.
602See Also
603--------
604psd
605 Returns the power spectral density.
606complex_spectrum
607 Returns the complex-valued frequency spectrum.
608magnitude_spectrum
609 Returns the absolute value of the `complex_spectrum`.
610angle_spectrum
611 Returns the angle of the `complex_spectrum`.
612phase_spectrum
613 Returns the phase (unwrapped angle) of the `complex_spectrum`.
614specgram
615 Can return the complex spectrum of segments within the signal.
616"""
619complex_spectrum = functools.partial(_single_spectrum_helper, "complex")
620complex_spectrum.__doc__ = _single_spectrum_docs.format(
621 quantity="complex-valued frequency spectrum",
622 **_docstring.interpd.params)
623magnitude_spectrum = functools.partial(_single_spectrum_helper, "magnitude")
624magnitude_spectrum.__doc__ = _single_spectrum_docs.format(
625 quantity="magnitude (absolute value) of the frequency spectrum",
626 **_docstring.interpd.params)
627angle_spectrum = functools.partial(_single_spectrum_helper, "angle")
628angle_spectrum.__doc__ = _single_spectrum_docs.format(
629 quantity="angle of the frequency spectrum (wrapped phase spectrum)",
630 **_docstring.interpd.params)
631phase_spectrum = functools.partial(_single_spectrum_helper, "phase")
632phase_spectrum.__doc__ = _single_spectrum_docs.format(
633 quantity="phase of the frequency spectrum (unwrapped phase spectrum)",
634 **_docstring.interpd.params)
637@_docstring.dedent_interpd
638def specgram(x, NFFT=None, Fs=None, detrend=None, window=None,
639 noverlap=None, pad_to=None, sides=None, scale_by_freq=None,
640 mode=None):
641 """
642 Compute a spectrogram.
644 Compute and plot a spectrogram of data in *x*. Data are split into
645 *NFFT* length segments and the spectrum of each section is
646 computed. The windowing function *window* is applied to each
647 segment, and the amount of overlap of each segment is
648 specified with *noverlap*.
650 Parameters
651 ----------
652 x : array-like
653 1-D array or sequence.
655 %(Spectral)s
657 %(PSD)s
659 noverlap : int, default: 128
660 The number of points of overlap between blocks.
661 mode : str, default: 'psd'
662 What sort of spectrum to use:
663 'psd'
664 Returns the power spectral density.
665 'complex'
666 Returns the complex-valued frequency spectrum.
667 'magnitude'
668 Returns the magnitude spectrum.
669 'angle'
670 Returns the phase spectrum without unwrapping.
671 'phase'
672 Returns the phase spectrum with unwrapping.
674 Returns
675 -------
676 spectrum : array-like
677 2D array, columns are the periodograms of successive segments.
679 freqs : array-like
680 1-D array, frequencies corresponding to the rows in *spectrum*.
682 t : array-like
683 1-D array, the times corresponding to midpoints of segments
684 (i.e the columns in *spectrum*).
686 See Also
687 --------
688 psd : differs in the overlap and in the return values.
689 complex_spectrum : similar, but with complex valued frequencies.
690 magnitude_spectrum : similar single segment when *mode* is 'magnitude'.
691 angle_spectrum : similar to single segment when *mode* is 'angle'.
692 phase_spectrum : similar to single segment when *mode* is 'phase'.
694 Notes
695 -----
696 *detrend* and *scale_by_freq* only apply when *mode* is set to 'psd'.
698 """
699 if noverlap is None:
700 noverlap = 128 # default in _spectral_helper() is noverlap = 0
701 if NFFT is None:
702 NFFT = 256 # same default as in _spectral_helper()
703 if len(x) <= NFFT:
704 _api.warn_external("Only one segment is calculated since parameter "
705 f"NFFT (={NFFT}) >= signal length (={len(x)}).")
707 spec, freqs, t = _spectral_helper(x=x, y=None, NFFT=NFFT, Fs=Fs,
708 detrend_func=detrend, window=window,
709 noverlap=noverlap, pad_to=pad_to,
710 sides=sides,
711 scale_by_freq=scale_by_freq,
712 mode=mode)
714 if mode != 'complex':
715 spec = spec.real # Needed since helper implements generically
717 return spec, freqs, t
720@_docstring.dedent_interpd
721def cohere(x, y, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning,
722 noverlap=0, pad_to=None, sides='default', scale_by_freq=None):
723 r"""
724 The coherence between *x* and *y*. Coherence is the normalized
725 cross spectral density:
727 .. math::
729 C_{xy} = \frac{|P_{xy}|^2}{P_{xx}P_{yy}}
731 Parameters
732 ----------
733 x, y
734 Array or sequence containing the data
736 %(Spectral)s
738 %(PSD)s
740 noverlap : int, default: 0 (no overlap)
741 The number of points of overlap between segments.
743 Returns
744 -------
745 Cxy : 1-D array
746 The coherence vector.
747 freqs : 1-D array
748 The frequencies for the elements in *Cxy*.
750 See Also
751 --------
752 :func:`psd`, :func:`csd` :
753 For information about the methods used to compute :math:`P_{xy}`,
754 :math:`P_{xx}` and :math:`P_{yy}`.
755 """
756 if len(x) < 2 * NFFT:
757 raise ValueError(
758 "Coherence is calculated by averaging over *NFFT* length "
759 "segments. Your signal is too short for your choice of *NFFT*.")
760 Pxx, f = psd(x, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
761 scale_by_freq)
762 Pyy, f = psd(y, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
763 scale_by_freq)
764 Pxy, f = csd(x, y, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
765 scale_by_freq)
766 Cxy = np.abs(Pxy) ** 2 / (Pxx * Pyy)
767 return Cxy, f
770class GaussianKDE:
771 """
772 Representation of a kernel-density estimate using Gaussian kernels.
774 Parameters
775 ----------
776 dataset : array-like
777 Datapoints to estimate from. In case of univariate data this is a 1-D
778 array, otherwise a 2D array with shape (# of dims, # of data).
779 bw_method : {'scott', 'silverman'} or float or callable, optional
780 The method used to calculate the estimator bandwidth. If a
781 float, this will be used directly as `kde.factor`. If a
782 callable, it should take a `GaussianKDE` instance as only
783 parameter and return a float. If None (default), 'scott' is used.
785 Attributes
786 ----------
787 dataset : ndarray
788 The dataset passed to the constructor.
789 dim : int
790 Number of dimensions.
791 num_dp : int
792 Number of datapoints.
793 factor : float
794 The bandwidth factor, obtained from `kde.covariance_factor`, with which
795 the covariance matrix is multiplied.
796 covariance : ndarray
797 The covariance matrix of *dataset*, scaled by the calculated bandwidth
798 (`kde.factor`).
799 inv_cov : ndarray
800 The inverse of *covariance*.
802 Methods
803 -------
804 kde.evaluate(points) : ndarray
805 Evaluate the estimated pdf on a provided set of points.
806 kde(points) : ndarray
807 Same as kde.evaluate(points)
808 """
810 # This implementation with minor modification was too good to pass up.
811 # from scipy: https://github.com/scipy/scipy/blob/master/scipy/stats/kde.py
813 def __init__(self, dataset, bw_method=None):
814 self.dataset = np.atleast_2d(dataset)
815 if not np.array(self.dataset).size > 1:
816 raise ValueError("`dataset` input should have multiple elements.")
818 self.dim, self.num_dp = np.array(self.dataset).shape
820 if bw_method is None:
821 pass
822 elif cbook._str_equal(bw_method, 'scott'):
823 self.covariance_factor = self.scotts_factor
824 elif cbook._str_equal(bw_method, 'silverman'):
825 self.covariance_factor = self.silverman_factor
826 elif isinstance(bw_method, Number):
827 self._bw_method = 'use constant'
828 self.covariance_factor = lambda: bw_method
829 elif callable(bw_method):
830 self._bw_method = bw_method
831 self.covariance_factor = lambda: self._bw_method(self)
832 else:
833 raise ValueError("`bw_method` should be 'scott', 'silverman', a "
834 "scalar or a callable")
836 # Computes the covariance matrix for each Gaussian kernel using
837 # covariance_factor().
839 self.factor = self.covariance_factor()
840 # Cache covariance and inverse covariance of the data
841 if not hasattr(self, '_data_inv_cov'):
842 self.data_covariance = np.atleast_2d(
843 np.cov(
844 self.dataset,
845 rowvar=1,
846 bias=False))
847 self.data_inv_cov = np.linalg.inv(self.data_covariance)
849 self.covariance = self.data_covariance * self.factor ** 2
850 self.inv_cov = self.data_inv_cov / self.factor ** 2
851 self.norm_factor = (np.sqrt(np.linalg.det(2 * np.pi * self.covariance))
852 * self.num_dp)
854 def scotts_factor(self):
855 return np.power(self.num_dp, -1. / (self.dim + 4))
857 def silverman_factor(self):
858 return np.power(
859 self.num_dp * (self.dim + 2.0) / 4.0, -1. / (self.dim + 4))
861 # Default method to calculate bandwidth, can be overwritten by subclass
862 covariance_factor = scotts_factor
864 def evaluate(self, points):
865 """
866 Evaluate the estimated pdf on a set of points.
868 Parameters
869 ----------
870 points : (# of dimensions, # of points)-array
871 Alternatively, a (# of dimensions,) vector can be passed in and
872 treated as a single point.
874 Returns
875 -------
876 (# of points,)-array
877 The values at each point.
879 Raises
880 ------
881 ValueError : if the dimensionality of the input points is different
882 than the dimensionality of the KDE.
884 """
885 points = np.atleast_2d(points)
887 dim, num_m = np.array(points).shape
888 if dim != self.dim:
889 raise ValueError(f"points have dimension {dim}, dataset has "
890 f"dimension {self.dim}")
892 result = np.zeros(num_m)
894 if num_m >= self.num_dp:
895 # there are more points than data, so loop over data
896 for i in range(self.num_dp):
897 diff = self.dataset[:, i, np.newaxis] - points
898 tdiff = np.dot(self.inv_cov, diff)
899 energy = np.sum(diff * tdiff, axis=0) / 2.0
900 result = result + np.exp(-energy)
901 else:
902 # loop over points
903 for i in range(num_m):
904 diff = self.dataset - points[:, i, np.newaxis]
905 tdiff = np.dot(self.inv_cov, diff)
906 energy = np.sum(diff * tdiff, axis=0) / 2.0
907 result[i] = np.sum(np.exp(-energy), axis=0)
909 result = result / self.norm_factor
911 return result
913 __call__ = evaluate