前言

最近研发部经理让研究一下寻峰算法，因为我们的软件在处理数据之前，单个设备每秒会从硬件读卡器读取上万条电弧放电的电流原始数据。这些原始数据里面有比较多的脏数据，需要基于某种规则把脏数据过滤掉，再对有效数据进行解析、存储和展示。由于我们的软件主要是用来监控变电站的设备是否存在电流过载导致跳闸停电的风险，所以针对电流原始数据，我们只需要关注峰值情况，对于波谷等数据可以直接过滤掉，同时由于存在毛刺现象导致误告警，我们也需要过滤掉毛刺，所以寻峰算法刚好可以解决我们遇到的问题。

算法原理

首先这些数值是一系列离散的点，把点连成线会得到一个上下波动的曲线，一个简单的想法就是拿当前的值跟前后两个值进行大小比较：
$F(x) > F(x - 1) && F(x) > F(x + 1)$
满足上述条件即为波峰。
但是上述算法有比较明显的问题，当曲线变平时，峰值点反而被过滤了。加上>= F(x + 1)是否就OK了呢，对于阶段状上升的曲线也会有问题。

经过调研，发现可以通过求导的方式寻找峰值点，对于峰值点来说，满足一阶导数为0 && 二阶导数为负。核心步骤：

  1. 计算一阶导数：首先，需要计算给定曲线的一阶导数。一阶导数反映了曲线在某点的斜率，即切线的倾斜程度。

  2. 符号化一阶导数：接下来，对计算得到的一阶导数进行符号化处理。这步操作的目的在于确定切线的方向，而不是关注斜率的具体数值大小。

  3. 确定一阶导数为零的点：在符号化的一阶导数中寻找值为零的点。这些点表示原曲线上的切线水平，可能是潜在的波峰或波谷点。

  4. 调整一阶导数符号：对于一阶导数为零的点，进一步分析它们周围的点，确定它们所在的坡面梯度方向（上坡或下坡）。根据坡面方向，调整这些点的一阶导数符号，以确保它们与所在坡面的梯度方向一致。

  5. 计算二阶导数：计算上述处理过的一阶导数的二阶导数。二阶导数反映了切线斜率的变化情况，即曲线的弯曲程度。

  6. 寻找二阶导数非零的点：在二阶导数中寻找值为正或负（通常是2或-2）的点。这些点表示原曲线在该处由凸变凹或由凹变凸，即波峰或波谷的转折点。

  7. 确定波峰和波谷：综合一阶导数和二阶导数的信息，确定原曲线中的波峰和波谷点。波峰点应满足一阶导数为零且二阶导数为负，而波谷点应满足一阶导数为零且二阶导数为正。


算法实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#include <stdio.h>  
#include <stdlib.h>  
  
void findPeaks(const int* v, int size, int** peakPositions, int* peakCount) {  
    // 为一阶差分数组分配内存  
    int* diff_v = (int*)malloc((size - 1) * sizeof(int));  
    if (!diff_v) {  
        perror("Failed to allocate memory for diff_v");  
        return;  
    }  
      
    // 计算V的一阶差分  
    for (int i = 0; i < size - 1; i++) {  
        if (v[i + 1] - v[i] > 0)  
            diff_v[i] = 1;  
        else if (v[i + 1] - v[i] < 0)  
            diff_v[i] = -1;  
        else  
            diff_v[i] = 0;  
    }  
      
    // 对diff_v进行处理，设置一阶导数为0的点的方向  
    for (int i = size - 2; i >= 0; i--) {  
        if (diff_v[i] == 0 && i == size - 2) {  
            diff_v[i] = 1;  
        } else if (diff_v[i] == 0) {  
            if (i + 1 < size - 1 && diff_v[i + 1] >= 0)  
                diff_v[i] = 1;  
            else  
                diff_v[i] = -1;  
        }  
    }  
      
    // 寻找二阶导数非零的点，即波峰或波谷  
    int* tempPeakPositions = NULL;  
    int tempPeakCount = 0;  
    for (int i = 0; i < size - 2; i++) {  
        if (diff_v[i + 1] - diff_v[i] == -2) {  
            // 动态扩展peakPositions数组  
            tempPeakPositions = (int*)realloc(tempPeakPositions, (tempPeakCount + 1) * sizeof(int));  
            if (!tempPeakPositions) {  
                perror("Failed to allocate memory for tempPeakPositions");  
                free(diff_v);  
                return;  
            }  
            tempPeakPositions[tempPeakCount++] = i + 1;  
        }  
    }  
      
    // 分配最终peakPositions数组内存并复制结果  
    *peakPositions = (int*)malloc(tempPeakCount * sizeof(int));  
    if (!*peakPositions) {  
        perror("Failed to allocate memory for peakPositions");  
        free(diff_v);  
        free(tempPeakPositions);  
        return;  
    }  
    for (int i = 0; i < tempPeakCount; i++) {  
        (*peakPositions)[i] = tempPeakPositions[i];  
    }  
      
    // 设置返回的波峰数量  
    *peakCount = tempPeakCount;  
      
    // 释放临时数组内存  
    free(tempPeakPositions);  
    free(diff_v);  
}  
  
int main() {  
    // 示例数组  
    int v[] = {1, 3, 7, 1, 0, -2, -4};  
    int size = sizeof(v) / sizeof(v[0]);  
      
    // 用于存储波峰位置的数组和波峰数量  
    int* peakPositions = NULL;  
    int peakCount = 0;  
      
    // 调用findPeaks函数  
    findPeaks(v, size, &peakPositions, &peakCount);  
      
    // 打印波峰位置  
    printf("Peak positions: ");  
    for (int i = 0; i < peakCount; i++) {  
        printf("%d ", peakPositions[i]);  
    }  
    printf("\n");  
      
    // 释放波峰位置数组内存  
    free(peakPositions);  
      
    return 0;  
}