一种峰值检测算法——AMPD算法（C语言实现）

2744557306

一种峰值检测算法——AMPD算法（C语言实现）
本文算法的原始论文出处：Algorithms | Free Full-Text | An Efficient Algorithm for Automatic Peak Detection in Noisy Periodic and Quasi-Periodic Signals | HTML (mdpi.com)
有位老哥在知乎写了Python代码：python代码

在数字信号处理中，经常涉及到波峰查找算法，如振动信号分析，样条插值法求包络等。对于周期信号或者准周期信号，文章介绍了一种名为Automatic multiscale-based peak detection (AMPD)，即自动多尺度峰值查找算法。
同时对非周期信号的效果也很nice，强烈安利！
其优势是：
（1）算法本身（几乎）对信号具有良好的自适应性，唯一的假设是信号是周期的或者准周期的；
（2）抗噪能力强，后面可以看到，对周期性的要求也不是很高。
原理不多讲，可以直接来看原文，也比较简单，就是用一个多尺度的滑动窗口去两侧进行比较，寻找局部最大值。

下面为该算法的C语言实现：

1. // 寻找数组最小值的下标
   
2. int argmin(int* index, int index_len)
   
3. {
   
4.   int min_index = 0;
   
5.   int min = index[0];
   
6.   for (int i = 1; i < index_len; i++)
   
7.   {
   
8.     if (index[i] < min)
   
9.     {
   
10.       min = index[i];
    
11.       min_index = i;
    
12.     }
    
13.   }
    
14.   return min_index;
    
15. }
    
16. 
    
17. //寻找极值点函数
    
18. // data是存放数据的数组
    
19. //index是存放峰值点下标的数组
    
20. //len_index是峰值个数，即index数组长度
    
21. void AMPD(double* data,int* index,int *len_index)
    
22. {
    
23.   int* p_data = (int*)malloc(sizeof(int) * size); //size可以最大为数组长度
    
24.   int* arr_rowsum = (int*)malloc(sizeof(int) * size);
    
25.   int min_index, max_window_length;
    
26.   for (int i = 0; i < size; i++)
    
27.   {
    
28.     p_data[i] = 0;
    
29.   }
    
30.   for (int k = 1; k <= size / 2 + 1; k++)
    
31.   {
    
32.     int row_sum = 0;
    
33.     for (int i = k; i <= size - k; i++)
    
34.     {
    
35.       if ((data[i] > data[i - k]) && (data[i] > data[i + k]))
    
36.         row_sum -= 1;
    
37.     }
    
38.     *(arr_rowsum + k - 1) = row_sum;
    
39.   }
    
40.   /*for (int i = 0; i < size/2; i++)
    
41.   {
    
42.     printf("%d\n", arr_rowsum[i]);
    
43.   }*/
    
44.   min_index = argmin(arr_rowsum, size/2); //此处为最大的窗口
    
45.   //printf("%d\n", min_index);
    
46.   max_window_length = min_index;
    
47.   for (int k = 1; k < max_window_length + 1;k++)
    
48.   {
    
49.     for (int i = 1; i < size - k; i++)
    
50.     {
    
51.       if ((data[i] > data[i - k]) && (data[i] > data[i + k]))
    
52.         p_data[i] += 1;
    
53.     }
    
54.   }
    
55.   for (int i_find = 0; i_find < size; i_find++)
    
56.   {
    
57.     if (p_data[i_find] == max_window_length)
    
58.     {
    
59.       index[*len_index] = i_find;
    
60.       (*len_index) += 1;
    
61.     }
    
62.   }
    
63.   free(p_data);
    
64.   free(arr_rowsum);
    
65. }
    

复制代码

周期信号效果原始文章内太多了，这里就不展示了。
下面来展示一波我自己使用的，非周期信号的峰值寻找情况

可以根据AMPD后得到的峰值数组，进行三次样条插值进行原信号包络的获取。（参考MATLAB内envelope函数）

寻找波谷的话直接将原始数据翻转一下，就可以得到波谷的下标了。