-
Notifications
You must be signed in to change notification settings - Fork 1
/
find_peaks.R
96 lines (92 loc) · 3.25 KB
/
find_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
#' A simple peak finding algorithm for ESR spectra
#'
#' This is a simple algorithm to determine local minima and maxima of ESR spectra, but
#' can also be used for various other kind of data.
#'
#' The algorithm characterises a peak by comparing the y-value at each x-value to its neighbouring y-values.
#' The threshold value specified by \code{th} determines to how many neighbours the value
#' is compared to, so that lower values of \code{th} yield more peaks than higher
#' threshold values.
#'
#' @param x \code{\link{data.frame}} (required): a two column \link{data.frame}
#' @param interval \code{\link{numeric}}: a vector of length two specifying
#' the range of x-values where peaks are searched
#' @param th \code{\link{numeric}}: an integer specifying the number of neighbouring values
#' to compare each x-value to
#' @return A \link{data.frame} containing the x- and y-values of each peak is returned.
#' @export
#' @author Christoph Burow, University of Cologne (Germany)
#' @examples
#'
#' # Import Bruker ELEXSYS500 spectrum (ASCII)
#' file <- system.file("extdata", "dpph.ASC", package = "ESR")
#' spec <- read_Spectrum(file)
#'
#' # Use the function
#' peaks1 <- find_Peaks(spec)
#'
#' # Since spec is an object of class ESR.Spectrum, you can also use
#' # its get-method
#' peaks2 <- spec$get_peaks()
#'
#' identical(peaks1, peaks2)
#'
#' @export find_Peaks
find_Peaks <- function(x, interval, th = 10) {
if (!inherits(x, "data.frame") && !inherits(x, "data.frame") && !inherits(x, "R6"))
stop("x must be a data.frame", call. = FALSE)
if (inherits(x, "R6") && inherits(x, "ESR.Spectrum"))
x <- as.data.frame(x)
if (inherits(x, "ESR.Spectrum"))
class(x) <- class(x)[which(class(x) != "ESR.Spectrum")]
if (ncol(x) != 2)
stop("x must have two columns.", call. = FALSE)
if (missing(interval))
interval <- range(x[, 1])
## Preparation
scan.width <- seq(from = 1, to = length(x[, 1]), by = 1)
peak.max.storage <- peak.min.storage <- matrix(nrow = length(scan.width), ncol = 2)
## FIND PEAKS
for (i in 1:length(x[, 1])) {
# find max peaks
if (any(abs(x[, 2][i:c(i + if (i + th > length(x[, 1])) {
length(x[, 1]) - i
} else {
th
})]) > abs(x[, 2][i])) == FALSE) {
if (any(abs(x[, 2][c(i - if (i < th) {
i - 1
} else {
th
}):i]) > abs(x[, 2][i])) == TRUE) {
} else {
if (x[, 1][i] > interval[1] && x[, 1][i] < interval[2]) {
peak.max.storage[i, ] <- as.matrix(c(x[, 1][i],
x[, 2][i]))
}
}
}
# find min peaks
if (any(abs(x[, 2][i:c(i + if (i + th > length(x[, 1])) {
length(x[, 1]) - i
} else {
th
})]) < abs(x[, 2][i])) == FALSE) {
if (any(abs(x[, 2][c(i - if (i < th) {
i - 1
} else {
th
}):i]) < abs(x[, 2][i])) == TRUE) {
} else {
if (x[, 1][i] > interval[1] && x[, 1][i] < interval[2]) {
peak.min.storage[i, ] <- as.matrix(c(x[, 1][i],
x[, 2][i]))
}
}
}
}
all.peaks <- as.data.frame(rbind(na.omit(peak.max.storage), na.omit(peak.min.storage)))
setorder(all.peaks)
setnames(all.peaks, c("magnetic.field", "ESR.intensity"))
return(all.peaks)
}