信号识别-波峰波谷二阶差分识别算法
- 前言
- 波峰波谷算法
- 实际上代码
前言
在图像分析里,投影曲线是我们经常要用到的一个图像特征,通过投影曲线我们可以看到在某一个方向上,图像灰度变化的规律,这在图像分割,文字提取方面应用比较广。一个投影曲线,它的关键信息就在于波峰与波谷,所以我们面临的第一个问题就是找到波峰与波谷。
第一次涉及到求波峰与波谷时,很多人都不以为意,觉得波谷波峰还不容易,无非是一些曲线变化为零的点,从离散的角度来说,也就是:
波峰:
波谷:
这么简单吗?显示不是,你首先就会遇到这样的曲线图,然后图上的波峰点并不满足上面的条件。
看到这种情况,你也许会考虑在上面的等式中把>>和<<改为≥≥和≤≤。波峰:
&&
&&
或者
&&
&&
波谷:
&&
&&F(x)<F(x+1)$
或者
&&
&&
这次是否就这样简单,答案显示不是,下面的这个图就会让你对一些非峰值点作出错误的判断。
上面这幅图真正的峰值只有一个,其他平台上的点,你如果按上面修改的公式,就会被错误的当成波峰点。
下面让我们看一下,到底如何能求得准确的曲线波峰与波谷。
波峰波谷算法
投影曲线实际上是一个一维的向量:
其中代表图像在第ii行或列上的灰度累积。当然不仅仅是投影曲线,T也可以是某一事件中变量的观测值,我们需要研究这个变量的变化规律。
1,假投影曲线可以表示为。
2,计算V的一阶差分向量:
3,对差分向量进行取符号函数运算,,即遍历
,若
大于0,则取1;如果小于0,则取-1,否则则值为0。
4,从尾部遍历TrendTrend向量,进行如下操作:
5,对向量进行一阶差分运算,如同步骤2,得到
。
6,遍历得到的差分向量RR,如果,则
为投影向量
的一个峰值位,对应的峰值为
;如果
,则
为投影向量
的一个波谷位,对应的波谷为
。下面我们来结合一个实际的向量值,给中中间结合的计算。
1,。
它的曲线图像如下把示,图中红色圈标出了曲线的峰值,而绿字圈标出了图像的波谷位置。
2,计算的一阶差分,我们得到
。
3,对进行取符号运算,得到向量
。
4,对作一次遍历,如步骤4。
。
5,对做一阶差分,得到向量
6,遍历向量,我们就得到了两个峰值点和一个波谷点。
实际上代码
我们用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;
}
运行结果如下: