请教 findPeaks() 函数的使用问题

caiwxl

用findPeaks（data）找出来的峰值好多都不对，请问是是什么原因？help里的说明挺少的

superdesolator

回复第1楼的 caiwxl：简单搜了下这个函数Matlab的吧，R怎么找这个函数，如果是一维找峰值自己实现下应该很简单，如果是2维3维的话就得给一个找峰值的算法，然后自己实现也不难，关键多维的时候峰值怎么定义？

superdesolator

回复第1楼的 caiwxl：我写了个丑陋的函数，输入为矩阵，峰值方式为cross,就是横竖4个值加中心，看中心是不是最大，边缘就是有多少数据考虑多少

findPeaks=function(data){<br />
  findoned=function(x){<br />
    xtmp=c(-Inf,x,-Inf)<br />
    tmp=sapply(seq_along(x)+1, function(y)   ifelse(xtmp[y-1]<xtmp[y]&&xtmp[y+1]<xtmp[y],1,-1))<br />
   return(which(tmp==1))<br />
  }<br />
  a=apply(data,1,findoned)<br />
  rowpos1=unlist(lapply(seq_along(a),function(x) rep(x,length(a[[x]]))))<br />
  colpos1=unlist(a)<br />
  a=apply(data,2,findoned)<br />
  colpos2=unlist(lapply(seq_along(a),function(x) rep(x,length(a[[x]]))))<br />
  rowpos2=unlist(a)<br />
 pos2=Map(function(x,y) c(x,y),rowpos2,colpos2)<br />
 pos1=Map(function(x,y) c(x,y),rowpos1,colpos1)<br />
 pos=intersect(pos1,pos2)<br />
 rowpos=sapply(pos,function(x) x[1])<br />
 colpos=sapply(pos,function(x) x[2])<br />
 Peaks=Map(function(x,y) data[x,y],rowpos,colpos)<br />
 result=list(Peaks=unlist(Peaks),rowpos=rowpos,colpos=colpos)<br />
 return(result)<br />
}<br />
R>data<br />
     [,1] [,2] [,3] [,4]<br />
[1,]    7   14   12   20<br />
[2,]   18    4   10   16<br />
[3,]   11   17   15    5<br />
[4,]    9    3    1    7<br />
[5,]    8   13    7   18<br />
[6,]   19    6    2   11<br />
R>findPeaks(data)<br />
$Peaks<br />
[1] 14 20 18 17 13 18 19</p>
<p>$rowpos<br />
[1] 1 1 2 3 5 5 6</p>
<p>$colpos<br />
[1] 2 4 1 2 2 4 1

</p>

superdesolator

想了下，增加了9宫格峰值的方法，并且结合了上面的代码，重新贴出

<br />
findPeaks=function(data,methods="cross"){<br />
    findoned=function(x){<br />
    xtmp=c(-Inf,x,-Inf)<br />
    tmp=sapply(seq_along(x)+1, function(y)   ifelse(xtmp[y-1]<xtmp[y]&&xtmp[y+1]<xtmp[y],1,-1))<br />
   return(which(tmp==1))<br />
  }<br />
  a=apply(data,1,findoned)<br />
  rowpos1=unlist(lapply(seq_along(a),function(x) rep(x,length(a[[x]]))))<br />
  colpos1=unlist(a)<br />
  a=apply(data,2,findoned)<br />
  colpos2=unlist(lapply(seq_along(a),function(x) rep(x,length(a[[x]]))))<br />
  rowpos2=unlist(a)<br />
 pos2=Map(function(x,y) c(x,y),rowpos2,colpos2)<br />
 pos1=Map(function(x,y) c(x,y),rowpos1,colpos1)<br />
 pos=intersect(pos1,pos2)<br />
 rowpos=sapply(pos,function(x) x[1])<br />
 colpos=sapply(pos,function(x) x[2])<br />
 Peaks=Map(function(x,y) data[x,y],rowpos,colpos)<br />
 if(methods=="cross") {<br />
 result=list(Peaks=unlist(Peaks),rowpos=rowpos,colpos=colpos)<br />
 return(result)  }<br />
 if(methods=="diagandcross"){<br />
  data1=apply(data,2,function(x) c(-Inf,x,-Inf) )<br />
  data2=apply(data1,1,function(x) c(-Inf,x,-Inf))<br />
  data3=t(data2)<br />
  resultlogical=unlist(Map(function(x,y){<br />
    tmp=data3[x,y]<br />
    return(data3[x-1,y-1]<tmp&&data3[x-1,y+1]<tmp&&data3[x+1,y-1]<tmp&&data3[x+1,y+1]<tmp)<br />
  },rowpos+1,colpos+1))<br />
  result=list(Peaks=unlist(Peaks)[resultlogical],rowpos=rowpos[resultlogical],colpos=colpos[resultlogical])<br />
  return(result)<br />
 }<br />
}<br />
R>data<br />
     [,1] [,2] [,3] [,4]<br />
[1,]    7   14   12   20<br />
[2,]   18    4   10   16<br />
[3,]   11   17   15    5<br />
[4,]    9    3    1    7<br />
[5,]    8   13    7   18<br />
[6,]   19    6    2   11<br />
R>findPeaks(data,"cross")<br />
$Peaks<br />
[1] 14 20 18 17 13 18 19</p>
<p>$rowpos<br />
[1] 1 1 2 3 5 5 6</p>
<p>$colpos<br />
[1] 2 4 1 2 2 4 1</p>
<p>R>findPeaks(data,"diagandcross")<br />
$Peaks<br />
[1] 20 18 18 19</p>
<p>$rowpos<br />
[1] 1 2 5 6</p>
<p>$colpos<br />
[1] 4 1 4 1

主要思路是：

1.写个1维寻找峰值的函数，findoned,考虑边缘的情况,就是两端

2.利用这个函数，对矩阵进行行和列的寻找峰值的位置，并且记录对应的行号和列号

3.考虑这2堆点的重合点就是cross情况的峰值

4.增加diagandcross的时候，在3返回的点中，在继续考察以它为中心的九宫格对角是不是它仍然最大，为了简化边缘的情况，我在原矩阵data外套了一圈都是-Inf的套子</p>

PS：当然，最原始的方法是：

直接在矩阵外套个-Inf的套子，然后用两层循环，来一个个判断，稍微简化的方式是判断过的，可以把对应的未判断过的设置为不用比较，因为那个九宫格中已经有峰值了，其它不可能是峰值了

superdesolator

用另两层循环写的，中间用了个简化的方式就是用一个矩阵作判断要不要继续判断的指标,牺牲一个data的内存来增加速度，结果是个两列的矩阵，第一列是横坐标第二列是纵坐标

findPeaks2=function(data,methods="cross"){<br />
  data=data+abs(min(data))+1<br />
  data1=apply(data,2,function(x) c(-Inf,x,-Inf) )<br />
  data2=apply(data1,1,function(x) c(-Inf,x,-Inf))<br />
  data3=t(data2)<br />
  result=NULL<br />
  if(methods=="cross"){<br />
   for(i in 2:(nrow(data3)-1)){<br />
     for(j in 2:(ncol(data3)-1)) {<br />
        data4=data3<br />
        tmp1=data4[i,j]<br />
        tmp=data3[i,j]<br />
       if(tmp1>0){<br />
       if(tmp>data3[i-1,j]&&tmp>data3[i,j-1]&&tmp>data3[i+1,j]&&tmp>data3[i,j+1])     {<br />
        data4[i+1,j]=-data3[i+1,j]<br />
        data4[i,j+1]=-data3[i,j+1]<br />
        result=rbind(result,c(i-1,j-1))<br />
        }<br />
        }<br />
      }<br />
   }<br />
   return(result)<br />
 }<br />
 if(methods=="diagandcross"){<br />
   for(i in 2:(nrow(data3)-1)){<br />
     for(j in 2:(ncol(data3)-1)) {<br />
       data4=data3<br />
       tmp1=data4[i,j]<br />
       tmp=data3[i,j]<br />
       if(tmp1>0){<br />
         if(tmp>data3[i-1,j]&&tmp>data3[i,j-1]&&tmp>data3[i+1,j]&&tmp>data3[i,j+1]&&tmp>data3[i-1,j-1]&&tmp>data3[i-1,j+1]&&tmp>data3[i+1,j-1]&&tmp>data3[i+1,j+1])     {<br />
           data4[i+1,j]=-data3[i+1,j]<br />
           data4[i,j+1]=-data3[i,j+1]<br />
           data4[i+1,j-1]=-data3[i+1,j-1]<br />
           data4[i+1,j+1]=-data3[i+1,j+1]<br />
           result=rbind(result,c(i-1,j-1))<br />
         }<br />
       }<br />
     }<br />
   }<br />
   return(result)<br />
 }</p>
<p>}

</p>