一种峰值检测算法——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语言实现:
// 寻找数组最小值的下标 int argmin(int* index, int index_len) { int min_index = 0; int min = index[0]; for (int i = 1; i < index_len; i++) { if (index[i] < min) { min = index[i]; min_index = i; } } return min_index; } //寻找极值点函数 // data是存放数据的数组 //index是存放峰值点下标的数组 //len_index是峰值个数,即index数组长度 void AMPD(double* data,int* index,int *len_index) { int* p_data = (int*)malloc(sizeof(int) * size); //size可以最大为数组长度 int* arr_rowsum = (int*)malloc(sizeof(int) * size); int min_index, max_window_length; for (int i = 0; i < size; i++) { p_data[i] = 0; } for (int k = 1; k <= size / 2 + 1; k++) { int row_sum = 0; for (int i = k; i <= size - k; i++) { if ((data[i] > data[i - k]) && (data[i] > data[i + k])) row_sum -= 1; } *(arr_rowsum + k - 1) = row_sum; } /*for (int i = 0; i < size/2; i++) { printf("%d\n", arr_rowsum[i]); }*/ min_index = argmin(arr_rowsum, size/2); //此处为最大的窗口 //printf("%d\n", min_index); max_window_length = min_index; for (int k = 1; k < max_window_length + 1;k++) { for (int i = 1; i < size - k; i++) { if ((data[i] > data[i - k]) && (data[i] > data[i + k])) p_data[i] += 1; } } for (int i_find = 0; i_find < size; i_find++) { if (p_data[i_find] == max_window_length) { index[*len_index] = i_find; (*len_index) += 1; } } free(p_data); free(arr_rowsum); }
周期信号效果原始文章内太多了,这里就不展示了。
下面来展示一波我自己使用的,非周期信号的峰值寻找情况
可以根据AMPD后得到的峰值数组,进行三次样条插值进行原信号包络的获取。(参考MATLAB内envelope函数)
寻找波谷的话直接将原始数据翻转一下,就可以得到波谷的下标了。