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

1"""Peak detection algorithm borrowed from PeakUtils 

2 

3original author: Lucas Hermann Negri <lucashnegri@gmail.com> 

4 

5""" 

6import numpy as np 

7 

8def peak_indices(y, threshold=0.1, min_dist=1): 

9 """Peak detection routine. 

10 

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. 

14 

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. 

26 

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) 

34 

35 # compute first order difference 

36 dy = np.diff(y) 

37 

38 # propagate left and right values successively to fill all 

39 # plateau pixels (0-value) 

40 zeros, = np.where(dy == 0) 

41 

42 # check if the signal is totally flat 

43 if len(zeros) == len(y) - 1: 

44 return np.array([]) 

45 

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) 

53 

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) 

58 

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) 

63 

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] 

71 

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] 

76 

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 

82 

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 

88 

89 peaks = np.arange(y.size)[~rem] 

90 return peaks