Coverage for larch/math/peaks.py: 6%
34 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1"""Peak detection algorithm borrowed from PeakUtils
3original author: Lucas Hermann Negri <lucashnegri@gmail.com>
5"""
6import numpy as np
8def peak_indices(y, threshold=0.1, min_dist=1):
9 """Peak detection routine.
11 Finds the indices of peaks in `y` data by taking its first order
12 derivative. and using the `thres` and `min_dist` parameters, to
13 control the number of peaks detected.
15 Parameters
16 ----------
17 y : 1D ndarray
18 1D amplitude data to search for peaks.
19 threshold : float [0.1]
20 Normalized peak threshold, as fraction of peak-to-peak of `y`.
21 Peaks with amplitude higher than the threshold will be detected.
22 min_dist: int [1]
23 Minimum index distance between detected peaks. When multiple
24 peaks are within this distance, the peak with highest amplitude
25 is preferred.
27 Returns
28 -------
29 ndarray
30 Array with numeric indexes of the peaks that were detected.
31 """
32 thres = threshold * np.ptp(y) + y.min()
33 min_dist = int(min_dist)
35 # compute first order difference
36 dy = np.diff(y)
38 # propagate left and right values successively to fill all
39 # plateau pixels (0-value)
40 zeros, = np.where(dy == 0)
42 # check if the signal is totally flat
43 if len(zeros) == len(y) - 1:
44 return np.array([])
46 if len(zeros):
47 # compute first order difference of zero indexes
48 zeros_diff = np.diff(zeros)
49 # check when zeros are not chained together
50 zeros_diff_not_one, = np.add(np.where(zeros_diff != 1), 1)
51 # make an array of the chained zero indexes
52 zero_plateaus = np.split(zeros, zeros_diff_not_one)
54 # fix if leftmost value in dy is zero
55 if zero_plateaus[0][0] == 0:
56 dy[zero_plateaus[0]] = dy[zero_plateaus[0][-1] + 1]
57 zero_plateaus.pop(0)
59 # fix if rightmost value of dy is zero
60 if len(zero_plateaus) and zero_plateaus[-1][-1] == len(dy) - 1:
61 dy[zero_plateaus[-1]] = dy[zero_plateaus[-1][0] - 1]
62 zero_plateaus.pop(-1)
64 # for each chain of zero indexes
65 for plateau in zero_plateaus:
66 median = np.median(plateau)
67 # set leftmost values to leftmost non zero values
68 dy[plateau[plateau < median]] = dy[plateau[0] - 1]
69 # set rightmost and middle values to rightmost non zero values
70 dy[plateau[plateau >= median]] = dy[plateau[-1] + 1]
72 # find the peaks by using the first order difference
73 peaks = np.where((np.hstack([dy, 0.0]) < 0.0) &
74 (np.hstack([0.0, dy]) > 0.0) &
75 (np.greater(y, thres)) )[0]
77 # handle multiple peaks, respecting the minimum distance
78 if peaks.size > 1 and min_dist > 1:
79 highest = peaks[np.argsort(y[peaks])][::-1]
80 rem = np.ones(y.size, dtype=bool)
81 rem[peaks] = False
83 for peak in highest:
84 if not rem[peak]:
85 sl = slice(max(0, peak - min_dist), peak + min_dist + 1)
86 rem[sl] = True
87 rem[peak] = False
89 peaks = np.arange(y.size)[~rem]
90 return peaks