C# 实现寻峰算法的简单优化（包含边峰，最小峰值，峰距）
核心寻峰算法的原理参考Ronny，链接：投影曲线的波峰查找，
C#翻译原理代码参考sowhat4999，链接：C#翻译Matlab中findpeaks方法
前人种树，后人乘凉。感谢原作者详细的解释说明。
这里先把翻译代码贴一下(略微的修改了sowhat4999代码中的几个参数）
//调用方法 List<double> data = new List<double>{25, 8, 15, 5, 6, 10, 10, 3, 1, 20, 7}; List<int> index = getPeaksIndex(trendSign(oneDiff(data)));
//第一次寻峰（基本峰距为1）算法 private double[] oneDiff(List<double> data) { double[] result = new double[data.Count - 1]; for (int i = 0; i < result.Length; i++) { result[i] = data[i + 1] - data[i]; } return result; } private int[] trendSign(double[] data) { int[] sign = new int[data.Length]; for (int i = 0; i < sign.Length; i++) { if (data[i] > 0) sign[i] = 1; else if (data[i] == 0) sign[i] = 0; else sign[i] = -1; } for (int i = sign.Length - 1; i >= 0; i--) { if (sign[i] == 0 && i == sign.Length - 1) { sign[i] = 1; } else if (sign[i] == 0) { if (sign[i + 1] >= 0) { sign[i] = 1; } else { sign[i] = -1; } } } return sign; } private List<int> getPeaksIndex(int[] diff) { List<int> data = new List<int>(); for (int i = 0; i != diff.Length - 1; i++) { if (diff[i + 1] - diff[i] == -2) { data.Add(i + 1); } } 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> { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 }; List<double> Yaxis = new List<double> { 25, 8, 15, 5, 6, 10, 10, 3, 1, 20, 7 }; // 峰距 int DisPeak = 3; // 峰距为3时得到的脚标 List<int> index = getPeaksIndex(trendSign(oneDiff(Yaxis))); // 已进行的判断 int level = 1; // 扩大峰距范围范围算法 while (DisPeak > level) { level++; List<int> result = DoPeakInstance(Yaxis, index, level); index = null; index = result; } // 获取两侧满足条件的边峰序列 index = GetBothSidePeakIndex(Xaxis, Yaxis, 1, 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 = 0; i < index.Count; i++) { //判断是否超出下界和上界 if (index[i] - level >= 0 && 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 = 0; int BottomIndex = Yaxis.Count - 1; for (int i = 0; i < FJ; i++) { if (Yaxis[i] >= Yaxis[TopIndex]) { TopIndex = i; } if (Yaxis[Yaxis.Count - 1 - i] >= Yaxis[BottomIndex]) { BottomIndex = Yaxis.Count - 1 - i; } } //判断是否满足条件检索条件 int newTopIndex = TopIndex; int newBottomIndex = BottomIndex; for (int i = 0; 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 != 0) { index.Insert(0, TopIndex); } if (BottomIndex >= BottomIndex - FJ && BottomIndex != Xaxis.Count - 1) { 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 = 0; i < index.Count; i++) { if (Yaxis[index[i]] >= minFZ) { finalresult.Add(index[i]); } } index = null; index = finalresult; return index; }
