Maybe it's already out there in the R, but I haven't seen it before.

is.localmax = function(i,context, k, FUN=max, ...){

context = c(rep(NA,k), context, rep(NA,k)) # pad with NAs

i = i + k # account for the NA padding

contextMatrix = sapply(i, function(x) {context[(x-k):(x+k)]})

context[i] == apply(contextMatrix, 2, FUN, ...)

}getPeaks = function(x,threshold,k) {

i = which(x>threshold)

if (length(i)>0)

i[is.localmax(i,context=x, k=k,FUN=max)]

else

NA

}

Basically, it finds all the indices above some threshold (quickly), and then only checks those to see if they are local maxima for some local window of length 2*k+1.

NOTE: this function fails if local maxima aren't unique. If you're likely to get two local maxima (e.g. from a low-resolution ADC) in a single peak (e.g. the vector [0 0 1 2 3 3 2 1 0 ], it will return both indices 5 and 6), which may be undesirable. To work around this, add a very very small amount of noise to the signal (x) to force the code to select one of the indices at random.