Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/signal/_peak_finding.py : 10%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Functions for identifying peaks in signals.
3"""
4import math
5import numpy as np
7from scipy.signal.wavelets import cwt, ricker
8from scipy.stats import scoreatpercentile
10from ._peak_finding_utils import (
11 _local_maxima_1d,
12 _select_by_peak_distance,
13 _peak_prominences,
14 _peak_widths
15)
18__all__ = ['argrelmin', 'argrelmax', 'argrelextrema', 'peak_prominences',
19 'peak_widths', 'find_peaks', 'find_peaks_cwt']
22def _boolrelextrema(data, comparator, axis=0, order=1, mode='clip'):
23 """
24 Calculate the relative extrema of `data`.
26 Relative extrema are calculated by finding locations where
27 ``comparator(data[n], data[n+1:n+order+1])`` is True.
29 Parameters
30 ----------
31 data : ndarray
32 Array in which to find the relative extrema.
33 comparator : callable
34 Function to use to compare two data points.
35 Should take two arrays as arguments.
36 axis : int, optional
37 Axis over which to select from `data`. Default is 0.
38 order : int, optional
39 How many points on each side to use for the comparison
40 to consider ``comparator(n,n+x)`` to be True.
41 mode : str, optional
42 How the edges of the vector are treated. 'wrap' (wrap around) or
43 'clip' (treat overflow as the same as the last (or first) element).
44 Default 'clip'. See numpy.take.
46 Returns
47 -------
48 extrema : ndarray
49 Boolean array of the same shape as `data` that is True at an extrema,
50 False otherwise.
52 See also
53 --------
54 argrelmax, argrelmin
56 Examples
57 --------
58 >>> testdata = np.array([1,2,3,2,1])
59 >>> _boolrelextrema(testdata, np.greater, axis=0)
60 array([False, False, True, False, False], dtype=bool)
62 """
63 if((int(order) != order) or (order < 1)):
64 raise ValueError('Order must be an int >= 1')
66 datalen = data.shape[axis]
67 locs = np.arange(0, datalen)
69 results = np.ones(data.shape, dtype=bool)
70 main = data.take(locs, axis=axis, mode=mode)
71 for shift in range(1, order + 1):
72 plus = data.take(locs + shift, axis=axis, mode=mode)
73 minus = data.take(locs - shift, axis=axis, mode=mode)
74 results &= comparator(main, plus)
75 results &= comparator(main, minus)
76 if(~results.any()):
77 return results
78 return results
81def argrelmin(data, axis=0, order=1, mode='clip'):
82 """
83 Calculate the relative minima of `data`.
85 Parameters
86 ----------
87 data : ndarray
88 Array in which to find the relative minima.
89 axis : int, optional
90 Axis over which to select from `data`. Default is 0.
91 order : int, optional
92 How many points on each side to use for the comparison
93 to consider ``comparator(n, n+x)`` to be True.
94 mode : str, optional
95 How the edges of the vector are treated.
96 Available options are 'wrap' (wrap around) or 'clip' (treat overflow
97 as the same as the last (or first) element).
98 Default 'clip'. See numpy.take.
100 Returns
101 -------
102 extrema : tuple of ndarrays
103 Indices of the minima in arrays of integers. ``extrema[k]`` is
104 the array of indices of axis `k` of `data`. Note that the
105 return value is a tuple even when `data` is 1-D.
107 See Also
108 --------
109 argrelextrema, argrelmax, find_peaks
111 Notes
112 -----
113 This function uses `argrelextrema` with np.less as comparator. Therefore, it
114 requires a strict inequality on both sides of a value to consider it a
115 minimum. This means flat minima (more than one sample wide) are not detected.
116 In case of 1-D `data` `find_peaks` can be used to detect all
117 local minima, including flat ones, by calling it with negated `data`.
119 .. versionadded:: 0.11.0
121 Examples
122 --------
123 >>> from scipy.signal import argrelmin
124 >>> x = np.array([2, 1, 2, 3, 2, 0, 1, 0])
125 >>> argrelmin(x)
126 (array([1, 5]),)
127 >>> y = np.array([[1, 2, 1, 2],
128 ... [2, 2, 0, 0],
129 ... [5, 3, 4, 4]])
130 ...
131 >>> argrelmin(y, axis=1)
132 (array([0, 2]), array([2, 1]))
134 """
135 return argrelextrema(data, np.less, axis, order, mode)
138def argrelmax(data, axis=0, order=1, mode='clip'):
139 """
140 Calculate the relative maxima of `data`.
142 Parameters
143 ----------
144 data : ndarray
145 Array in which to find the relative maxima.
146 axis : int, optional
147 Axis over which to select from `data`. Default is 0.
148 order : int, optional
149 How many points on each side to use for the comparison
150 to consider ``comparator(n, n+x)`` to be True.
151 mode : str, optional
152 How the edges of the vector are treated.
153 Available options are 'wrap' (wrap around) or 'clip' (treat overflow
154 as the same as the last (or first) element).
155 Default 'clip'. See `numpy.take`.
157 Returns
158 -------
159 extrema : tuple of ndarrays
160 Indices of the maxima in arrays of integers. ``extrema[k]`` is
161 the array of indices of axis `k` of `data`. Note that the
162 return value is a tuple even when `data` is 1-D.
164 See Also
165 --------
166 argrelextrema, argrelmin, find_peaks
168 Notes
169 -----
170 This function uses `argrelextrema` with np.greater as comparator. Therefore,
171 it requires a strict inequality on both sides of a value to consider it a
172 maximum. This means flat maxima (more than one sample wide) are not detected.
173 In case of 1-D `data` `find_peaks` can be used to detect all
174 local maxima, including flat ones.
176 .. versionadded:: 0.11.0
178 Examples
179 --------
180 >>> from scipy.signal import argrelmax
181 >>> x = np.array([2, 1, 2, 3, 2, 0, 1, 0])
182 >>> argrelmax(x)
183 (array([3, 6]),)
184 >>> y = np.array([[1, 2, 1, 2],
185 ... [2, 2, 0, 0],
186 ... [5, 3, 4, 4]])
187 ...
188 >>> argrelmax(y, axis=1)
189 (array([0]), array([1]))
190 """
191 return argrelextrema(data, np.greater, axis, order, mode)
194def argrelextrema(data, comparator, axis=0, order=1, mode='clip'):
195 """
196 Calculate the relative extrema of `data`.
198 Parameters
199 ----------
200 data : ndarray
201 Array in which to find the relative extrema.
202 comparator : callable
203 Function to use to compare two data points.
204 Should take two arrays as arguments.
205 axis : int, optional
206 Axis over which to select from `data`. Default is 0.
207 order : int, optional
208 How many points on each side to use for the comparison
209 to consider ``comparator(n, n+x)`` to be True.
210 mode : str, optional
211 How the edges of the vector are treated. 'wrap' (wrap around) or
212 'clip' (treat overflow as the same as the last (or first) element).
213 Default is 'clip'. See `numpy.take`.
215 Returns
216 -------
217 extrema : tuple of ndarrays
218 Indices of the maxima in arrays of integers. ``extrema[k]`` is
219 the array of indices of axis `k` of `data`. Note that the
220 return value is a tuple even when `data` is 1-D.
222 See Also
223 --------
224 argrelmin, argrelmax
226 Notes
227 -----
229 .. versionadded:: 0.11.0
231 Examples
232 --------
233 >>> from scipy.signal import argrelextrema
234 >>> x = np.array([2, 1, 2, 3, 2, 0, 1, 0])
235 >>> argrelextrema(x, np.greater)
236 (array([3, 6]),)
237 >>> y = np.array([[1, 2, 1, 2],
238 ... [2, 2, 0, 0],
239 ... [5, 3, 4, 4]])
240 ...
241 >>> argrelextrema(y, np.less, axis=1)
242 (array([0, 2]), array([2, 1]))
244 """
245 results = _boolrelextrema(data, comparator,
246 axis, order, mode)
247 return np.nonzero(results)
250def _arg_x_as_expected(value):
251 """Ensure argument `x` is a 1-D C-contiguous array of dtype('float64').
253 Used in `find_peaks`, `peak_prominences` and `peak_widths` to make `x`
254 compatible with the signature of the wrapped Cython functions.
256 Returns
257 -------
258 value : ndarray
259 A 1-D C-contiguous array with dtype('float64').
260 """
261 value = np.asarray(value, order='C', dtype=np.float64)
262 if value.ndim != 1:
263 raise ValueError('`x` must be a 1-D array')
264 return value
267def _arg_peaks_as_expected(value):
268 """Ensure argument `peaks` is a 1-D C-contiguous array of dtype('intp').
270 Used in `peak_prominences` and `peak_widths` to make `peaks` compatible
271 with the signature of the wrapped Cython functions.
273 Returns
274 -------
275 value : ndarray
276 A 1-D C-contiguous array with dtype('intp').
277 """
278 value = np.asarray(value)
279 if value.size == 0:
280 # Empty arrays default to np.float64 but are valid input
281 value = np.array([], dtype=np.intp)
282 try:
283 # Safely convert to C-contiguous array of type np.intp
284 value = value.astype(np.intp, order='C', casting='safe',
285 subok=False, copy=False)
286 except TypeError:
287 raise TypeError("cannot safely cast `peaks` to dtype('intp')")
288 if value.ndim != 1:
289 raise ValueError('`peaks` must be a 1-D array')
290 return value
293def _arg_wlen_as_expected(value):
294 """Ensure argument `wlen` is of type `np.intp` and larger than 1.
296 Used in `peak_prominences` and `peak_widths`.
298 Returns
299 -------
300 value : np.intp
301 The original `value` rounded up to an integer or -1 if `value` was
302 None.
303 """
304 if value is None:
305 # _peak_prominences expects an intp; -1 signals that no value was
306 # supplied by the user
307 value = -1
308 elif 1 < value:
309 # Round up to a positive integer
310 if not np.can_cast(value, np.intp, "safe"):
311 value = math.ceil(value)
312 value = np.intp(value)
313 else:
314 raise ValueError('`wlen` must be larger than 1, was {}'
315 .format(value))
316 return value
319def peak_prominences(x, peaks, wlen=None):
320 """
321 Calculate the prominence of each peak in a signal.
323 The prominence of a peak measures how much a peak stands out from the
324 surrounding baseline of the signal and is defined as the vertical distance
325 between the peak and its lowest contour line.
327 Parameters
328 ----------
329 x : sequence
330 A signal with peaks.
331 peaks : sequence
332 Indices of peaks in `x`.
333 wlen : int, optional
334 A window length in samples that optionally limits the evaluated area for
335 each peak to a subset of `x`. The peak is always placed in the middle of
336 the window therefore the given length is rounded up to the next odd
337 integer. This parameter can speed up the calculation (see Notes).
339 Returns
340 -------
341 prominences : ndarray
342 The calculated prominences for each peak in `peaks`.
343 left_bases, right_bases : ndarray
344 The peaks' bases as indices in `x` to the left and right of each peak.
345 The higher base of each pair is a peak's lowest contour line.
347 Raises
348 ------
349 ValueError
350 If a value in `peaks` is an invalid index for `x`.
352 Warns
353 -----
354 PeakPropertyWarning
355 For indices in `peaks` that don't point to valid local maxima in `x`,
356 the returned prominence will be 0 and this warning is raised. This
357 also happens if `wlen` is smaller than the plateau size of a peak.
359 Warnings
360 --------
361 This function may return unexpected results for data containing NaNs. To
362 avoid this, NaNs should either be removed or replaced.
364 See Also
365 --------
366 find_peaks
367 Find peaks inside a signal based on peak properties.
368 peak_widths
369 Calculate the width of peaks.
371 Notes
372 -----
373 Strategy to compute a peak's prominence:
375 1. Extend a horizontal line from the current peak to the left and right
376 until the line either reaches the window border (see `wlen`) or
377 intersects the signal again at the slope of a higher peak. An
378 intersection with a peak of the same height is ignored.
379 2. On each side find the minimal signal value within the interval defined
380 above. These points are the peak's bases.
381 3. The higher one of the two bases marks the peak's lowest contour line. The
382 prominence can then be calculated as the vertical difference between the
383 peaks height itself and its lowest contour line.
385 Searching for the peak's bases can be slow for large `x` with periodic
386 behavior because large chunks or even the full signal need to be evaluated
387 for the first algorithmic step. This evaluation area can be limited with the
388 parameter `wlen` which restricts the algorithm to a window around the
389 current peak and can shorten the calculation time if the window length is
390 short in relation to `x`.
391 However, this may stop the algorithm from finding the true global contour
392 line if the peak's true bases are outside this window. Instead, a higher
393 contour line is found within the restricted window leading to a smaller
394 calculated prominence. In practice, this is only relevant for the highest set
395 of peaks in `x`. This behavior may even be used intentionally to calculate
396 "local" prominences.
398 .. versionadded:: 1.1.0
400 References
401 ----------
402 .. [1] Wikipedia Article for Topographic Prominence:
403 https://en.wikipedia.org/wiki/Topographic_prominence
405 Examples
406 --------
407 >>> from scipy.signal import find_peaks, peak_prominences
408 >>> import matplotlib.pyplot as plt
410 Create a test signal with two overlayed harmonics
412 >>> x = np.linspace(0, 6 * np.pi, 1000)
413 >>> x = np.sin(x) + 0.6 * np.sin(2.6 * x)
415 Find all peaks and calculate prominences
417 >>> peaks, _ = find_peaks(x)
418 >>> prominences = peak_prominences(x, peaks)[0]
419 >>> prominences
420 array([1.24159486, 0.47840168, 0.28470524, 3.10716793, 0.284603 ,
421 0.47822491, 2.48340261, 0.47822491])
423 Calculate the height of each peak's contour line and plot the results
425 >>> contour_heights = x[peaks] - prominences
426 >>> plt.plot(x)
427 >>> plt.plot(peaks, x[peaks], "x")
428 >>> plt.vlines(x=peaks, ymin=contour_heights, ymax=x[peaks])
429 >>> plt.show()
431 Let's evaluate a second example that demonstrates several edge cases for
432 one peak at index 5.
434 >>> x = np.array([0, 1, 0, 3, 1, 3, 0, 4, 0])
435 >>> peaks = np.array([5])
436 >>> plt.plot(x)
437 >>> plt.plot(peaks, x[peaks], "x")
438 >>> plt.show()
439 >>> peak_prominences(x, peaks) # -> (prominences, left_bases, right_bases)
440 (array([3.]), array([2]), array([6]))
442 Note how the peak at index 3 of the same height is not considered as a
443 border while searching for the left base. Instead, two minima at 0 and 2
444 are found in which case the one closer to the evaluated peak is always
445 chosen. On the right side, however, the base must be placed at 6 because the
446 higher peak represents the right border to the evaluated area.
448 >>> peak_prominences(x, peaks, wlen=3.1)
449 (array([2.]), array([4]), array([6]))
451 Here, we restricted the algorithm to a window from 3 to 7 (the length is 5
452 samples because `wlen` was rounded up to the next odd integer). Thus, the
453 only two candidates in the evaluated area are the two neighboring samples
454 and a smaller prominence is calculated.
455 """
456 x = _arg_x_as_expected(x)
457 peaks = _arg_peaks_as_expected(peaks)
458 wlen = _arg_wlen_as_expected(wlen)
459 return _peak_prominences(x, peaks, wlen)
462def peak_widths(x, peaks, rel_height=0.5, prominence_data=None, wlen=None):
463 """
464 Calculate the width of each peak in a signal.
466 This function calculates the width of a peak in samples at a relative
467 distance to the peak's height and prominence.
469 Parameters
470 ----------
471 x : sequence
472 A signal with peaks.
473 peaks : sequence
474 Indices of peaks in `x`.
475 rel_height : float, optional
476 Chooses the relative height at which the peak width is measured as a
477 percentage of its prominence. 1.0 calculates the width of the peak at
478 its lowest contour line while 0.5 evaluates at half the prominence
479 height. Must be at least 0. See notes for further explanation.
480 prominence_data : tuple, optional
481 A tuple of three arrays matching the output of `peak_prominences` when
482 called with the same arguments `x` and `peaks`. This data are calculated
483 internally if not provided.
484 wlen : int, optional
485 A window length in samples passed to `peak_prominences` as an optional
486 argument for internal calculation of `prominence_data`. This argument
487 is ignored if `prominence_data` is given.
489 Returns
490 -------
491 widths : ndarray
492 The widths for each peak in samples.
493 width_heights : ndarray
494 The height of the contour lines at which the `widths` where evaluated.
495 left_ips, right_ips : ndarray
496 Interpolated positions of left and right intersection points of a
497 horizontal line at the respective evaluation height.
499 Raises
500 ------
501 ValueError
502 If `prominence_data` is supplied but doesn't satisfy the condition
503 ``0 <= left_base <= peak <= right_base < x.shape[0]`` for each peak,
504 has the wrong dtype, is not C-contiguous or does not have the same
505 shape.
507 Warns
508 -----
509 PeakPropertyWarning
510 Raised if any calculated width is 0. This may stem from the supplied
511 `prominence_data` or if `rel_height` is set to 0.
513 Warnings
514 --------
515 This function may return unexpected results for data containing NaNs. To
516 avoid this, NaNs should either be removed or replaced.
518 See Also
519 --------
520 find_peaks
521 Find peaks inside a signal based on peak properties.
522 peak_prominences
523 Calculate the prominence of peaks.
525 Notes
526 -----
527 The basic algorithm to calculate a peak's width is as follows:
529 * Calculate the evaluation height :math:`h_{eval}` with the formula
530 :math:`h_{eval} = h_{Peak} - P \\cdot R`, where :math:`h_{Peak}` is the
531 height of the peak itself, :math:`P` is the peak's prominence and
532 :math:`R` a positive ratio specified with the argument `rel_height`.
533 * Draw a horizontal line at the evaluation height to both sides, starting at
534 the peak's current vertical position until the lines either intersect a
535 slope, the signal border or cross the vertical position of the peak's
536 base (see `peak_prominences` for an definition). For the first case,
537 intersection with the signal, the true intersection point is estimated
538 with linear interpolation.
539 * Calculate the width as the horizontal distance between the chosen
540 endpoints on both sides. As a consequence of this the maximal possible
541 width for each peak is the horizontal distance between its bases.
543 As shown above to calculate a peak's width its prominence and bases must be
544 known. You can supply these yourself with the argument `prominence_data`.
545 Otherwise, they are internally calculated (see `peak_prominences`).
547 .. versionadded:: 1.1.0
549 Examples
550 --------
551 >>> from scipy.signal import chirp, find_peaks, peak_widths
552 >>> import matplotlib.pyplot as plt
554 Create a test signal with two overlayed harmonics
556 >>> x = np.linspace(0, 6 * np.pi, 1000)
557 >>> x = np.sin(x) + 0.6 * np.sin(2.6 * x)
559 Find all peaks and calculate their widths at the relative height of 0.5
560 (contour line at half the prominence height) and 1 (at the lowest contour
561 line at full prominence height).
563 >>> peaks, _ = find_peaks(x)
564 >>> results_half = peak_widths(x, peaks, rel_height=0.5)
565 >>> results_half[0] # widths
566 array([ 64.25172825, 41.29465463, 35.46943289, 104.71586081,
567 35.46729324, 41.30429622, 181.93835853, 45.37078546])
568 >>> results_full = peak_widths(x, peaks, rel_height=1)
569 >>> results_full[0] # widths
570 array([181.9396084 , 72.99284945, 61.28657872, 373.84622694,
571 61.78404617, 72.48822812, 253.09161876, 79.36860878])
573 Plot signal, peaks and contour lines at which the widths where calculated
575 >>> plt.plot(x)
576 >>> plt.plot(peaks, x[peaks], "x")
577 >>> plt.hlines(*results_half[1:], color="C2")
578 >>> plt.hlines(*results_full[1:], color="C3")
579 >>> plt.show()
580 """
581 x = _arg_x_as_expected(x)
582 peaks = _arg_peaks_as_expected(peaks)
583 if prominence_data is None:
584 # Calculate prominence if not supplied and use wlen if supplied.
585 wlen = _arg_wlen_as_expected(wlen)
586 prominence_data = _peak_prominences(x, peaks, wlen)
587 return _peak_widths(x, peaks, rel_height, *prominence_data)
590def _unpack_condition_args(interval, x, peaks):
591 """
592 Parse condition arguments for `find_peaks`.
594 Parameters
595 ----------
596 interval : number or ndarray or sequence
597 Either a number or ndarray or a 2-element sequence of the former. The
598 first value is always interpreted as `imin` and the second, if supplied,
599 as `imax`.
600 x : ndarray
601 The signal with `peaks`.
602 peaks : ndarray
603 An array with indices used to reduce `imin` and / or `imax` if those are
604 arrays.
606 Returns
607 -------
608 imin, imax : number or ndarray or None
609 Minimal and maximal value in `argument`.
611 Raises
612 ------
613 ValueError :
614 If interval border is given as array and its size does not match the size
615 of `x`.
617 Notes
618 -----
620 .. versionadded:: 1.1.0
621 """
622 try:
623 imin, imax = interval
624 except (TypeError, ValueError):
625 imin, imax = (interval, None)
627 # Reduce arrays if arrays
628 if isinstance(imin, np.ndarray):
629 if imin.size != x.size:
630 raise ValueError('array size of lower interval border must match x')
631 imin = imin[peaks]
632 if isinstance(imax, np.ndarray):
633 if imax.size != x.size:
634 raise ValueError('array size of upper interval border must match x')
635 imax = imax[peaks]
637 return imin, imax
640def _select_by_property(peak_properties, pmin, pmax):
641 """
642 Evaluate where the generic property of peaks confirms to an interval.
644 Parameters
645 ----------
646 peak_properties : ndarray
647 An array with properties for each peak.
648 pmin : None or number or ndarray
649 Lower interval boundary for `peak_properties`. ``None`` is interpreted as
650 an open border.
651 pmax : None or number or ndarray
652 Upper interval boundary for `peak_properties`. ``None`` is interpreted as
653 an open border.
655 Returns
656 -------
657 keep : bool
658 A boolean mask evaluating to true where `peak_properties` confirms to the
659 interval.
661 See Also
662 --------
663 find_peaks
665 Notes
666 -----
668 .. versionadded:: 1.1.0
669 """
670 keep = np.ones(peak_properties.size, dtype=bool)
671 if pmin is not None:
672 keep &= (pmin <= peak_properties)
673 if pmax is not None:
674 keep &= (peak_properties <= pmax)
675 return keep
678def _select_by_peak_threshold(x, peaks, tmin, tmax):
679 """
680 Evaluate which peaks fulfill the threshold condition.
682 Parameters
683 ----------
684 x : ndarray
685 A 1-D array which is indexable by `peaks`.
686 peaks : ndarray
687 Indices of peaks in `x`.
688 tmin, tmax : scalar or ndarray or None
689 Minimal and / or maximal required thresholds. If supplied as ndarrays
690 their size must match `peaks`. ``None`` is interpreted as an open
691 border.
693 Returns
694 -------
695 keep : bool
696 A boolean mask evaluating to true where `peaks` fulfill the threshold
697 condition.
698 left_thresholds, right_thresholds : ndarray
699 Array matching `peak` containing the thresholds of each peak on
700 both sides.
702 Notes
703 -----
705 .. versionadded:: 1.1.0
706 """
707 # Stack thresholds on both sides to make min / max operations easier:
708 # tmin is compared with the smaller, and tmax with the greater thresold to
709 # each peak's side
710 stacked_thresholds = np.vstack([x[peaks] - x[peaks - 1],
711 x[peaks] - x[peaks + 1]])
712 keep = np.ones(peaks.size, dtype=bool)
713 if tmin is not None:
714 min_thresholds = np.min(stacked_thresholds, axis=0)
715 keep &= (tmin <= min_thresholds)
716 if tmax is not None:
717 max_thresholds = np.max(stacked_thresholds, axis=0)
718 keep &= (max_thresholds <= tmax)
720 return keep, stacked_thresholds[0], stacked_thresholds[1]
723def find_peaks(x, height=None, threshold=None, distance=None,
724 prominence=None, width=None, wlen=None, rel_height=0.5,
725 plateau_size=None):
726 """
727 Find peaks inside a signal based on peak properties.
729 This function takes a 1-D array and finds all local maxima by
730 simple comparison of neighboring values. Optionally, a subset of these
731 peaks can be selected by specifying conditions for a peak's properties.
733 Parameters
734 ----------
735 x : sequence
736 A signal with peaks.
737 height : number or ndarray or sequence, optional
738 Required height of peaks. Either a number, ``None``, an array matching
739 `x` or a 2-element sequence of the former. The first element is
740 always interpreted as the minimal and the second, if supplied, as the
741 maximal required height.
742 threshold : number or ndarray or sequence, optional
743 Required threshold of peaks, the vertical distance to its neighboring
744 samples. Either a number, ``None``, an array matching `x` or a
745 2-element sequence of the former. The first element is always
746 interpreted as the minimal and the second, if supplied, as the maximal
747 required threshold.
748 distance : number, optional
749 Required minimal horizontal distance (>= 1) in samples between
750 neighbouring peaks. Smaller peaks are removed first until the condition
751 is fulfilled for all remaining peaks.
752 prominence : number or ndarray or sequence, optional
753 Required prominence of peaks. Either a number, ``None``, an array
754 matching `x` or a 2-element sequence of the former. The first
755 element is always interpreted as the minimal and the second, if
756 supplied, as the maximal required prominence.
757 width : number or ndarray or sequence, optional
758 Required width of peaks in samples. Either a number, ``None``, an array
759 matching `x` or a 2-element sequence of the former. The first
760 element is always interpreted as the minimal and the second, if
761 supplied, as the maximal required width.
762 wlen : int, optional
763 Used for calculation of the peaks prominences, thus it is only used if
764 one of the arguments `prominence` or `width` is given. See argument
765 `wlen` in `peak_prominences` for a full description of its effects.
766 rel_height : float, optional
767 Used for calculation of the peaks width, thus it is only used if `width`
768 is given. See argument `rel_height` in `peak_widths` for a full
769 description of its effects.
770 plateau_size : number or ndarray or sequence, optional
771 Required size of the flat top of peaks in samples. Either a number,
772 ``None``, an array matching `x` or a 2-element sequence of the former.
773 The first element is always interpreted as the minimal and the second,
774 if supplied as the maximal required plateau size.
776 .. versionadded:: 1.2.0
778 Returns
779 -------
780 peaks : ndarray
781 Indices of peaks in `x` that satisfy all given conditions.
782 properties : dict
783 A dictionary containing properties of the returned peaks which were
784 calculated as intermediate results during evaluation of the specified
785 conditions:
787 * 'peak_heights'
788 If `height` is given, the height of each peak in `x`.
789 * 'left_thresholds', 'right_thresholds'
790 If `threshold` is given, these keys contain a peaks vertical
791 distance to its neighbouring samples.
792 * 'prominences', 'right_bases', 'left_bases'
793 If `prominence` is given, these keys are accessible. See
794 `peak_prominences` for a description of their content.
795 * 'width_heights', 'left_ips', 'right_ips'
796 If `width` is given, these keys are accessible. See `peak_widths`
797 for a description of their content.
798 * 'plateau_sizes', left_edges', 'right_edges'
799 If `plateau_size` is given, these keys are accessible and contain
800 the indices of a peak's edges (edges are still part of the
801 plateau) and the calculated plateau sizes.
803 .. versionadded:: 1.2.0
805 To calculate and return properties without excluding peaks, provide the
806 open interval ``(None, None)`` as a value to the appropriate argument
807 (excluding `distance`).
809 Warns
810 -----
811 PeakPropertyWarning
812 Raised if a peak's properties have unexpected values (see
813 `peak_prominences` and `peak_widths`).
815 Warnings
816 --------
817 This function may return unexpected results for data containing NaNs. To
818 avoid this, NaNs should either be removed or replaced.
820 See Also
821 --------
822 find_peaks_cwt
823 Find peaks using the wavelet transformation.
824 peak_prominences
825 Directly calculate the prominence of peaks.
826 peak_widths
827 Directly calculate the width of peaks.
829 Notes
830 -----
831 In the context of this function, a peak or local maximum is defined as any
832 sample whose two direct neighbours have a smaller amplitude. For flat peaks
833 (more than one sample of equal amplitude wide) the index of the middle
834 sample is returned (rounded down in case the number of samples is even).
835 For noisy signals the peak locations can be off because the noise might
836 change the position of local maxima. In those cases consider smoothing the
837 signal before searching for peaks or use other peak finding and fitting
838 methods (like `find_peaks_cwt`).
840 Some additional comments on specifying conditions:
842 * Almost all conditions (excluding `distance`) can be given as half-open or
843 closed intervals, e.g., ``1`` or ``(1, None)`` defines the half-open
844 interval :math:`[1, \\infty]` while ``(None, 1)`` defines the interval
845 :math:`[-\\infty, 1]`. The open interval ``(None, None)`` can be specified
846 as well, which returns the matching properties without exclusion of peaks.
847 * The border is always included in the interval used to select valid peaks.
848 * For several conditions the interval borders can be specified with
849 arrays matching `x` in shape which enables dynamic constrains based on
850 the sample position.
851 * The conditions are evaluated in the following order: `plateau_size`,
852 `height`, `threshold`, `distance`, `prominence`, `width`. In most cases
853 this order is the fastest one because faster operations are applied first
854 to reduce the number of peaks that need to be evaluated later.
855 * While indices in `peaks` are guaranteed to be at least `distance` samples
856 apart, edges of flat peaks may be closer than the allowed `distance`.
857 * Use `wlen` to reduce the time it takes to evaluate the conditions for
858 `prominence` or `width` if `x` is large or has many local maxima
859 (see `peak_prominences`).
861 .. versionadded:: 1.1.0
863 Examples
864 --------
865 To demonstrate this function's usage we use a signal `x` supplied with
866 SciPy (see `scipy.misc.electrocardiogram`). Let's find all peaks (local
867 maxima) in `x` whose amplitude lies above 0.
869 >>> import matplotlib.pyplot as plt
870 >>> from scipy.misc import electrocardiogram
871 >>> from scipy.signal import find_peaks
872 >>> x = electrocardiogram()[2000:4000]
873 >>> peaks, _ = find_peaks(x, height=0)
874 >>> plt.plot(x)
875 >>> plt.plot(peaks, x[peaks], "x")
876 >>> plt.plot(np.zeros_like(x), "--", color="gray")
877 >>> plt.show()
879 We can select peaks below 0 with ``height=(None, 0)`` or use arrays matching
880 `x` in size to reflect a changing condition for different parts of the
881 signal.
883 >>> border = np.sin(np.linspace(0, 3 * np.pi, x.size))
884 >>> peaks, _ = find_peaks(x, height=(-border, border))
885 >>> plt.plot(x)
886 >>> plt.plot(-border, "--", color="gray")
887 >>> plt.plot(border, ":", color="gray")
888 >>> plt.plot(peaks, x[peaks], "x")
889 >>> plt.show()
891 Another useful condition for periodic signals can be given with the
892 `distance` argument. In this case, we can easily select the positions of
893 QRS complexes within the electrocardiogram (ECG) by demanding a distance of
894 at least 150 samples.
896 >>> peaks, _ = find_peaks(x, distance=150)
897 >>> np.diff(peaks)
898 array([186, 180, 177, 171, 177, 169, 167, 164, 158, 162, 172])
899 >>> plt.plot(x)
900 >>> plt.plot(peaks, x[peaks], "x")
901 >>> plt.show()
903 Especially for noisy signals peaks can be easily grouped by their
904 prominence (see `peak_prominences`). E.g., we can select all peaks except
905 for the mentioned QRS complexes by limiting the allowed prominence to 0.6.
907 >>> peaks, properties = find_peaks(x, prominence=(None, 0.6))
908 >>> properties["prominences"].max()
909 0.5049999999999999
910 >>> plt.plot(x)
911 >>> plt.plot(peaks, x[peaks], "x")
912 >>> plt.show()
914 And, finally, let's examine a different section of the ECG which contains
915 beat forms of different shape. To select only the atypical heart beats, we
916 combine two conditions: a minimal prominence of 1 and width of at least 20
917 samples.
919 >>> x = electrocardiogram()[17000:18000]
920 >>> peaks, properties = find_peaks(x, prominence=1, width=20)
921 >>> properties["prominences"], properties["widths"]
922 (array([1.495, 2.3 ]), array([36.93773946, 39.32723577]))
923 >>> plt.plot(x)
924 >>> plt.plot(peaks, x[peaks], "x")
925 >>> plt.vlines(x=peaks, ymin=x[peaks] - properties["prominences"],
926 ... ymax = x[peaks], color = "C1")
927 >>> plt.hlines(y=properties["width_heights"], xmin=properties["left_ips"],
928 ... xmax=properties["right_ips"], color = "C1")
929 >>> plt.show()
930 """
931 # _argmaxima1d expects array of dtype 'float64'
932 x = _arg_x_as_expected(x)
933 if distance is not None and distance < 1:
934 raise ValueError('`distance` must be greater or equal to 1')
936 peaks, left_edges, right_edges = _local_maxima_1d(x)
937 properties = {}
939 if plateau_size is not None:
940 # Evaluate plateau size
941 plateau_sizes = right_edges - left_edges + 1
942 pmin, pmax = _unpack_condition_args(plateau_size, x, peaks)
943 keep = _select_by_property(plateau_sizes, pmin, pmax)
944 peaks = peaks[keep]
945 properties["plateau_sizes"] = plateau_sizes
946 properties["left_edges"] = left_edges
947 properties["right_edges"] = right_edges
948 properties = {key: array[keep] for key, array in properties.items()}
950 if height is not None:
951 # Evaluate height condition
952 peak_heights = x[peaks]
953 hmin, hmax = _unpack_condition_args(height, x, peaks)
954 keep = _select_by_property(peak_heights, hmin, hmax)
955 peaks = peaks[keep]
956 properties["peak_heights"] = peak_heights
957 properties = {key: array[keep] for key, array in properties.items()}
959 if threshold is not None:
960 # Evaluate threshold condition
961 tmin, tmax = _unpack_condition_args(threshold, x, peaks)
962 keep, left_thresholds, right_thresholds = _select_by_peak_threshold(
963 x, peaks, tmin, tmax)
964 peaks = peaks[keep]
965 properties["left_thresholds"] = left_thresholds
966 properties["right_thresholds"] = right_thresholds
967 properties = {key: array[keep] for key, array in properties.items()}
969 if distance is not None:
970 # Evaluate distance condition
971 keep = _select_by_peak_distance(peaks, x[peaks], distance)
972 peaks = peaks[keep]
973 properties = {key: array[keep] for key, array in properties.items()}
975 if prominence is not None or width is not None:
976 # Calculate prominence (required for both conditions)
977 wlen = _arg_wlen_as_expected(wlen)
978 properties.update(zip(
979 ['prominences', 'left_bases', 'right_bases'],
980 _peak_prominences(x, peaks, wlen=wlen)
981 ))
983 if prominence is not None:
984 # Evaluate prominence condition
985 pmin, pmax = _unpack_condition_args(prominence, x, peaks)
986 keep = _select_by_property(properties['prominences'], pmin, pmax)
987 peaks = peaks[keep]
988 properties = {key: array[keep] for key, array in properties.items()}
990 if width is not None:
991 # Calculate widths
992 properties.update(zip(
993 ['widths', 'width_heights', 'left_ips', 'right_ips'],
994 _peak_widths(x, peaks, rel_height, properties['prominences'],
995 properties['left_bases'], properties['right_bases'])
996 ))
997 # Evaluate width condition
998 wmin, wmax = _unpack_condition_args(width, x, peaks)
999 keep = _select_by_property(properties['widths'], wmin, wmax)
1000 peaks = peaks[keep]
1001 properties = {key: array[keep] for key, array in properties.items()}
1003 return peaks, properties
1006def _identify_ridge_lines(matr, max_distances, gap_thresh):
1007 """
1008 Identify ridges in the 2-D matrix.
1010 Expect that the width of the wavelet feature increases with increasing row
1011 number.
1013 Parameters
1014 ----------
1015 matr : 2-D ndarray
1016 Matrix in which to identify ridge lines.
1017 max_distances : 1-D sequence
1018 At each row, a ridge line is only connected
1019 if the relative max at row[n] is within
1020 `max_distances`[n] from the relative max at row[n+1].
1021 gap_thresh : int
1022 If a relative maximum is not found within `max_distances`,
1023 there will be a gap. A ridge line is discontinued if
1024 there are more than `gap_thresh` points without connecting
1025 a new relative maximum.
1027 Returns
1028 -------
1029 ridge_lines : tuple
1030 Tuple of 2 1-D sequences. `ridge_lines`[ii][0] are the rows of the
1031 ii-th ridge-line, `ridge_lines`[ii][1] are the columns. Empty if none
1032 found. Each ridge-line will be sorted by row (increasing), but the
1033 order of the ridge lines is not specified.
1035 References
1036 ----------
1037 .. [1] Bioinformatics (2006) 22 (17): 2059-2065.
1038 :doi:`10.1093/bioinformatics/btl355`
1040 Examples
1041 --------
1042 >>> data = np.random.rand(5,5)
1043 >>> ridge_lines = _identify_ridge_lines(data, 1, 1)
1045 Notes
1046 -----
1047 This function is intended to be used in conjunction with `cwt`
1048 as part of `find_peaks_cwt`.
1050 """
1051 if(len(max_distances) < matr.shape[0]):
1052 raise ValueError('Max_distances must have at least as many rows '
1053 'as matr')
1055 all_max_cols = _boolrelextrema(matr, np.greater, axis=1, order=1)
1056 # Highest row for which there are any relative maxima
1057 has_relmax = np.nonzero(all_max_cols.any(axis=1))[0]
1058 if(len(has_relmax) == 0):
1059 return []
1060 start_row = has_relmax[-1]
1061 # Each ridge line is a 3-tuple:
1062 # rows, cols,Gap number
1063 ridge_lines = [[[start_row],
1064 [col],
1065 0] for col in np.nonzero(all_max_cols[start_row])[0]]
1066 final_lines = []
1067 rows = np.arange(start_row - 1, -1, -1)
1068 cols = np.arange(0, matr.shape[1])
1069 for row in rows:
1070 this_max_cols = cols[all_max_cols[row]]
1072 # Increment gap number of each line,
1073 # set it to zero later if appropriate
1074 for line in ridge_lines:
1075 line[2] += 1
1077 # XXX These should always be all_max_cols[row]
1078 # But the order might be different. Might be an efficiency gain
1079 # to make sure the order is the same and avoid this iteration
1080 prev_ridge_cols = np.array([line[1][-1] for line in ridge_lines])
1081 # Look through every relative maximum found at current row
1082 # Attempt to connect them with existing ridge lines.
1083 for ind, col in enumerate(this_max_cols):
1084 # If there is a previous ridge line within
1085 # the max_distance to connect to, do so.
1086 # Otherwise start a new one.
1087 line = None
1088 if(len(prev_ridge_cols) > 0):
1089 diffs = np.abs(col - prev_ridge_cols)
1090 closest = np.argmin(diffs)
1091 if diffs[closest] <= max_distances[row]:
1092 line = ridge_lines[closest]
1093 if(line is not None):
1094 # Found a point close enough, extend current ridge line
1095 line[1].append(col)
1096 line[0].append(row)
1097 line[2] = 0
1098 else:
1099 new_line = [[row],
1100 [col],
1101 0]
1102 ridge_lines.append(new_line)
1104 # Remove the ridge lines with gap_number too high
1105 # XXX Modifying a list while iterating over it.
1106 # Should be safe, since we iterate backwards, but
1107 # still tacky.
1108 for ind in range(len(ridge_lines) - 1, -1, -1):
1109 line = ridge_lines[ind]
1110 if line[2] > gap_thresh:
1111 final_lines.append(line)
1112 del ridge_lines[ind]
1114 out_lines = []
1115 for line in (final_lines + ridge_lines):
1116 sortargs = np.array(np.argsort(line[0]))
1117 rows, cols = np.zeros_like(sortargs), np.zeros_like(sortargs)
1118 rows[sortargs] = line[0]
1119 cols[sortargs] = line[1]
1120 out_lines.append([rows, cols])
1122 return out_lines
1125def _filter_ridge_lines(cwt, ridge_lines, window_size=None, min_length=None,
1126 min_snr=1, noise_perc=10):
1127 """
1128 Filter ridge lines according to prescribed criteria. Intended
1129 to be used for finding relative maxima.
1131 Parameters
1132 ----------
1133 cwt : 2-D ndarray
1134 Continuous wavelet transform from which the `ridge_lines` were defined.
1135 ridge_lines : 1-D sequence
1136 Each element should contain 2 sequences, the rows and columns
1137 of the ridge line (respectively).
1138 window_size : int, optional
1139 Size of window to use to calculate noise floor.
1140 Default is ``cwt.shape[1] / 20``.
1141 min_length : int, optional
1142 Minimum length a ridge line needs to be acceptable.
1143 Default is ``cwt.shape[0] / 4``, ie 1/4-th the number of widths.
1144 min_snr : float, optional
1145 Minimum SNR ratio. Default 1. The signal is the value of
1146 the cwt matrix at the shortest length scale (``cwt[0, loc]``), the
1147 noise is the `noise_perc`th percentile of datapoints contained within a
1148 window of `window_size` around ``cwt[0, loc]``.
1149 noise_perc : float, optional
1150 When calculating the noise floor, percentile of data points
1151 examined below which to consider noise. Calculated using
1152 scipy.stats.scoreatpercentile.
1154 References
1155 ----------
1156 .. [1] Bioinformatics (2006) 22 (17): 2059-2065.
1157 :doi:`10.1093/bioinformatics/btl355`
1159 """
1160 num_points = cwt.shape[1]
1161 if min_length is None:
1162 min_length = np.ceil(cwt.shape[0] / 4)
1163 if window_size is None:
1164 window_size = np.ceil(num_points / 20)
1166 window_size = int(window_size)
1167 hf_window, odd = divmod(window_size, 2)
1169 # Filter based on SNR
1170 row_one = cwt[0, :]
1171 noises = np.zeros_like(row_one)
1172 for ind, val in enumerate(row_one):
1173 window_start = max(ind - hf_window, 0)
1174 window_end = min(ind + hf_window + odd, num_points)
1175 noises[ind] = scoreatpercentile(row_one[window_start:window_end],
1176 per=noise_perc)
1178 def filt_func(line):
1179 if len(line[0]) < min_length:
1180 return False
1181 snr = abs(cwt[line[0][0], line[1][0]] / noises[line[1][0]])
1182 if snr < min_snr:
1183 return False
1184 return True
1186 return list(filter(filt_func, ridge_lines))
1189def find_peaks_cwt(vector, widths, wavelet=None, max_distances=None,
1190 gap_thresh=None, min_length=None,
1191 min_snr=1, noise_perc=10, window_size=None):
1192 """
1193 Find peaks in a 1-D array with wavelet transformation.
1195 The general approach is to smooth `vector` by convolving it with
1196 `wavelet(width)` for each width in `widths`. Relative maxima which
1197 appear at enough length scales, and with sufficiently high SNR, are
1198 accepted.
1200 Parameters
1201 ----------
1202 vector : ndarray
1203 1-D array in which to find the peaks.
1204 widths : sequence
1205 1-D array of widths to use for calculating the CWT matrix. In general,
1206 this range should cover the expected width of peaks of interest.
1207 wavelet : callable, optional
1208 Should take two parameters and return a 1-D array to convolve
1209 with `vector`. The first parameter determines the number of points
1210 of the returned wavelet array, the second parameter is the scale
1211 (`width`) of the wavelet. Should be normalized and symmetric.
1212 Default is the ricker wavelet.
1213 max_distances : ndarray, optional
1214 At each row, a ridge line is only connected if the relative max at
1215 row[n] is within ``max_distances[n]`` from the relative max at
1216 ``row[n+1]``. Default value is ``widths/4``.
1217 gap_thresh : float, optional
1218 If a relative maximum is not found within `max_distances`,
1219 there will be a gap. A ridge line is discontinued if there are more
1220 than `gap_thresh` points without connecting a new relative maximum.
1221 Default is the first value of the widths array i.e. widths[0].
1222 min_length : int, optional
1223 Minimum length a ridge line needs to be acceptable.
1224 Default is ``cwt.shape[0] / 4``, ie 1/4-th the number of widths.
1225 min_snr : float, optional
1226 Minimum SNR ratio. Default 1. The signal is the value of
1227 the cwt matrix at the shortest length scale (``cwt[0, loc]``), the
1228 noise is the `noise_perc`th percentile of datapoints contained within a
1229 window of `window_size` around ``cwt[0, loc]``.
1230 noise_perc : float, optional
1231 When calculating the noise floor, percentile of data points
1232 examined below which to consider noise. Calculated using
1233 `stats.scoreatpercentile`. Default is 10.
1234 window_size : int, optional
1235 Size of window to use to calculate noise floor.
1236 Default is ``cwt.shape[1] / 20``.
1238 Returns
1239 -------
1240 peaks_indices : ndarray
1241 Indices of the locations in the `vector` where peaks were found.
1242 The list is sorted.
1244 See Also
1245 --------
1246 cwt
1247 Continuous wavelet transform.
1248 find_peaks
1249 Find peaks inside a signal based on peak properties.
1251 Notes
1252 -----
1253 This approach was designed for finding sharp peaks among noisy data,
1254 however with proper parameter selection it should function well for
1255 different peak shapes.
1257 The algorithm is as follows:
1258 1. Perform a continuous wavelet transform on `vector`, for the supplied
1259 `widths`. This is a convolution of `vector` with `wavelet(width)` for
1260 each width in `widths`. See `cwt`.
1261 2. Identify "ridge lines" in the cwt matrix. These are relative maxima
1262 at each row, connected across adjacent rows. See identify_ridge_lines
1263 3. Filter the ridge_lines using filter_ridge_lines.
1265 .. versionadded:: 0.11.0
1267 References
1268 ----------
1269 .. [1] Bioinformatics (2006) 22 (17): 2059-2065.
1270 :doi:`10.1093/bioinformatics/btl355`
1272 Examples
1273 --------
1274 >>> from scipy import signal
1275 >>> xs = np.arange(0, np.pi, 0.05)
1276 >>> data = np.sin(xs)
1277 >>> peakind = signal.find_peaks_cwt(data, np.arange(1,10))
1278 >>> peakind, xs[peakind], data[peakind]
1279 ([32], array([ 1.6]), array([ 0.9995736]))
1281 """
1282 widths = np.asarray(widths)
1284 if gap_thresh is None:
1285 gap_thresh = np.ceil(widths[0])
1286 if max_distances is None:
1287 max_distances = widths / 4.0
1288 if wavelet is None:
1289 wavelet = ricker
1291 cwt_dat = cwt(vector, wavelet, widths, window_size=window_size)
1292 ridge_lines = _identify_ridge_lines(cwt_dat, max_distances, gap_thresh)
1293 filtered = _filter_ridge_lines(cwt_dat, ridge_lines, min_length=min_length,
1294 window_size=window_size, min_snr=min_snr,
1295 noise_perc=noise_perc)
1296 max_locs = np.asarray([x[1][0] for x in filtered])
1297 max_locs.sort()
1299 return max_locs