一种寻找曲线峰值并统计峰的个数的python代码

前些时间我需要用python写一段代码来判别一个曲线图中有多少个峰，曲线图类型大概如下图所示：图1原始曲线图如果用肉眼鉴别，很容易就知道图中峰的个数为2。然而，计算机不是人脑的思维方式，它靠的是数据和算法，因此需要将形象的东西数值化并将判别公式化。当时我在百度以各种关键词查询如何写一段python程序来判别峰值并进行个数的计算。然而都没有好的方法。最简单粗暴的回答是，寻找这个曲线的最大值和最小值，显然，这个方法是行不通的，从上面的图可以看到，第二个峰并不是最大值，但它确实是一个峰。图2寻找最大值和最小值确定波峰和波谷然后我又查找了一下，有这么个方法，所谓峰值就是前后两个值都比当前的值小，这样子的值便可以定义为峰值，并设计了代码如下：图3寻找极值作为峰值这方法看上去很合理，可是在不平滑的曲线上会接连遇到问题，首选我们看一下图1，红色箭头指示的两个位置，一个是噪声（箭头1），它符合前后的值都小于它的条件，会被判别为峰，而另一个（箭头2），由于最高峰位置是个平顶，不符合前后都大于它，会被认为不是峰，这些均导致了上述方法也行不通！图4有噪声不规则曲线判别存在的问题为了消除噪声被误判为峰，我们需要设定一个阈值，大于该阈值的才会被确认为峰值，避免噪声也被判为峰。图5阈值设置而对于平顶的情况，处理起来比较复杂，我们先一步一步说明。除了平顶的情况，按照上面的判别算法，即便设置了阈值，如果在一个峰的顶部有尖峰，也会被判为一个峰，如果有几个尖峰，就一个就变成了几个。图6尖峰导致错判面对平顶或者有尖峰的情况，最好的办法就是将粗糙的曲线进行平滑。下图可以看到，经过平滑后情况便好了很多。图7进行曲线平滑的情况基本平滑以后可以防止上述错误的出现，但为了保险起见（主要防止平顶和尖峰在平滑后还是存在），我们可以认为大于阈值的比前后都大的值以及与前或后都值相同的都是峰，然后设置峰的宽度，目的是保证在这个宽度内的所有尖峰都视同是同一个峰。总的来说，就是先进行曲线平滑，然后设置阈值和峰宽度，然后寻找符合阈值的极值和平顶值，在根据峰的宽度看这些值是否在宽度内，如果在宽度内将被认为是同一个峰。我试了一下基本都可以。以下是程序（注意：输入的曲线我当时是用了已经归一化的曲线，所以值是在0-1范围）：###对信号(signal)进行均值滤波，滤波窗口大小为window_size###要求信号(signal)是一个列表(list)defgen(l,window_size):index=0ans=0times=0whileTrue:whileindexwindow_size:ans+=l[times+index]index+=1yieldfloat(ans)/float(window_size)###Resetindex=0ans=0times+=1defmean_filter(signal,window_size):window_size=8temp=gen(signal,window_size)filtered=[]foriinrange(len(signal)-window_size):filtered.append(next(temp))returnfiltereddeffind_peak(filtered_signal,length_data,thre,peak_width):l=[];foriinrange(1,length_data-1):#在整个B的长度内找出极值iffiltered_signal[i-1]filtered_signal[i]andfiltered_signal[i]filtered_signal[i+1]andfiltered_signal[i]thre:#找出极值，并设置阈值，这里阈值设为20l.append(i);#找出极值的位置eliffiltered_signal[i]==filtered_signal[i-1]andfiltered_signal[i]thre:l.append(i);#最高点前后可能有相等的情况CC=len(l)#统计极值有几个，如果有两条线，就会有两个，如果有一条线就只有一个#print(CC)#print(l)cou=0forjinrange(1,CC):ifl[j]-l[j-1]peak_width:#此判断用于将位于同一个峰内的极值点去除cou=cou+1rcou=CC-coureturnrcoudefmain():signal=Dwindow_size=8filtered_signal=mean_filter(signal,window_size)plt.plot(filtered_signal)plt.show()length_data=len(filtered_signal)thre=0.3#峰的阈值高度，只有高于此值才算一个峰，否则算是噪声波动，建议这个值在均值之上peak_width=20#为了避免一个峰内出现几个尖峰而导致被判别为多个峰，这里设定峰的宽度，在此宽度内均认为属于同一个峰rcou=find_peak(filtered_signal,length_data,thre,peak_width)print(rcou)if__name__==__main__:main()