信号识别-波峰波谷二阶差分识别算法

---

前言

在图像分析里，投影曲线是我们经常要用到的一个图像特征，通过投影曲线我们可以看到在某一个方向上，图像灰度变化的规律，这在图像分割，文字提取方面应用比较广。一个投影曲线，它的关键信息就在于波峰与波谷，所以我们面临的第一个问题就是找到波峰与波谷。
第一次涉及到求波峰与波谷时，很多人都不以为意，觉得波谷波峰还不容易，无非是一些曲线变化为零的点，从离散的角度来说，也就是:

波峰：$F(x)>F(x-1)且F(x)>F(x+1)$
波谷：$F(x)<F(x-1)且F(x)<F(x+1)$

这么简单吗？显示不是，你首先就会遇到这样的曲线图，然后图上的波峰点并不满足上面的条件。

看到这种情况，你也许会考虑在上面的等式中把>>和<<改为≥≥和≤≤。

波峰：
$F(x)\ge F(x-1)$&&$F(x)>F(x+1)F(x)\ge F(x-1)$&&$F(x)>F(x+1)$
或者
$F(x)>F(x-1)$&&$F(x)\ge F(x+1)F(x)>F(x-1)$&&$F(x)\ge F(x+1)$
波谷：
$F(x)\le F(x-1)$&&$F(x)<F(x+1)F(x)\le F(x-1)$&&F(x)<F(x+1)$
或者
$F(x)<F(x-1)$&&$F(x)\le F(x+1)F(x)<F(x-1)$&&$F(x)\le F(x+1)$
这次是否就这样简单，答案显示不是，下面的这个图就会让你对一些非峰值点作出错误的判断。

上面这幅图真正的峰值只有一个，其他平台上的点，你如果按上面修改的公式，就会被错误的当成波峰点。

下面让我们看一下，到底如何能求得准确的曲线波峰与波谷。

---

波峰波谷算法

投影曲线实际上是一个一维的向量：
$V=[v1,v2,\dots ,vn]$
其中$vi,i\in [1,2,\dots ,N]$代表图像在第ii行或列上的灰度累积。当然不仅仅是投影曲线，T也可以是某一事件中变量的观测值，我们需要研究这个变量的变化规律。
1，假投影曲线可以表示为$V=[v1,v2,\dots ,vn]$。

2，计算V的一阶差分向量$DiffV$:

$Diffv(i)=V(i+1)-V(i),其中i\in 1,2,\dots ,N-1$
3,对差分向量进行取符号函数运算，$Trend=sign(Diffv)$,即遍历$Diffv$，若$Diffv(i)$大于0，则取1；如果小于0，则取-1，否则则值为0。

4，从尾部遍历TrendTrend向量,进行如下操作：
$ifTrend(i)=0且Trend(i+1)\ge 0，则Trend(i)=1$
$ifTrend(i)=0且Trend(i+1)<0，则Trend(i)=-1$
5，对$Trend$向量进行一阶差分运算，如同步骤2，得到$R=diff(Trend)$。
6，遍历得到的差分向量RR，如果$R(i)=-2$，则$i+1$为投影向量$V$的一个峰值位，对应的峰值为$V(i+1)$；如果$R(i)=2$，则$i+1$为投影向量$V$的一个波谷位，对应的波谷为$V(i+1)$。

下面我们来结合一个实际的向量值，给中中间结合的计算。
1，$V=[-5,10,10,14,14,8,8,6,6,-3,2,2,2,2,-3]$。
它的曲线图像如下把示，图中红色圈标出了曲线的峰值，而绿字圈标出了图像的波谷位置。

2，计算$V$的一阶差分，我们得到$Diff(V)=[15,0,4,0,,-6,0,-2,0,-9,5,0,0,0,-5]$。

3，对$Diffv$进行取符号运算，得到向量$Trend=[1,0,1,0,-1,0,-1,0,-1,1,0,0,0,-1]$。

4，对$Trend$作一次遍历，如步骤4。$Trend=[1,1,1,-1,-1,-1,-1,-1,-1,1,-1,-1,-1,-1]$。

5，对$Trend$做一阶差分，得到向量$R=Diff(Trend)=[0,0,-2,,0,0,0,0,0,2,-2,0,0,0]$

6，遍历向量$R$,我们就得到了两个峰值点和一个波谷点。

实际上代码

我们用matlab生成一个100个数组的波形如下

matlab随机生成的代码如下。

x = [];
y = [];
for c=1:100
    x(c)=c;
    y(c)=rand*32000;
    
end
disp(x) 
format shortG
disp(y) 
for c=1:100
    x(c)=c;
    y(c)=rand*32000;
    X =sprintf('%f,',y(c));
    disp(X) 
end
plot(x,y)

把打印出来的随机的y左边的代码放到我们的程序中去

#include <stdio.h>
#include <stdlib.h>
#define SAMPLE_MAX  1300

#define PV_MAX      1300

//float Sample[SAMPLE_MAX]={1,2,3,4,4,4,5,2,1,0,0,5,1,0,0,1,2,3,4,0};
float Sample[SAMPLE_MAX]={
13648.136043,
2439.662155,
9298.772911,
17962.721514,
20266.670986,
29784.827490,
31288.613046,
2995.089472,
21175.518221,
19288.732784,
15162.171623,
11400.202960,
15218.510132,
21472.704585,
30708.646304,
2850.696496,
25527.758715,
18904.823877,
29190.296136,
3236.138204,
9385.428575,
1650.822572,
16132.101779,
24588.024061,
9055.532039,
7211.527865,
10601.268564,
14504.043245,
23596.308422,
16316.340951,
12240.454139,
28975.462939,
30888.247967,
20104.555282,
4224.997442,
19785.656907,
12256.643330,
31718.199922,
9178.473112,
22598.125747,
17126.588720,
6182.753235,
22061.923230,
1614.560101,
5901.891702,
1461.066351,
28321.328065,
26873.420346,
3780.967845,
13133.275318,
3847.314880,
18306.960712,
30380.484308,
8204.325114,
31675.693925,
11193.857313,
6672.715965,
21306.469121,
31147.030615,
19926.394289,
2033.212562,
11952.305064,
5320.049134,
7400.898499,
1670.676265,
28856.202135,
25385.334397,
11936.458395,
26625.750572,
24122.705811,
19899.619112,
12610.981787,
11496.894766,
2843.278749,
10933.664116,
17557.459010,
14737.518755,
20654.475950,
16432.666921,
26061.646994,
3109.862592,
14838.833941,
18874.159667,
5989.506257,
19562.564692,
1662.143404,
18423.262351,
26955.030459,
15991.231156,
14048.796405,
4769.826664,
904.942429,
24213.433215,
25475.399550,
9393.779021,
3686.615873,
12002.932192,
26524.599842,
26936.850623,
21287.630450,
30724.479228,
30179.772551,
3606.383756,
20745.198375,
15385.729654,
2128.662142,
28728.679217,
15911.364126,
24681.709723,
1931.591342,
8398.638877,
20834.219798,
4275.325134,
20433.463578,
12318.179039,
24502.340022,
20893.301829,
12207.650089,
9600.594648,
10884.465776,
29405.649175,
14600.538026,
14159.902483,
14533.946228,
30249.024294,
7011.806482,
28236.889470,
636.013530,
10936.476264,
24512.878262,
10969.718449,
19801.804933,
14496.677187,
325.203284,
19170.597359,
19250.187147,
20781.358431,
10967.060243,
15785.576178,
22456.766248,
28409.683246,
1761.851077,
3147.582417,
20793.049135,
24450.268520,
31614.686131,};
  
float SampleDiff[SAMPLE_MAX]={0};

typedef struct _tag_FindPV
{
int Pos_Peak[PV_MAX];//波峰位置存储
int Pos_Valley[PV_MAX];//波谷位置存储
int Pcnt;//所识别的波峰计数
int Vcnt;//所识别的波谷计数
}SFindPV;

 
SFindPV stFindPV;
 
 /********************************************
  *  Fuction : initialFindPV
 *  Note    : 初始化相关数据 
  *******************************************/
void initialFindPV(void)
{
	int Index=0;
	
	for(Index=0;Index<SAMPLE_MAX;Index++)
	{
	SampleDiff[Index]=0;
	}
	
	for(Index=0;Index<PV_MAX;Index++)
	{
	stFindPV.Pos_Peak[Index]=-1;
	stFindPV.Pos_Valley[Index]=-1;
	}
	stFindPV.Pcnt=0;
	stFindPV.Vcnt = 0;
 }

//找波峰波谷 
void FindPV(SFindPV *pFindPV,float*Sample)
 {
	int i=0;
	//step 1 :首先进行前向差分，并归一化
	for(i=0;i<SAMPLE_MAX-1;i++)
	{
		int samplei1=Sample[i+1]/1000;
		int samplei=Sample[i]/1000;
		//printf("%d   %d \n",samplei1,samplei);
		if(samplei1-samplei>0)
		SampleDiff[i]=1;
		else if(samplei1-samplei<0)
		SampleDiff[i]=-1;
		else
		SampleDiff[i]=0;
	}

	//step 2 :对相邻相等的点进行领边坡度处理
	for(i=0;i<SAMPLE_MAX-1;i++)
	{
		if(SampleDiff[i]==0)
		{
			if(i==(SAMPLE_MAX-2))
			{

				if (SampleDiff[i-1]>=0)
				{

			
					SampleDiff[i]=1;	}
				else	{
					SampleDiff[i]=-1;	}
			}
			else
			{
				if(SampleDiff[i+1]>=0)
				SampleDiff[i]=1;
				else
				SampleDiff[i]=-1;
			}
		}
	}

	//step 3 :对相邻相等的点进行领边坡度处理
	for(i=0;i<SAMPLE_MAX-1;i++)
	{
	if(SampleDiff[i+1]-SampleDiff[i]==-2)//波峰识别
	{
	pFindPV->Pos_Peak[pFindPV->Pcnt]=i+1;
	pFindPV->Pcnt++;
	}
	else if(SampleDiff[i+1]-SampleDiff[i]==2)//波谷识别
	{
	pFindPV->Pos_Valley[pFindPV->Vcnt]=i+1;
	pFindPV->Vcnt++;
	}
}
}

/********************************************
 *  Fuction : main
 *  Note    : 模拟查找波峰波谷 
 *******************************************/
int main(int argc, char *argv[]) {
int i=0;
initialFindPV();
FindPV(&stFindPV,Sample);

printf("Peak cnt %d \n",stFindPV.Pcnt);
for(i=0;i<stFindPV.Pcnt;i++)
{
printf("-%d",stFindPV.Pos_Peak[i]+1);//加1是为了与上图横坐标一致 
}

printf("\nValley cnt %d \n",stFindPV.Vcnt);
for(i=0;i<stFindPV.Vcnt;i++)
{
printf("-%d",stFindPV.Pos_Valley[i]+1);
}

printf("\n\n");
return 0;
}

运行结果如下：

这就是找到的所有的波峰波谷！！！
大功告成！！！

代码地址
本文的代码地址