findpeaks: Find local maxima
In gjmvanboxtel/DSPutils: Some Useful Signal Processing Utilities in R

Description Usage Arguments Details Value Author(s) References Examples

Find local maxima (peaks) in an input signal vector.

1	findpeaks (data, minh = .Machine$double.eps, mind = 1, minw = 1, maxw = Inf, ds = FALSE)

`data`	The input signal vector.
`minh`	Minimum peak height (non-negative scalar). Only peaks that exceed this value will be returned. For data taking positive and negative values use the option `ds`. Default: machine precision (.Machine$double.eps).
`mind`	Minimum separation between peaks (positive integer). Peaks separated by less than this distance are considered a single peak. This distance is also used to fit a second order polynomial to the peaks to estimate their width, (see `details`, therefore it acts as a smoothing parameter. Default value 1.
`minw`	Minimum width of peaks (positive integer). Default value 1.
`maxw`	Maximum width of peaks (positive integer). Default value `Inf`.
`ds`	Double-sided. Tells the function that data takes positive and negative values. The base-line for the peaks is taken as the mean value of the function. This is equivalent as passing the absolute value of the data after removing the mean.

This function searches for peaks in a signal vector. Peaks of a positive array of data are defined as local maxima. For double-sided data, they are maxima of the positive part and minima of the negative part. The function provides various options to search for peaks in noisy data, such specifying a minimum peak height (minh), a minimum distance between peaks (mind), and a minimum or maximum width of the peaks (minw and maxw).

The width of the peaks is estimated using a parabola fitted to the neighborhood of each peak. The width is caulculated with the formula a * (width - x0)^2 = 1, where a is the the concavity of the parabola and x0 its vertex. The neighborhood size is equal to the value of mind.

When called with minw = 0 and maxw = Inf, this function returns a list containing two values:

`$pks:`	array containing the value of `data` at the peaks
`$idx:`	array containing the peak indices

When called with either minw > 0 or maxw < Inf, then the returned list contains these additional variables:

`$parabol:`	a `list` containing additional information about the parabol fitted to the peak. The `list` `$pp` contains the coefficients of the 2nd degree polynomial (`a`, `b`, and `b2`), and `$x` the extrema of the interval where it was fitted (`$from`, `to`).
`$height:`	The estimated height of the returned peaks (in units of `data`).
`$baseline:`	The height at which the roots of the returned peaks were calculated (in units of `data`).
`$roots:`	The abscissa values (in index units) at which the parabola fitted to each of the returned peaks realizes its width.

Geert van Boxtel

Octave Signal package

# Example 1: Finding the peaks of smooth data is not a big deal

t <- 2*pi*seq(0,1,length=1024)
y <- sin(3.14*t) + 0.5*cos(6.09*t) + 0.1*sin(10.11*t+1/6) + 0.1*sin(15.3*t+1/3)

data1 <- abs(y) # Positive values
peaks1 <- findpeaks(data1)

data2 <- y # Double-sided
peaks2 <- findpeaks(data2, ds=TRUE)
peaks3 <- findpeaks (data2, ds=TRUE, minh=0.5)

## Not run: 
op <- par(mfrow=c(1,2))
plot(t,data1,type="l", xlab="", ylab="")
points (t[peaks1$idx],peaks1$pks,col="red", pch=1)
plot(t,data2,type="l", xlab="", ylab="")
points (t[peaks2$idx],peaks2$pks,col="red", pch=1)
points (t[peaks3$idx],peaks3$pks,col="red", pch=4)
legend ("topleft", '0: >2*sd, x: >0.5', bty="n", text.col="red")
par (op)
## End(Not run)

# Example 2: Noisy data may need tuning of the parameters. In this example,
# "mind" is used as a smoother of the peaks.

t <- 2*pi*seq(0,1,length=1024)
y <- sin(3.14*t) + 0.5*cos(6.09*t) + 0.1*sin(10.11*t+1/6) + 0.1*sin(15.3*t+1/3)
data <- abs(y + 0.1*rnorm(length(y),1)); # Positive values + noise
peaks1 <- findpeaks(data, minh=1)
dt <- t[2]-t[1]
peaks2 <- findpeaks(data, minh=1, mind=round(0.5/dt))

## Not run: 
op <- par(mfrow=c(1,2))
plot(t, data, type="l", xlab="", ylab="")
points (t[peaks1$idx],peaks1$pks,col="red", pch=1)
plot(t, data, type="l", xlab="", ylab="")
points (t[peaks2$idx],peaks2$pks,col="red", pch=1)
par (op)
## End(Not run)