前段时间接到了个需求，数据是几千个散点，绘制出来波形图，有多个波峰，需求是求出各个波峰与第一个波峰的面积比值。乍一看需求有点难以实现，其实细化成各个模块之后也算比较简单，在这里总结一下经验方法

概要

首先分清要做什么，既然是散点，那肯定需要平滑。得出了平滑曲线之后，怎么求波峰？当然是取极大值，就是散点求导，求完导可能需要对导函数散点再次平滑。求面积就积分，至于积分取值区间，就是前一个极小值到后一个极小值的差。（我们的需求比较特殊，要取值区间固定大小，所以在这里略有不同），总结一下：

1. 散点平滑
2. 求散点导
3. 散点导平滑
4. 求极值
5. 积分求面积

下面逐条分析

散点平滑

平滑就比较简单了，这里用了均值滤波，说白了就是前后求均值。

直接上代码吧，原理简单效果显著

std::vector<double> blur(const std::vector<double>& data, int kernel_size, double center_value)
{
    std::vector<double> dst;
    for (int i = 0; i != data.size(); ++i) {
        // 均值滤波
        double sum = data.at(i) * center_value;
        double count = center_value;
        for (int j = 1; j != kernel_size; ++j) {
            if (i - j >= 0) {
                sum += data.at(i - j);
                ++count;
            }
            if (i + j < data.size()) {
                sum += data.at(i + j);
                ++count;
            }
        }
        dst.push_back(sum / count);
    }
    return dst;
}

当然也可以尝试更高级的高斯滤波或者其他方法，有兴趣可以再研究一下。

求散点导

曲线拟合比较复杂，主要波形也不固定需要分段拟合，运行环境也不一定能上ml，所有还是手工用算法操作，这里使用的是中间差分法，具体的数学原理就不是我们研究的范围了，公式可以参考一下

std::vector<double> derivative(const std::vector<double>& data, const int h_value)
{
    std::vector<double> dst;
    for (int x = 2 * h_value; x != data.size() - (2 * h_value); ++x) {
        double fx = (-data.at(x + 2 * h_value) + 8 * data.at(x + h_value) - 8 * data.at(x - h_value) + data.at(x - 2 * h_value)) / (12 * h_value);
        dst.push_back(fx);
    }
    return dst;
}

这里没有对前后2*h_value进行计算，可以使用向前差分和向后差分进行完整的散点导数计算，由于我们的项目不需要，所以不做赘述

上面是重新调整参数后的平滑曲线，下面是已经平滑过的导函数曲线

散点导平滑

没什么好说的，将导数的数据再次进行平滑即可

求极值

由于我们求出来的是散点导数，而不是导数方程，所以直接通过f'(x) == 0判断并不合适（因为散点并不一定落在x轴上），但是导数点一定是两边异号，这样也就很好判断了，通过f'(x-x0) * f'(x) <= 0即可。

std::vector<double> extremum(const std::vector<double>& data, int flag = 0)
{
    std::vector<double> dst;
    for (int i = 1; i != data.size(); ++i) {
        // 左右异号，则为极值
        if (data.at(i - 1) * data.at(i) <= 0) {
            if (flag > 0) {
                if (data.at(i - 1) < 0 || data.at(i) > 0) {
                    continue;
                }
            } else if (flag < 0) {
                if (data.at(i - 1) > 0 || data.at(i) < 0) {
                    continue;
                }
            }
            else {  // flag == 0
                ;
            }
            // 取离0近的那个
            if (std::abs(data.at(i - 1)) < std::abs(data.at(i))) {
                dst.push_back(i - 1);
            } else {
                dst.push_back(i);
            }
        }
    }
    return dst;
}

代码中的flag参数用于仅取极大值或极小值。

积分求面积

Σf(x-x0)，x0即为要选择的底部，积分范围从上一个极小值到下一个极小值即可。

我们的需求比较特殊，需要取值区间宽度固定，由于以第一个波峰作为基准，所以取其上一个极小值到下一个极小值的宽度的比值（实际取了三分之一）

double crest_area(const std::vector<double>& data, const std::vector<double>& extrem, int crest_index, int width_ratio = 3)
{
    int width = (extrem.at(1) - extrem.at(0)) / width_ratio;
    double bottom = data.at(extrem.at(1));
    // 积分求面积
    double area = 0;
    for (int i = -width; i != width; ++i) {
        area += data.at(extrem.at(i * 2)) - bottom;
    }
    return area;
}

调用

依次调各个函数吧

struct CrestRatioArg;

std::vector<double> calc_crest_ratio(const std::vector<double>& data, const CrestRatioArg & arg)
{
    auto blur_dst = blur(data, arg.kernel_size1, arg.center_value1);
    auto dv_dst = derivative(blur_dst, arg.dv_h_value);
    auto dv_blur_dst = blur(dv_dst, arg.kernel_size2, arg.center_value2);
    auto full_extrem_dst = extremum(dv_blur_dst);

    std::vector<double> result;
    if (full_extrem_dst.empty()) {
        return result;
    }
    double first_crest_area = crest_area(blur_dst, full_extrem_dst, 0, arg.integral_width_ratio);
    if (first_crest_area == 0) {
        return result;
    }
    for (int i = 1; i != full_extrem_dst.size() / 2 + 1; ++i) {
        double current_area = crest_area(blur_dst, full_extrem_dst, i, arg.integral_width_ratio);
        result.push_back(current_area / first_crest_area);
    }
    return result;
}

struct CrestRatioArg {
    int kernel_size1 = 100;
    int center_value1 = 1;
    int dv_h_value = 100;
    int kernel_size2 = 20;
    int center_value2 = 1;
    int integral_width_ratio = 3;
};