-
Notifications
You must be signed in to change notification settings - Fork 0
/
detect_peaks.R
295 lines (269 loc) · 11.6 KB
/
detect_peaks.R
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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
# Code retrieved from package TagTools of DeRuiter et al.
# [https://github.com/stacyderuiter/TagTools/blob/master/R/tagtools/R/detect_peaks.R]
#' Detect peaks in signal vector data
#'
#' This function detects peaks in time series data that exceed a specified threshold and returns each peak's start time, end time, maximum peak value, time of the maximum peak value, threshold level, and blanking time.
#' @param data A vector (of all positive values) or matrix of data to be used in peak detection. If data is a matrix, you must specify a FUN to be applied to data.
#' @param FUN A function to be applied to data before the data is run through the peak detector. Only specify the function name (i.e. njerk). If left blank, the data input will be immediately passed through the peak detector.
#' @param sr The sampling rate in Hz of the date. This is the same as fs in other tagtools functions. This is used to calculate the bktime in the case that the input for bktime is missing.
#' @param thresh The threshold level above which peaks in signal are detected. Inputs must be in the same units as the signal. If the input for thresh is missing/empty, the default level is the 0.99 quantile
#' @param bktime The specified length of time (seconds) between signal values detected above the threshold value (from the moment the first peak recedes below the threshold level to the moment the second peak surpasses the threshold level) that is required for each value to be considered a separate and unique peak. If the input for bktime is missing/empty the default value for the blanking time is set as the .80 quantile of the vector of time differences for signal values above the specified threshold.
#' @param plot_peaks A conditional input. If the input is TRUE or missing, an interactive plot is generated, allowing the user to manipulate the thresh and bktime values and observe the changes in peak detection. If the input is FALSE, the interactive plot is not generated. Look to the console for help on how to use the plot upon running of this function.
#' @param ... Additional inputs to be passed to FUN
#' @export
#' @return A list containing vectors for the start times, end times, peak times, peak maxima, thresh, and bktime. All times are presented as the sampling value.
#' @note As specified above under the description for the input of plot_peaks, an interactive plot can be generated, allowing the user to manipulate the thresh and bktime values and observe the changes in peak detection. The plot output is only given if the input for plot_peaks is specified as true or if the input is left missing/empty.
#' @examples \dontrun{
#' BW = beaked_whale
#' detect_peaks(data = BW$A$data, sr = BW$A$sampling_rate, FUN = njerk, thresh = NULL, bktime = NULL, plot_peaks = NULL, fs = BW$A$sampling_rate)
#' }
detect_peaks = function(data, sr, FUN = NULL, thresh = NULL, bktime = NULL, plot_peaks = NULL, ...)
{
if (missing(data) | missing(sr)) {
stop("inputs for data and sr are both required")
}
#apply function specified in the inputs to data
if (!is.null(FUN)) {
dnew = FUN(data, ...)
} else {
dnew = data
}
if ('depid' %in% names(data)) {
dnew = dnew$data
}
#set default threshold level
if (is.null(thresh) == TRUE) {
thresh = stats::quantile(dnew, 0.99)
}
if (is.null(plot_peaks) == TRUE) {
plot_peaks = TRUE
}
if (thresh > max(dnew)) {
start_time = NA
end_time = NA
peak_time = NA
peak_max = NA
thresh = thresh
if (is.null(bktime) == TRUE) {
bktime = NA
} else {
bktime = bktime
}
warning("Threshold level is greater the the maximum of the signal. No peaks are detected.")
} else {
#create matrix for time-series and corresponding sampling number
d1 = matrix(c(1:length(dnew)), ncol = 1)
d2 = matrix(dnew, ncol = 1)
d = cbind(d1, d2)
#determine peaks that are above the threshold
pt = d[, 2] >= thresh
pk = d[pt, ]
#is there more than one peak?
if (length(pk) == 2) {
start_time = pk[1]
end_time = pk[1]
peak_time = pk[1]
peak_max = pk[2]
thresh = thresh
bktime = as.numeric(bktime)
} else {
#set default blanking time
if (is.null(bktime)) {
dpk = diff(pk[, 1])
bktime = stats::quantile(dpk, c(.8))
} else {
bktime = as.numeric(bktime * sr)
}
#determine start times for each peak
dt = diff(pk[, 1])
pkst = c(1, (dt >= bktime))
start_time = pk[(pkst == 1), 1]
#determine the end times for each peak
if (sum(pkst) == 1) {
if (dnew[length(dnew)] > thresh) {
start_time = c()
end_time = c()
} else {
if (dnew[length(dnew)] <= thresh) {
end_time = pk[nrow(pk), 1]
}
}
}
if (sum(pkst) > 1) {
if (pkst[length(pkst)] == 0) {
if (dnew[length(dnew)] <= thresh) {
ending = which(pkst == 1) - 1
end_time = c(pk[ending[2:length(ending)], 1], pk[nrow(pk), 1])
} else {
if (dnew[length(dnew)] > thresh) {
ending = which(pkst == 1) - 1
end_time = c(pk[ending[2:length(ending)], 1], pk[nrow(pk), 1])
#if the last peak does not end before the end of recording, the peak is removed from analysis
start_time = start_time[1:length(start_time - 1)]
end_time = end_time[1:length(end_time - 1)]
}
}
} else {
if (pkst[length(pkst)] == 1) {
ending = which(pkst == 1) - 1
end_time = c(pk[ending[2:length(ending)], 1], pk[nrow(pk), 1])
}
}
}
#determine the time and maximum of each peak
peak_time = matrix(0, length(start_time), 1)
peak_max = matrix(0, length(start_time), 1)
if (is.null(start_time) & is.null(end_time)) {
peak_time = c()
peak_max = c()
} else {
for (a in 1:length(start_time)) {
td = dnew[start_time[a]:end_time[a]]
m = max(td)
mindex = which.max(td)
peak_time[a] = mindex + start_time[a] - 1
peak_max[a] = m
}
}
bktime = bktime / sr
}
}
#create a list of start times, end times, peak times, peak maxima, thresh, and bktime
peaks = list(start_time = start_time, end_time = end_time,
peak_time = peak_time, peak_max = peak_max, thresh = thresh, bktime = bktime)
if (plot_peaks == TRUE) {
#script for second run of detect_peaks that doesn't change the bktime from seconds to samples
detect_peaks2 = function(data, sr, thresh = NULL, bktime = NULL) {
if (thresh > max(dnew)) {
start_time = NA
end_time = NA
peak_time = NA
peak_max = NA
thresh = thresh
warning("Threshold level is greater the the maximum of the signal. No peaks are detected.")
} else {
#create matrix for time-series and corresponding sampling number
d1 = matrix(c(1:length(dnew)), ncol = 1)
d2 = matrix(dnew, ncol = 1)
d = cbind(d1, d2)
#determine peaks that are above the threshold
pt = d[, 2] >= thresh
pk = d[pt, ]
#is there more than one peak?
if (length(pk) == 2) {
start_time = pk[1]
end_time = pk[1]
peak_time = pk[1]
peak_max = pk[2]
thresh = thresh
bktime = as.numeric(bktime)
} else {
#blanking time is already known from interactive
#determine start times for each peak
dt = diff(pk[, 1])
pkst = c(1, (dt >= bktime))
start_time = pk[(pkst == 1), 1]
#determine the end times for each peak
if (sum(pkst) == 1) {
if (dnew[length(dnew)] > thresh) {
start_time = c()
end_time = c()
} else {
if (dnew[length(dnew)] <= thresh) {
end_time = pk[nrow(pk), 1]
}
}
}
if (sum(pkst) > 1) {
if (pkst[length(pkst)] == 0) {
if (dnew[length(dnew)] <= thresh) {
ending = which(pkst == 1) - 1
end_time = c(pk[ending[2:length(ending)], 1], pk[nrow(pk), 1])
} else {
if (dnew[length(dnew)] > thresh) {
ending = which(pkst == 1) - 1
end_time = c(pk[ending[2:length(ending)], 1], pk[nrow(pk), 1])
#if the last peak does not end before the end of recording, the peak is removed from analysis
start_time = start_time[1:length(start_time - 1)]
end_time = end_time[1:length(end_time - 1)]
}
}
} else {
if (pkst[length(pkst)] == 1) {
ending = which(pkst == 1) - 1
end_time = c(pk[ending[2:length(ending)], 1], pk[nrow(pk), 1])
}
}
}
#determine the time and maximum of each peak
peak_time = matrix(0, length(start_time), 1)
peak_max = matrix(0, length(start_time), 1)
if (is.null(start_time) & is.null(end_time)) {
peak_time = c()
peak_max = c()
} else {
for (a in 1:length(start_time)) {
td = dnew[start_time[a]:end_time[a]]
m = max(td)
mindex = which.max(td)
peak_time[a] = mindex + start_time[a] - 1
peak_max[a] = m
}
}
}
}
#create a list of start times, end times, peak times, peak maxima, thresh, and bktime
peaks = list(start_time = start_time, end_time = end_time, peak_time = peak_time, peak_max = peak_max, thresh = thresh, bktime = bktime)
graphics::plot(dnew, type = "l", col = "blue", xlim = c(0, length(dnew)), ylim = c(0, max(dnew)), ylab = "Signal Power", xlab = "Time (1/sampling_rate)")
x = peaks$peak_time
y = peaks$peak_max
graphics::par(new = TRUE)
graphics::plot(x, y, pch = 9, type = "p", col = "orange", xlim = c(0, length(dnew)), ylim = c(0, max(dnew)), cex = .75, ylab = "Signal Power", xlab = "Time (1/sampling_rate)")
graphics::abline(a = thresh, b = 0, col = "red", lty=2)
return(peaks)
}
#create a plot which allows for the thresh and bktime to be manipulated
graphics::plot(dnew, type = "l", col = "blue", xlim = c(0, length(dnew)), ylim = c(0, max(dnew)), ylab = "Signal Power", xlab = "Time (1/sampling_rate)")
print('GRAPH HELP:')
print('For changing only the thresh level, click once within the plot and then push enter')
print(' to specify the y-value at which your new thresh level will be.')
print('For changing just the bktime value, click twice within the plot and then push enter')
print(' to specify the length for which your bktime will be.')
print('To change both the bktime and the thresh, click three times within the plot:')
print(' the first click will change the thresh level,')
print(' the second and third clicks will change the bktime.')
print('To return your results without changing the thresh and bktime from their default')
print(' values, simply push enter.')
x = peaks$peak_time
y = peaks$peak_max
graphics::par(new = TRUE)
graphics::plot(x, y, pch = 9, type = "p", col = "orange", xlim = c(0, length(dnew)), ylim = c(0, max(dnew)), cex = .75, ylab = "Signal Power", xlab = "Time (1/sampling_rate)")
graphics::abline(a = thresh, b = 0, col = "red", lty=2)
pts = graphics::locator(n = 3)
if (length(pts$x) == 3) {
thresh = pts$y[1]
bktime = max(pts$x[2:3]) - min(pts$x[2:3])
peaks = detect_peaks2(dnew, sr, thresh = thresh, bktime = bktime)
} else {
if (length(pts$x) == 1) {
thresh = pts$y[1]
peaks = detect_peaks2(dnew, sr, thresh = thresh, bktime = bktime)
} else {
if (length(pts$x) == 2) {
bktime = max(pts$x) - min(pts$x)
peaks = detect_peaks2(dnew, sr, thresh = thresh, bktime = bktime)
} else {
peaks = detect_peaks2(dnew, sr, thresh = thresh, bktime = bktime)
}
}
}
} else {
# graphics::plot(dnew, type = "l", col = "blue", xlim = c(0, length(dnew)), ylim = c(0, max(dnew)), ylab = "Signal Power", xlab = "Time (1/sampling_rate)")
# x = peaks$peak_time
# y = peaks$peak_max
# graphics::par(new = TRUE)
# graphics::plot(x, y, pch = 9, type = "p", col = "orange", xlim = c(0, length(dnew)), ylim = c(0, max(dnew)), cex = .75, ylab = "Signal Power", xlab = "Time (1/sampling_rate)")
# graphics::abline(a = thresh, b = 0, col = "red", lty=2)
}
return(peaks)
}