我同意 whuber 的回应，但只是想添加代码的“+2”部分，它试图移动索引以匹配新发现的峰值实际上“过冲”并且应该是“+1”。例如在手头的例子中，我们得到：

> findPeaks(cc)
[1]  3 22 41 59 78 96

当我们在图表上突出显示这些发现的峰值时（红色粗体）：

我们看到它们始终与实际峰值相差 1 点。

结果

pks[x[pks - 1] - x[pks] > thresh]

应该是pks[x[pks] - x[pks + 1] > thresh]或pks[x[pks] - x[pks - 1] > thresh]

大更新

在我自己寻求找到足够的峰值查找功能之后，我写了这个：

find_peaks <- function (x, m = 3){
    shape <- diff(sign(diff(x, na.pad = FALSE)))
    pks <- sapply(which(shape < 0), FUN = function(i){
       z <- i - m + 1
       z <- ifelse(z > 0, z, 1)
       w <- i + m + 1
       w <- ifelse(w < length(x), w, length(x))
       if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
    })
     pks <- unlist(pks)
     pks
}

“峰值”被定义为局部最大值，m其两侧的点都小于它。因此，参数 越大m，峰值筹资程序越严格。所以：

find_peaks(cc, m = 1)
[1]  2 21 40 58 77 95

该函数还可用于通过 找到任何序列向量的局部x最小值find_peaks(-x)。

注意：如果有人需要，我现在已将功能放在 gitHub 上：https ://github.com/stas-g/findPeaks

此代码的来源是通过在 R 提示符下键入其名称来获得的。输出是

function (x, thresh = 0) 
{
    pks <- which(diff(sign(diff(x, na.pad = FALSE)), na.pad = FALSE) < 0) + 2
    if (!missing(thresh)) {
        pks[x[pks - 1] - x[pks] > thresh]
    }
    else pks
}

该测试x[pks - 1] - x[pks] > thresh将每个峰值与系列中紧随其后的值（而不是系列中的下一个低谷）进行比较。它在峰值之后立即使用函数斜率大小的（粗略）估计，并仅选择斜率超过thresh大小的那些峰值。在您的情况下，只有前三个峰足够尖锐以通过测试。您将使用默认值检测所有峰值：

> findPeaks(cc)
[1]  3 22 41 59 78 96

Eek：小更新。我必须更改两行代码，即边界（添加 -1 和 +1）以达到与 Stas_G 函数的等效性（它在实际数据集中发现了太多的“额外峰值”）。为任何人道歉，我的原始帖子使我误入歧途。

我已经使用 Stas_g 的 find peaks 算法有一段时间了。由于它的简单性，它对我后来的一个项目很有帮助。然而，我需要在计算中使用它数百万次，所以我在 Rcpp 中重写了它（参见 Rcpp 包）。在简单的测试中，它比 R 版本快大约 6 倍。如果有人感兴趣，我在下面添加了代码。希望我能帮助别人，干杯！

一些小警告。此函数以 R 代码的相反顺序返回峰值索引。它需要一个内部 C++ Sign 函数，我包括在内。它尚未完全优化，但预计不会有任何进一步的性能提升。

//This function returns the sign of a given real valued double.
// [[Rcpp::export]]
double signDblCPP (double x){
  double ret = 0;
  if(x > 0){ret = 1;}
  if(x < 0){ret = -1;}
  return(ret);
}

//Tested to be 6x faster(37 us vs 207 us). This operation is done from 200x per layer
//Original R function by Stas_G
// [[Rcpp::export]]
NumericVector findPeaksCPP( NumericVector vY, int m = 3) {
  int sze = vY.size();
  int i = 0;//generic iterator
  int q = 0;//second generic iterator

  int lb = 0;//left bound
  int rb = 0;//right bound

  bool isGreatest = true;//flag to state whether current index is greatest known value

  NumericVector ret(1);
  int pksFound = 0;

  for(i = 0; i < (sze-2); ++i){
    //Find all regions with negative laplacian between neighbors
    //following expression is identical to diff(sign(diff(xV, na.pad = FALSE)))
    if(signDblCPP( vY(i + 2)  - vY( i + 1 ) ) - signDblCPP( vY( i + 1 )  - vY( i ) ) < 0){
      //Now assess all regions with negative laplacian between neighbors...
      lb = i - m - 1;// define left bound of vector
      if(lb < 0){lb = 0;}//ensure our neighbor comparison is bounded by vector length
      rb = i + m + 1;// define right bound of vector
      if(rb >= (sze-2)){rb = (sze-3);}//ensure our neighbor comparison is bounded by vector length
      //Scan through loop and ensure that the neighbors are smaller in magnitude
      for(q = lb; q < rb; ++q){
        if(vY(q) > vY(i+1)){ isGreatest = false; }
      }

      //We have found a peak by our criterion
      if(isGreatest){
        if(pksFound > 0){//Check vector size.
         ret.insert( 0, double(i + 2) );
       }else{
         ret(0) = double(i + 2);
        }
        pksFound = pksFound + 1;
      }else{ // we did not find a peak, reset location is peak max flag.
        isGreatest = true;
      }//End if found peak
    }//End if laplace condition
  }//End loop
  return(ret);
}//End Fn

首先：该算法还错误地调用平坦高原右侧的下降，因为sign(diff(x, na.pad = FALSE)) 它将是 0 然后 -1，因此它的差异也将是 -1。一个简单的解决方法是确保负输入之前的符号差异不是零而是正：

    n <- length(x)
    dx.1 <- sign(diff(x, na.pad = FALSE))
    pks <- which(diff(dx.1, na.pad = FALSE) < 0 & dx.1[-(n-1)] > 0) + 1

第二：该算法给出了非常局部的结果，例如，在序列中三个连续项的任何运行中，一个“向上”，然后是一个“向下”。如果有人对噪声连续函数的局部最大值感兴趣，那么 - 可能还有其他更好的东西，但这是我的廉价而直接的解决方案

1. 首先使用 3 个连续点的运行平均值来识别峰值，以
   使数据变得如此轻微。也采用上述控制来防止平坦然后下降。

2. 对于黄土平滑版本，通过比较以每个峰值为中心的窗口内的平均值与外部局部项的平​​均值来过滤这些候选。

   "myfindPeaks" <- 
   function (x, thresh=0.05, span=0.25, lspan=0.05, noisey=TRUE)
   {
     n <- length(x)
     y <- x
     mu.y.loc <- y
     if(noisey)
     {
       mu.y.loc <- (x[1:(n-2)] + x[2:(n-1)] + x[3:n])/3
       mu.y.loc <- c(mu.y.loc[1], mu.y.loc, mu.y.loc[n-2])
     }
     y.loess <- loess(x~I(1:n), span=span)
     y <- y.loess[[2]]
     sig.y <- var(y.loess$resid, na.rm=TRUE)^0.5
     DX.1 <- sign(diff(mu.y.loc, na.pad = FALSE))
     pks <- which(diff(DX.1, na.pad = FALSE) < 0 & DX.1[-(n-1)] > 0) + 1
     out <- pks
     if(noisey)
     {
       n.w <- floor(lspan*n/2)
       out <- NULL
       for(pk in pks)
       {
         inner <- (pk-n.w):(pk+n.w)
         outer <- c((pk-2*n.w):(pk-n.w),(pk+2*n.w):(pk+n.w))
         mu.y.outer <- mean(y[outer])
         if(!is.na(mu.y.outer)) 
           if (mean(y[inner])-mu.y.outer > thresh*sig.y) out <- c(out, pk)
       }
     }
     out
   }