核心寻峰算法的原理参考Ronny,链接:投影曲线的波峰查找,
C#翻译原理代码参考sowhat4999,链接:C#翻译Matlab中findpeaks方法
前人种树,后人乘凉。感谢原作者详细的解释说明。
这里先把翻译代码贴一下(略微的修改了sowhat4999代码中的几个参数)
//调用方法
List<double> data = new List<double>{, , , , , , , , , , };
List<int> index = getPeaksIndex(trendSign(oneDiff(data)));
//第一次寻峰(基本峰距为1)算法
private double[] oneDiff(List<double> data)
{
double[] result = new double[data.Count - ];
for (int i = ; i < result.Length; i++)
{
result[i] = data[i + ] - data[i];
}
return result;
} private int[] trendSign(double[] data)
{
int[] sign = new int[data.Length];
for (int i = ; i < sign.Length; i++)
{
if (data[i] > ) sign[i] = ;
else if (data[i] == ) sign[i] = ;
else sign[i] = -;
} for (int i = sign.Length - ; i >= ; i--)
{
if (sign[i] == && i == sign.Length - )
{
sign[i] = ;
}
else if (sign[i] == )
{
if (sign[i + ] >= )
{
sign[i] = ;
}
else
{
sign[i] = -;
}
}
}
return sign;
} private List<int> getPeaksIndex(int[] diff)
{
List<int> data = new List<int>();
for (int i = ; i != diff.Length - ; i++)
{
if (diff[i + ] - diff[i] == -)
{
data.Add(i + );
}
}
return data;//相当于原数组的下标
}
以上方法并没有将峰距、边锋、峰值情况考虑在内,但已经给与我们后人一个完整的思路。
峰距情况分析:
我们可以将上述方法理解为峰距1的寻峰算法,当我们需要完成峰距为2的寻峰情况时我们需要判断
data[i]是否大于data[i+1],data[i+2],data[i-1],data[i-2]
同理按照此方法完成点数为100000,峰距为1000的寻峰,则需要进行100000的1000次方次运算,这显然需要花费大量的时间进行运算。
优化过程中,我们并不能改变峰距(即幂指数1000),但我们可以改变点数(即底数100000)的大小。从而实现运算量的降低。
以上峰距为1的寻峰方法此时已经完成判断
data[i]是否大于data[i+1],data[i-1]
并返还峰值对应的索引列
峰距为2时,我们只需要再次对索引列中内容进行判断即可(只有在峰距为1的判断中胜出的点,才有可能在峰距为2的判断中胜出)
data[i]是否大于data[i+2],data[i-2]
此时你会发现我们需要遍历的底数已经并不是原点数100000,而是上次返还的寻峰序列个数
// 调用方法
List<double> Xaxis = new List<double> { , , , , , , , , , , };
List<double> Yaxis = new List<double> { , , , , , , , , , , };
// 峰距
int DisPeak = ;
// 峰距为3时得到的脚标
List<int> index = getPeaksIndex(trendSign(oneDiff(Yaxis)));
// 已进行的判断
int level = ;
// 扩大峰距范围范围算法
while (DisPeak > level)
{
level++;
List<int> result = DoPeakInstance(Yaxis, index, level);
index = null;
index = result;
} // 获取两侧满足条件的边峰序列
index = GetBothSidePeakIndex(Xaxis, Yaxis, , index); double minFZ = 10.0;
// 根据最小峰值序列进行筛选
index = FindMinPeakValue(minFZ, Yaxis, index);
//扩大寻峰范围算法
private List<int> DoPeakInstance(List<double> data, List<int> index, int level)
{
//相当于原数组的下标
List<int> result = new List<int>();
for (int i = ; i < index.Count; i++)
{
//判断是否超出下界和上界
if (index[i] - level >= && index[i] + level < data.Count)
{
if (data[index[i] + level] <= data[index[i]] && data[index[i] - level] <= data[index[i]])
{
result.Add(index[i]);
}
}
}
return result;
}
边锋情况分析:
仔细阅读上述两算法,你会发现该算法存在一个无法避免的问题 如:
峰距是3,此时峰首部点序(点0,点1,点2)因无法向前比较,导致并没有参与到峰值计算中。 尾部点则因无法向后比较没有参与到峰值计算中。
此情况我们首先要清楚,因上述情况未参与比较的点序中,首部最多仅有一个峰值,尾部最多仅有一个峰值。
那我们把它加上就好了,美滋滋。
//获取两侧满足条件的边峰序列
private static List<int> GetBothSidePeakIndex(List<double> Xaxis, List<double> Yaxis, int FJ, List<int> index)
{
//获取数据首尾两侧最大峰值(0,FJ)点序和(Date.CountFJ-FJ,Data.Count)点序
int TopIndex = ;
int BottomIndex = Yaxis.Count - ;
for (int i = ; i < FJ; i++)
{
if (Yaxis[i] >= Yaxis[TopIndex])
{
TopIndex = i;
}
if (Yaxis[Yaxis.Count - - i] >= Yaxis[BottomIndex])
{
BottomIndex = Yaxis.Count - - i;
}
}
//判断是否满足条件检索条件
int newTopIndex = TopIndex;
int newBottomIndex = BottomIndex;
for (int i = ; i <= FJ; i++)
{ if (Yaxis[TopIndex + i] >= Yaxis[TopIndex])
{
newTopIndex = TopIndex + i;
}
if (Yaxis[BottomIndex - i] >= Yaxis[BottomIndex])
{
newBottomIndex = BottomIndex - i;
}
}
TopIndex = newTopIndex;
BottomIndex = newBottomIndex; //添加到结果序列
if (TopIndex <= FJ && TopIndex != )
{
index.Insert(, TopIndex);
}
if (BottomIndex >= BottomIndex - FJ && BottomIndex != Xaxis.Count - )
{
index.Add(BottomIndex);
}
return index;
}
最后,也就是最简单的峰值判断了。比一下就好了。
//根据最小峰值序列进行筛选
private static List<int> FindMinPeakValue(double minFZ, List<double> Yaxis, List<int> index)
{
List<int> finalresult = new List<int>();
for (int i = ; i < index.Count; i++)
{
if (Yaxis[index[i]] >= minFZ)
{
finalresult.Add(index[i]);
}
}
index = null;
index = finalresult;
return index;
}