Coverage for /usr/lib/python3/dist-packages/matplotlib/mlab.py: 15%
275 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
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.
49`stride_windows`
50 Get all windows in an array in a memory-efficient manner
51"""
53import functools
54from numbers import Number
56import numpy as np
58from matplotlib import _api, _docstring, cbook
61def window_hanning(x):
62 """
63 Return *x* times the Hanning (or Hann) window of len(*x*).
65 See Also
66 --------
67 window_none : Another window algorithm.
68 """
69 return np.hanning(len(x))*x
72def window_none(x):
73 """
74 No window function; simply return *x*.
76 See Also
77 --------
78 window_hanning : Another window algorithm.
79 """
80 return x
83def detrend(x, key=None, axis=None):
84 """
85 Return *x* with its trend removed.
87 Parameters
88 ----------
89 x : array or sequence
90 Array or sequence containing the data.
92 key : {'default', 'constant', 'mean', 'linear', 'none'} or function
93 The detrending algorithm to use. 'default', 'mean', and 'constant' are
94 the same as `detrend_mean`. 'linear' is the same as `detrend_linear`.
95 'none' is the same as `detrend_none`. The default is 'mean'. See the
96 corresponding functions for more details regarding the algorithms. Can
97 also be a function that carries out the detrend operation.
99 axis : int
100 The axis along which to do the detrending.
102 See Also
103 --------
104 detrend_mean : Implementation of the 'mean' algorithm.
105 detrend_linear : Implementation of the 'linear' algorithm.
106 detrend_none : Implementation of the 'none' algorithm.
107 """
108 if key is None or key in ['constant', 'mean', 'default']:
109 return detrend(x, key=detrend_mean, axis=axis)
110 elif key == 'linear':
111 return detrend(x, key=detrend_linear, axis=axis)
112 elif key == 'none':
113 return detrend(x, key=detrend_none, axis=axis)
114 elif callable(key):
115 x = np.asarray(x)
116 if axis is not None and axis + 1 > x.ndim:
117 raise ValueError(f'axis(={axis}) out of bounds')
118 if (axis is None and x.ndim == 0) or (not axis and x.ndim == 1):
119 return key(x)
120 # try to use the 'axis' argument if the function supports it,
121 # otherwise use apply_along_axis to do it
122 try:
123 return key(x, axis=axis)
124 except TypeError:
125 return np.apply_along_axis(key, axis=axis, arr=x)
126 else:
127 raise ValueError(
128 f"Unknown value for key: {key!r}, must be one of: 'default', "
129 f"'constant', 'mean', 'linear', or a function")
132def detrend_mean(x, axis=None):
133 """
134 Return *x* minus the mean(*x*).
136 Parameters
137 ----------
138 x : array or sequence
139 Array or sequence containing the data
140 Can have any dimensionality
142 axis : int
143 The axis along which to take the mean. See `numpy.mean` for a
144 description of this argument.
146 See Also
147 --------
148 detrend_linear : Another detrend algorithm.
149 detrend_none : Another detrend algorithm.
150 detrend : A wrapper around all the detrend algorithms.
151 """
152 x = np.asarray(x)
154 if axis is not None and axis+1 > x.ndim:
155 raise ValueError('axis(=%s) out of bounds' % axis)
157 return x - x.mean(axis, keepdims=True)
160def detrend_none(x, axis=None):
161 """
162 Return *x*: no detrending.
164 Parameters
165 ----------
166 x : any object
167 An object containing the data
169 axis : int
170 This parameter is ignored.
171 It is included for compatibility with detrend_mean
173 See Also
174 --------
175 detrend_mean : Another detrend algorithm.
176 detrend_linear : Another detrend algorithm.
177 detrend : A wrapper around all the detrend algorithms.
178 """
179 return x
182def detrend_linear(y):
183 """
184 Return *x* minus best fit line; 'linear' detrending.
186 Parameters
187 ----------
188 y : 0-D or 1-D array or sequence
189 Array or sequence containing the data
191 See Also
192 --------
193 detrend_mean : Another detrend algorithm.
194 detrend_none : Another detrend algorithm.
195 detrend : A wrapper around all the detrend algorithms.
196 """
197 # This is faster than an algorithm based on linalg.lstsq.
198 y = np.asarray(y)
200 if y.ndim > 1:
201 raise ValueError('y cannot have ndim > 1')
203 # short-circuit 0-D array.
204 if not y.ndim:
205 return np.array(0., dtype=y.dtype)
207 x = np.arange(y.size, dtype=float)
209 C = np.cov(x, y, bias=1)
210 b = C[0, 1]/C[0, 0]
212 a = y.mean() - b*x.mean()
213 return y - (b*x + a)
216@_api.deprecated("3.6")
217def stride_windows(x, n, noverlap=None, axis=0):
218 """
219 Get all windows of *x* with length *n* as a single array,
220 using strides to avoid data duplication.
222 .. warning::
224 It is not safe to write to the output array. Multiple
225 elements may point to the same piece of memory,
226 so modifying one value may change others.
228 Parameters
229 ----------
230 x : 1D array or sequence
231 Array or sequence containing the data.
232 n : int
233 The number of data points in each window.
234 noverlap : int, default: 0 (no overlap)
235 The overlap between adjacent windows.
236 axis : int
237 The axis along which the windows will run.
239 References
240 ----------
241 `stackoverflow: Rolling window for 1D arrays in Numpy?
242 <https://stackoverflow.com/a/6811241>`_
243 `stackoverflow: Using strides for an efficient moving average filter
244 <https://stackoverflow.com/a/4947453>`_
245 """
246 if noverlap is None:
247 noverlap = 0
248 if np.ndim(x) != 1:
249 raise ValueError('only 1-dimensional arrays can be used')
250 return _stride_windows(x, n, noverlap, axis)
253def _stride_windows(x, n, noverlap=0, axis=0):
254 # np>=1.20 provides sliding_window_view, and we only ever use axis=0.
255 if hasattr(np.lib.stride_tricks, "sliding_window_view") and axis == 0:
256 if noverlap >= n:
257 raise ValueError('noverlap must be less than n')
258 return np.lib.stride_tricks.sliding_window_view(
259 x, n, axis=0)[::n - noverlap].T
261 if noverlap >= n:
262 raise ValueError('noverlap must be less than n')
263 if n < 1:
264 raise ValueError('n cannot be less than 1')
266 x = np.asarray(x)
268 if n == 1 and noverlap == 0:
269 if axis == 0:
270 return x[np.newaxis]
271 else:
272 return x[np.newaxis].T
273 if n > x.size:
274 raise ValueError('n cannot be greater than the length of x')
276 # np.lib.stride_tricks.as_strided easily leads to memory corruption for
277 # non integer shape and strides, i.e. noverlap or n. See #3845.
278 noverlap = int(noverlap)
279 n = int(n)
281 step = n - noverlap
282 if axis == 0:
283 shape = (n, (x.shape[-1]-noverlap)//step)
284 strides = (x.strides[0], step*x.strides[0])
285 else:
286 shape = ((x.shape[-1]-noverlap)//step, n)
287 strides = (step*x.strides[0], x.strides[0])
288 return np.lib.stride_tricks.as_strided(x, shape=shape, strides=strides)
291def _spectral_helper(x, y=None, NFFT=None, Fs=None, detrend_func=None,
292 window=None, noverlap=None, pad_to=None,
293 sides=None, scale_by_freq=None, mode=None):
294 """
295 Private helper implementing the common parts between the psd, csd,
296 spectrogram and complex, magnitude, angle, and phase spectrums.
297 """
298 if y is None:
299 # if y is None use x for y
300 same_data = True
301 else:
302 # The checks for if y is x are so that we can use the same function to
303 # implement the core of psd(), csd(), and spectrogram() without doing
304 # extra calculations. We return the unaveraged Pxy, freqs, and t.
305 same_data = y is x
307 if Fs is None:
308 Fs = 2
309 if noverlap is None:
310 noverlap = 0
311 if detrend_func is None:
312 detrend_func = detrend_none
313 if window is None:
314 window = window_hanning
316 # if NFFT is set to None use the whole signal
317 if NFFT is None:
318 NFFT = 256
320 if mode is None or mode == 'default':
321 mode = 'psd'
322 _api.check_in_list(
323 ['default', 'psd', 'complex', 'magnitude', 'angle', 'phase'],
324 mode=mode)
326 if not same_data and mode != 'psd':
327 raise ValueError("x and y must be equal if mode is not 'psd'")
329 # Make sure we're dealing with a numpy array. If y and x were the same
330 # object to start with, keep them that way
331 x = np.asarray(x)
332 if not same_data:
333 y = np.asarray(y)
335 if sides is None or sides == 'default':
336 if np.iscomplexobj(x):
337 sides = 'twosided'
338 else:
339 sides = 'onesided'
340 _api.check_in_list(['default', 'onesided', 'twosided'], sides=sides)
342 # zero pad x and y up to NFFT if they are shorter than NFFT
343 if len(x) < NFFT:
344 n = len(x)
345 x = np.resize(x, NFFT)
346 x[n:] = 0
348 if not same_data and len(y) < NFFT:
349 n = len(y)
350 y = np.resize(y, NFFT)
351 y[n:] = 0
353 if pad_to is None:
354 pad_to = NFFT
356 if mode != 'psd':
357 scale_by_freq = False
358 elif scale_by_freq is None:
359 scale_by_freq = True
361 # For real x, ignore the negative frequencies unless told otherwise
362 if sides == 'twosided':
363 numFreqs = pad_to
364 if pad_to % 2:
365 freqcenter = (pad_to - 1)//2 + 1
366 else:
367 freqcenter = pad_to//2
368 scaling_factor = 1.
369 elif sides == 'onesided':
370 if pad_to % 2:
371 numFreqs = (pad_to + 1)//2
372 else:
373 numFreqs = pad_to//2 + 1
374 scaling_factor = 2.
376 if not np.iterable(window):
377 window = window(np.ones(NFFT, x.dtype))
378 if len(window) != NFFT:
379 raise ValueError(
380 "The window length must match the data's first dimension")
382 result = _stride_windows(x, NFFT, noverlap)
383 result = detrend(result, detrend_func, axis=0)
384 result = result * window.reshape((-1, 1))
385 result = np.fft.fft(result, n=pad_to, axis=0)[:numFreqs, :]
386 freqs = np.fft.fftfreq(pad_to, 1/Fs)[:numFreqs]
388 if not same_data:
389 # if same_data is False, mode must be 'psd'
390 resultY = _stride_windows(y, NFFT, noverlap)
391 resultY = detrend(resultY, detrend_func, axis=0)
392 resultY = resultY * window.reshape((-1, 1))
393 resultY = np.fft.fft(resultY, n=pad_to, axis=0)[:numFreqs, :]
394 result = np.conj(result) * resultY
395 elif mode == 'psd':
396 result = np.conj(result) * result
397 elif mode == 'magnitude':
398 result = np.abs(result) / np.abs(window).sum()
399 elif mode == 'angle' or mode == 'phase':
400 # we unwrap the phase later to handle the onesided vs. twosided case
401 result = np.angle(result)
402 elif mode == 'complex':
403 result /= np.abs(window).sum()
405 if mode == 'psd':
407 # Also include scaling factors for one-sided densities and dividing by
408 # the sampling frequency, if desired. Scale everything, except the DC
409 # component and the NFFT/2 component:
411 # if we have a even number of frequencies, don't scale NFFT/2
412 if not NFFT % 2:
413 slc = slice(1, -1, None)
414 # if we have an odd number, just don't scale DC
415 else:
416 slc = slice(1, None, None)
418 result[slc] *= scaling_factor
420 # MATLAB divides by the sampling frequency so that density function
421 # has units of dB/Hz and can be integrated by the plotted frequency
422 # values. Perform the same scaling here.
423 if scale_by_freq:
424 result /= Fs
425 # Scale the spectrum by the norm of the window to compensate for
426 # windowing loss; see Bendat & Piersol Sec 11.5.2.
427 result /= (np.abs(window)**2).sum()
428 else:
429 # In this case, preserve power in the segment, not amplitude
430 result /= np.abs(window).sum()**2
432 t = np.arange(NFFT/2, len(x) - NFFT/2 + 1, NFFT - noverlap)/Fs
434 if sides == 'twosided':
435 # center the frequency range at zero
436 freqs = np.roll(freqs, -freqcenter, axis=0)
437 result = np.roll(result, -freqcenter, axis=0)
438 elif not pad_to % 2:
439 # get the last value correctly, it is negative otherwise
440 freqs[-1] *= -1
442 # we unwrap the phase here to handle the onesided vs. twosided case
443 if mode == 'phase':
444 result = np.unwrap(result, axis=0)
446 return result, freqs, t
449def _single_spectrum_helper(
450 mode, x, Fs=None, window=None, pad_to=None, sides=None):
451 """
452 Private helper implementing the commonality between the complex, magnitude,
453 angle, and phase spectrums.
454 """
455 _api.check_in_list(['complex', 'magnitude', 'angle', 'phase'], mode=mode)
457 if pad_to is None:
458 pad_to = len(x)
460 spec, freqs, _ = _spectral_helper(x=x, y=None, NFFT=len(x), Fs=Fs,
461 detrend_func=detrend_none, window=window,
462 noverlap=0, pad_to=pad_to,
463 sides=sides,
464 scale_by_freq=False,
465 mode=mode)
466 if mode != 'complex':
467 spec = spec.real
469 if spec.ndim == 2 and spec.shape[1] == 1:
470 spec = spec[:, 0]
472 return spec, freqs
475# Split out these keyword docs so that they can be used elsewhere
476_docstring.interpd.update(
477 Spectral="""\
478Fs : float, default: 2
479 The sampling frequency (samples per time unit). It is used to calculate
480 the Fourier frequencies, *freqs*, in cycles per time unit.
482window : callable or ndarray, default: `.window_hanning`
483 A function or a vector of length *NFFT*. To create window vectors see
484 `.window_hanning`, `.window_none`, `numpy.blackman`, `numpy.hamming`,
485 `numpy.bartlett`, `scipy.signal`, `scipy.signal.get_window`, etc. If a
486 function is passed as the argument, it must take a data segment as an
487 argument and return the windowed version of the segment.
489sides : {'default', 'onesided', 'twosided'}, optional
490 Which sides of the spectrum to return. 'default' is one-sided for real
491 data and two-sided for complex data. 'onesided' forces the return of a
492 one-sided spectrum, while 'twosided' forces two-sided.""",
494 Single_Spectrum="""\
495pad_to : int, optional
496 The number of points to which the data segment is padded when performing
497 the FFT. While not increasing the actual resolution of the spectrum (the
498 minimum distance between resolvable peaks), this can give more points in
499 the plot, allowing for more detail. This corresponds to the *n* parameter
500 in the call to `~numpy.fft.fft`. The default is None, which sets *pad_to*
501 equal to the length of the input signal (i.e. no padding).""",
503 PSD="""\
504pad_to : int, optional
505 The number of points to which the data segment is padded when performing
506 the FFT. This can be different from *NFFT*, which specifies the number
507 of data points used. While not increasing the actual resolution of the
508 spectrum (the minimum distance between resolvable peaks), this can give
509 more points in the plot, allowing for more detail. This corresponds to
510 the *n* parameter in the call to `~numpy.fft.fft`. The default is None,
511 which sets *pad_to* equal to *NFFT*
513NFFT : int, default: 256
514 The number of data points used in each block for the FFT. A power 2 is
515 most efficient. This should *NOT* be used to get zero padding, or the
516 scaling of the result will be incorrect; use *pad_to* for this instead.
518detrend : {'none', 'mean', 'linear'} or callable, default: 'none'
519 The function applied to each segment before fft-ing, designed to remove
520 the mean or linear trend. Unlike in MATLAB, where the *detrend* parameter
521 is a vector, in Matplotlib it is a function. The :mod:`~matplotlib.mlab`
522 module defines `.detrend_none`, `.detrend_mean`, and `.detrend_linear`,
523 but you can use a custom function as well. You can also use a string to
524 choose one of the functions: 'none' calls `.detrend_none`. 'mean' calls
525 `.detrend_mean`. 'linear' calls `.detrend_linear`.
527scale_by_freq : bool, default: True
528 Whether the resulting density values should be scaled by the scaling
529 frequency, which gives density in units of 1/Hz. This allows for
530 integration over the returned frequency values. The default is True for
531 MATLAB compatibility.""")
534@_docstring.dedent_interpd
535def psd(x, NFFT=None, Fs=None, detrend=None, window=None,
536 noverlap=None, pad_to=None, sides=None, scale_by_freq=None):
537 r"""
538 Compute the power spectral density.
540 The power spectral density :math:`P_{xx}` by Welch's average
541 periodogram method. The vector *x* is divided into *NFFT* length
542 segments. Each segment is detrended by function *detrend* and
543 windowed by function *window*. *noverlap* gives the length of
544 the overlap between segments. The :math:`|\mathrm{fft}(i)|^2`
545 of each segment :math:`i` are averaged to compute :math:`P_{xx}`.
547 If len(*x*) < *NFFT*, it will be zero padded to *NFFT*.
549 Parameters
550 ----------
551 x : 1-D array or sequence
552 Array or sequence containing the data
554 %(Spectral)s
556 %(PSD)s
558 noverlap : int, default: 0 (no overlap)
559 The number of points of overlap between segments.
561 Returns
562 -------
563 Pxx : 1-D array
564 The values for the power spectrum :math:`P_{xx}` (real valued)
566 freqs : 1-D array
567 The frequencies corresponding to the elements in *Pxx*
569 References
570 ----------
571 Bendat & Piersol -- Random Data: Analysis and Measurement Procedures, John
572 Wiley & Sons (1986)
574 See Also
575 --------
576 specgram
577 `specgram` differs in the default overlap; in not returning the mean of
578 the segment periodograms; and in returning the times of the segments.
580 magnitude_spectrum : returns the magnitude spectrum.
582 csd : returns the spectral density between two signals.
583 """
584 Pxx, freqs = csd(x=x, y=None, NFFT=NFFT, Fs=Fs, detrend=detrend,
585 window=window, noverlap=noverlap, pad_to=pad_to,
586 sides=sides, scale_by_freq=scale_by_freq)
587 return Pxx.real, freqs
590@_docstring.dedent_interpd
591def csd(x, y, NFFT=None, Fs=None, detrend=None, window=None,
592 noverlap=None, pad_to=None, sides=None, scale_by_freq=None):
593 """
594 Compute the cross-spectral density.
596 The cross spectral density :math:`P_{xy}` by Welch's average
597 periodogram method. The vectors *x* and *y* are divided into
598 *NFFT* length segments. Each segment is detrended by function
599 *detrend* and windowed by function *window*. *noverlap* gives
600 the length of the overlap between segments. The product of
601 the direct FFTs of *x* and *y* are averaged over each segment
602 to compute :math:`P_{xy}`, with a scaling to correct for power
603 loss due to windowing.
605 If len(*x*) < *NFFT* or len(*y*) < *NFFT*, they will be zero
606 padded to *NFFT*.
608 Parameters
609 ----------
610 x, y : 1-D arrays or sequences
611 Arrays or sequences containing the data
613 %(Spectral)s
615 %(PSD)s
617 noverlap : int, default: 0 (no overlap)
618 The number of points of overlap between segments.
620 Returns
621 -------
622 Pxy : 1-D array
623 The values for the cross spectrum :math:`P_{xy}` before scaling (real
624 valued)
626 freqs : 1-D array
627 The frequencies corresponding to the elements in *Pxy*
629 References
630 ----------
631 Bendat & Piersol -- Random Data: Analysis and Measurement Procedures, John
632 Wiley & Sons (1986)
634 See Also
635 --------
636 psd : equivalent to setting ``y = x``.
637 """
638 if NFFT is None:
639 NFFT = 256
640 Pxy, freqs, _ = _spectral_helper(x=x, y=y, NFFT=NFFT, Fs=Fs,
641 detrend_func=detrend, window=window,
642 noverlap=noverlap, pad_to=pad_to,
643 sides=sides, scale_by_freq=scale_by_freq,
644 mode='psd')
646 if Pxy.ndim == 2:
647 if Pxy.shape[1] > 1:
648 Pxy = Pxy.mean(axis=1)
649 else:
650 Pxy = Pxy[:, 0]
651 return Pxy, freqs
654_single_spectrum_docs = """\
655Compute the {quantity} of *x*.
656Data is padded to a length of *pad_to* and the windowing function *window* is
657applied to the signal.
659Parameters
660----------
661x : 1-D array or sequence
662 Array or sequence containing the data
664{Spectral}
666{Single_Spectrum}
668Returns
669-------
670spectrum : 1-D array
671 The {quantity}.
672freqs : 1-D array
673 The frequencies corresponding to the elements in *spectrum*.
675See Also
676--------
677psd
678 Returns the power spectral density.
679complex_spectrum
680 Returns the complex-valued frequency spectrum.
681magnitude_spectrum
682 Returns the absolute value of the `complex_spectrum`.
683angle_spectrum
684 Returns the angle of the `complex_spectrum`.
685phase_spectrum
686 Returns the phase (unwrapped angle) of the `complex_spectrum`.
687specgram
688 Can return the complex spectrum of segments within the signal.
689"""
692complex_spectrum = functools.partial(_single_spectrum_helper, "complex")
693complex_spectrum.__doc__ = _single_spectrum_docs.format(
694 quantity="complex-valued frequency spectrum",
695 **_docstring.interpd.params)
696magnitude_spectrum = functools.partial(_single_spectrum_helper, "magnitude")
697magnitude_spectrum.__doc__ = _single_spectrum_docs.format(
698 quantity="magnitude (absolute value) of the frequency spectrum",
699 **_docstring.interpd.params)
700angle_spectrum = functools.partial(_single_spectrum_helper, "angle")
701angle_spectrum.__doc__ = _single_spectrum_docs.format(
702 quantity="angle of the frequency spectrum (wrapped phase spectrum)",
703 **_docstring.interpd.params)
704phase_spectrum = functools.partial(_single_spectrum_helper, "phase")
705phase_spectrum.__doc__ = _single_spectrum_docs.format(
706 quantity="phase of the frequency spectrum (unwrapped phase spectrum)",
707 **_docstring.interpd.params)
710@_docstring.dedent_interpd
711def specgram(x, NFFT=None, Fs=None, detrend=None, window=None,
712 noverlap=None, pad_to=None, sides=None, scale_by_freq=None,
713 mode=None):
714 """
715 Compute a spectrogram.
717 Compute and plot a spectrogram of data in *x*. Data are split into
718 *NFFT* length segments and the spectrum of each section is
719 computed. The windowing function *window* is applied to each
720 segment, and the amount of overlap of each segment is
721 specified with *noverlap*.
723 Parameters
724 ----------
725 x : array-like
726 1-D array or sequence.
728 %(Spectral)s
730 %(PSD)s
732 noverlap : int, default: 128
733 The number of points of overlap between blocks.
734 mode : str, default: 'psd'
735 What sort of spectrum to use:
736 'psd'
737 Returns the power spectral density.
738 'complex'
739 Returns the complex-valued frequency spectrum.
740 'magnitude'
741 Returns the magnitude spectrum.
742 'angle'
743 Returns the phase spectrum without unwrapping.
744 'phase'
745 Returns the phase spectrum with unwrapping.
747 Returns
748 -------
749 spectrum : array-like
750 2D array, columns are the periodograms of successive segments.
752 freqs : array-like
753 1-D array, frequencies corresponding to the rows in *spectrum*.
755 t : array-like
756 1-D array, the times corresponding to midpoints of segments
757 (i.e the columns in *spectrum*).
759 See Also
760 --------
761 psd : differs in the overlap and in the return values.
762 complex_spectrum : similar, but with complex valued frequencies.
763 magnitude_spectrum : similar single segment when *mode* is 'magnitude'.
764 angle_spectrum : similar to single segment when *mode* is 'angle'.
765 phase_spectrum : similar to single segment when *mode* is 'phase'.
767 Notes
768 -----
769 *detrend* and *scale_by_freq* only apply when *mode* is set to 'psd'.
771 """
772 if noverlap is None:
773 noverlap = 128 # default in _spectral_helper() is noverlap = 0
774 if NFFT is None:
775 NFFT = 256 # same default as in _spectral_helper()
776 if len(x) <= NFFT:
777 _api.warn_external("Only one segment is calculated since parameter "
778 f"NFFT (={NFFT}) >= signal length (={len(x)}).")
780 spec, freqs, t = _spectral_helper(x=x, y=None, NFFT=NFFT, Fs=Fs,
781 detrend_func=detrend, window=window,
782 noverlap=noverlap, pad_to=pad_to,
783 sides=sides,
784 scale_by_freq=scale_by_freq,
785 mode=mode)
787 if mode != 'complex':
788 spec = spec.real # Needed since helper implements generically
790 return spec, freqs, t
793@_docstring.dedent_interpd
794def cohere(x, y, NFFT=256, Fs=2, detrend=detrend_none, window=window_hanning,
795 noverlap=0, pad_to=None, sides='default', scale_by_freq=None):
796 r"""
797 The coherence between *x* and *y*. Coherence is the normalized
798 cross spectral density:
800 .. math::
802 C_{xy} = \frac{|P_{xy}|^2}{P_{xx}P_{yy}}
804 Parameters
805 ----------
806 x, y
807 Array or sequence containing the data
809 %(Spectral)s
811 %(PSD)s
813 noverlap : int, default: 0 (no overlap)
814 The number of points of overlap between segments.
816 Returns
817 -------
818 Cxy : 1-D array
819 The coherence vector.
820 freqs : 1-D array
821 The frequencies for the elements in *Cxy*.
823 See Also
824 --------
825 :func:`psd`, :func:`csd` :
826 For information about the methods used to compute :math:`P_{xy}`,
827 :math:`P_{xx}` and :math:`P_{yy}`.
828 """
829 if len(x) < 2 * NFFT:
830 raise ValueError(
831 "Coherence is calculated by averaging over *NFFT* length "
832 "segments. Your signal is too short for your choice of *NFFT*.")
833 Pxx, f = psd(x, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
834 scale_by_freq)
835 Pyy, f = psd(y, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
836 scale_by_freq)
837 Pxy, f = csd(x, y, NFFT, Fs, detrend, window, noverlap, pad_to, sides,
838 scale_by_freq)
839 Cxy = np.abs(Pxy) ** 2 / (Pxx * Pyy)
840 return Cxy, f
843class GaussianKDE:
844 """
845 Representation of a kernel-density estimate using Gaussian kernels.
847 Parameters
848 ----------
849 dataset : array-like
850 Datapoints to estimate from. In case of univariate data this is a 1-D
851 array, otherwise a 2D array with shape (# of dims, # of data).
852 bw_method : str, scalar or callable, optional
853 The method used to calculate the estimator bandwidth. This can be
854 'scott', 'silverman', a scalar constant or a callable. If a
855 scalar, this will be used directly as `kde.factor`. If a
856 callable, it should take a `GaussianKDE` instance as only
857 parameter and return a scalar. If None (default), 'scott' is used.
859 Attributes
860 ----------
861 dataset : ndarray
862 The dataset passed to the constructor.
863 dim : int
864 Number of dimensions.
865 num_dp : int
866 Number of datapoints.
867 factor : float
868 The bandwidth factor, obtained from `kde.covariance_factor`, with which
869 the covariance matrix is multiplied.
870 covariance : ndarray
871 The covariance matrix of *dataset*, scaled by the calculated bandwidth
872 (`kde.factor`).
873 inv_cov : ndarray
874 The inverse of *covariance*.
876 Methods
877 -------
878 kde.evaluate(points) : ndarray
879 Evaluate the estimated pdf on a provided set of points.
880 kde(points) : ndarray
881 Same as kde.evaluate(points)
882 """
884 # This implementation with minor modification was too good to pass up.
885 # from scipy: https://github.com/scipy/scipy/blob/master/scipy/stats/kde.py
887 def __init__(self, dataset, bw_method=None):
888 self.dataset = np.atleast_2d(dataset)
889 if not np.array(self.dataset).size > 1:
890 raise ValueError("`dataset` input should have multiple elements.")
892 self.dim, self.num_dp = np.array(self.dataset).shape
894 if bw_method is None:
895 pass
896 elif cbook._str_equal(bw_method, 'scott'):
897 self.covariance_factor = self.scotts_factor
898 elif cbook._str_equal(bw_method, 'silverman'):
899 self.covariance_factor = self.silverman_factor
900 elif isinstance(bw_method, Number):
901 self._bw_method = 'use constant'
902 self.covariance_factor = lambda: bw_method
903 elif callable(bw_method):
904 self._bw_method = bw_method
905 self.covariance_factor = lambda: self._bw_method(self)
906 else:
907 raise ValueError("`bw_method` should be 'scott', 'silverman', a "
908 "scalar or a callable")
910 # Computes the covariance matrix for each Gaussian kernel using
911 # covariance_factor().
913 self.factor = self.covariance_factor()
914 # Cache covariance and inverse covariance of the data
915 if not hasattr(self, '_data_inv_cov'):
916 self.data_covariance = np.atleast_2d(
917 np.cov(
918 self.dataset,
919 rowvar=1,
920 bias=False))
921 self.data_inv_cov = np.linalg.inv(self.data_covariance)
923 self.covariance = self.data_covariance * self.factor ** 2
924 self.inv_cov = self.data_inv_cov / self.factor ** 2
925 self.norm_factor = (np.sqrt(np.linalg.det(2 * np.pi * self.covariance))
926 * self.num_dp)
928 def scotts_factor(self):
929 return np.power(self.num_dp, -1. / (self.dim + 4))
931 def silverman_factor(self):
932 return np.power(
933 self.num_dp * (self.dim + 2.0) / 4.0, -1. / (self.dim + 4))
935 # Default method to calculate bandwidth, can be overwritten by subclass
936 covariance_factor = scotts_factor
938 def evaluate(self, points):
939 """
940 Evaluate the estimated pdf on a set of points.
942 Parameters
943 ----------
944 points : (# of dimensions, # of points)-array
945 Alternatively, a (# of dimensions,) vector can be passed in and
946 treated as a single point.
948 Returns
949 -------
950 (# of points,)-array
951 The values at each point.
953 Raises
954 ------
955 ValueError : if the dimensionality of the input points is different
956 than the dimensionality of the KDE.
958 """
959 points = np.atleast_2d(points)
961 dim, num_m = np.array(points).shape
962 if dim != self.dim:
963 raise ValueError("points have dimension {}, dataset has dimension "
964 "{}".format(dim, self.dim))
966 result = np.zeros(num_m)
968 if num_m >= self.num_dp:
969 # there are more points than data, so loop over data
970 for i in range(self.num_dp):
971 diff = self.dataset[:, i, np.newaxis] - points
972 tdiff = np.dot(self.inv_cov, diff)
973 energy = np.sum(diff * tdiff, axis=0) / 2.0
974 result = result + np.exp(-energy)
975 else:
976 # loop over points
977 for i in range(num_m):
978 diff = self.dataset - points[:, i, np.newaxis]
979 tdiff = np.dot(self.inv_cov, diff)
980 energy = np.sum(diff * tdiff, axis=0) / 2.0
981 result[i] = np.sum(np.exp(-energy), axis=0)
983 result = result / self.norm_factor
985 return result
987 __call__ = evaluate