1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
# Copyright (c) DataLab Platform Developers, BSD 3-Clause license, see LICENSE file.
"""
.. Peak Detection (see parent package :mod:`sigima.algorithms.signal`)
"""
from __future__ import annotations
import numpy as np
from sigima.tools.checks import check_1d_arrays
def peak_indices(
y, thres: float = 0.3, min_dist: int = 1, thres_abs: bool = False
) -> np.ndarray:
# Copyright (c) 2014 Lucas Hermann Negri
# Unmodified code snippet from PeakUtils 1.3.0
"""Peak detection routine.
Finds the numeric index of the peaks in *y* by taking its first order
difference. By using *thres* and *min_dist* parameters, it is possible
to reduce the number of detected peaks. *y* must be signed.
Parameters
----------
y : ndarray (signed)
1D amplitude data to search for peaks.
thres : float between [0., 1.]
Normalized threshold. Only the peaks with amplitude higher than the
threshold will be detected.
min_dist : int
Minimum distance between each detected peak. The peak with the highest
amplitude is preferred to satisfy this constraint.
thres_abs: boolean
If True, the thres value will be interpreted as an absolute value,
instead of a normalized threshold.
Returns
-------
ndarray
Array containing the numeric indices of the peaks that were detected
"""
if isinstance(y, np.ndarray) and np.issubdtype(y.dtype, np.unsignedinteger):
raise ValueError("y must be signed")
if not thres_abs:
thres = thres * (np.max(y) - np.min(y)) + np.min(y)
# compute first order difference
dy = np.diff(y)
# propagate left and right values successively to fill all plateau pixels
# (0-value)
(zeros,) = np.where(dy == 0)
# check if the signal is totally flat
if len(zeros) == len(y) - 1:
return np.array([])
if len(zeros):
# compute first order difference of zero indices
zeros_diff = np.diff(zeros)
# check when zeros are not chained together
(zeros_diff_not_one,) = np.add(np.where(zeros_diff != 1), 1)
# make an array of the chained zero indices
zero_plateaus = np.split(zeros, zeros_diff_not_one)
# fix if leftmost value in dy is zero
if zero_plateaus[0][0] == 0:
dy[zero_plateaus[0]] = dy[zero_plateaus[0][-1] + 1]
zero_plateaus.pop(0)
# fix if rightmost value of dy is zero
if len(zero_plateaus) > 0 and zero_plateaus[-1][-1] == len(dy) - 1:
dy[zero_plateaus[-1]] = dy[zero_plateaus[-1][0] - 1]
zero_plateaus.pop(-1)
# for each chain of zero indices
for plateau in zero_plateaus:
median = np.median(plateau)
# set leftmost values to leftmost non zero values
dy[plateau[plateau < median]] = dy[plateau[0] - 1]
# set rightmost and middle values to rightmost non zero values
dy[plateau[plateau >= median]] = dy[plateau[-1] + 1]
# find the peaks by using the first order difference
peaks = np.where(
(np.hstack([dy, 0.0]) < 0.0)
& (np.hstack([0.0, dy]) > 0.0)
& (np.greater(y, thres))
)[0]
# handle multiple peaks, respecting the minimum distance
if peaks.size > 1 and min_dist > 1:
highest = peaks[np.argsort(y[peaks])][::-1]
rem = np.ones(y.size, dtype=bool)
rem[peaks] = False
for peak in highest:
if not rem[peak]:
sl = slice(max(0, peak - min_dist), peak + min_dist + 1)
rem[sl] = True
rem[peak] = False
peaks = np.arange(y.size)[~rem]
return peaks
@check_1d_arrays
def xpeak(x: np.ndarray, y: np.ndarray) -> float:
"""Return default peak X-position (assuming a single peak).
Args:
x: X data
y: Y data
Returns:
Peak X-position
"""
peaks = peak_indices(y)
if peaks.size == 1:
return x[peaks[0]]
return np.average(x, weights=y)
|
