一种用于噪声周期性和准周期性信号的自动峰值检测高效算法

Felix Scholkmann *, Jens Boss 和 Martin Wolf
生物医学光学研究实验室，苏黎世大学医院新生儿科，瑞士苏黎世8091；电子邮件：bossj@student.ethz.ch (J.B.); martin.wolf@usz.ch (M.W.)

• 通信作者；电子邮件：felix.scholkmann@usz.ch; 电话：+41-(0)44-255-93-26。
  收到：2012年8月3日；修订：2012年10月29日 / 接受：2012年11月13日 / 发布：2012年11月21日

摘要：我们提出了一种新的方法用于自动检测噪声周期性和准周期性信号中的峰值。这种新方法被称为基于自动多尺度的峰值检测（AMPD），其基础是在局部极大值小波图的情境下进行计算与分析。该小波图是一个包含局部极大值的尺度依赖性矩阵。我们通过对模拟和现实信号应用AMPD算法来展示所提出方法的有效性。

关键词：峰值检测；局部极大值小波图；多尺度局部极大值检测；自动多尺度基础的峰值检测（AMPD）算法

1. 引言
   在信号处理中，检测信号中的峰值是一个关键步骤。迄今为止，已经开发了许多不同的方法，包括基于传统窗口阈值技术的 [1–4]，小波变换 [5–13]，希尔伯特变换 [14]，结合希尔伯特和小波变换 [15]，人工神经网络 [16,17]，使用模板的技术 [18,19]，形态过滤 [20–22]，非线性过滤 [23–26]，卡尔曼过滤 [27]，Gabor过滤 [28]，高斯二阶导数过滤 [29]，线性预测分析 [30]，高阶统计 [31]，K-均值聚类 [32]，模糊C-均值聚类 [33]，经验模态分解 [34]，隐马尔可夫模型 [35]，以及利用熵 [36]，动量 [37]，直方图/累积分布函数 [38]，强度加权方差 [39]，随机共鸣 [40]，或平滑非线性能量算子 [41]的技术。

现有的多数峰值检测算法的问题在于，算法越通用，适用性越广，则需要选择的自由参数（即窗口长度或阈值）就越多。而相对而言，参数较少的算法则通常只能应用于特定的情况，例如在心电图（ECG）信号中检测R峰 [1,7,8,10,12,14–17,20–24,30–34] 或在色谱数据中检测峰值 [3,4]。此外，分析信号中的噪声的存在也给许多峰值检测算法带来了挑战。

现有峰值检测方法的不令人满意的现状促使我们开发了一种简单而有效的新峰值检测算法，用于噪声周期性和准周期性信号。该算法的基础在于认识到，可以通过 (i) 使用多尺度技术检测信号的所有局部极大值，并 (ii) 通过自动分析应用多尺度技术的结果来找到“真实峰值”，而保持算法的通用适用性。具体而言，我们的工作的目标是开发一个具有以下特性的峰值检测算法：(i) 没有需要在分析前由用户选择的自由参数值；(ii) 能够检测周期性和准周期性信号中的峰值；(iii) 噪声对高低频峰值检测效率的影响相当稳健。

需要指出的是，本工作的目标是专门开发一个峰值检测算法，用于检测噪声周期性和准周期性信号中的峰值。因此，本研究并不打算开发一个适用于非周期性信号的通用峰值检测方法，例如在分析化学领域中通常获得的信号，包括质谱、光谱或气相色谱数据。

以上是论文的摘要以及现阶段研究状况。

下面是算法部分，我总结下：就是输入最大窗口（最大峰值区间），然后会根据波形数据计算出实际最大窗口，最后把符合这种峰值区间筛选出来，每个区间的最高点就是峰值点。

论文中的具体内容如下：

下面是结果图：

我根据以上论文提到的算法写了个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, std::vector<int>& index,int WindowsSize,int size)
{        
        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 / WindowsSize); //此处为最大的窗口
        //printf("%d\n", min_index);
        max_window_length = min_index;
        for (int k = 1; k < max_window_length + 1; k++)
        {
                for (int i = k; 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.push_back(i_find);
                }
        }
        free(p_data);
        free(arr_rowsum);
}