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

概要

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

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;
    int n = static_cast<int>(data.size());
    for (int x = 2 * h_value; x < n - 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.0 * h_value);
        dst.push_back(fx);
    }
    return dst;
}

这里没有对前后2*h_value进行计算，可以使用向前差分和向后差分补上。由于我们的项目不需要，所以不再赘述。导数数组因此比原散点短，下标i对应原数组的i + 2 * h_value。后面算面积要把这个偏移加回去，不能直接拿导数下标去原数组里取

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

散点导平滑

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

求极值

由于我们求出来的是散点导数，而不是导数方程，所以直接通过f'(x) == 0判断并不合适（因为散点并不一定落在x轴上）。过零点的两侧导数异号，看相邻两个点即可：f'(x - 1) * f'(x) <= 0。

// flag > 0 只要极大值（导数由正变负），flag < 0 只要极小值（导数由负变正），0 则都要
std::vector<int> extremum(const std::vector<double>& data, int flag = 0)
{
    std::vector<int> dst;
    for (int i = 1; i < static_cast<int>(data.size()); ++i) {
        // 左右异号，则为极值
        if (data.at(i - 1) * data.at(i) <= 0) {
            if (flag > 0 && (data.at(i - 1) < 0 || data.at(i) > 0)) {
                continue;
            }
            if (flag < 0 && (data.at(i - 1) > 0 || data.at(i) < 0)) {
                continue;
            }
            // 取离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) - baseline)，baseline 取这个波峰前一个极小值的高度。积分范围本来可以取上一个极小值到下一个极小值，我们的需求要固定窗口宽度：用第一个波峰「上一个极小值到下一个极小值」的距离除以 width_ratio（实际是三分之一），再以当前波峰为中心取这么宽。

double crest_area(const std::vector<double>& data, const std::vector<int>& minima,
                  const std::vector<int>& maxima, int crest_index, int index_offset, int width_ratio = 3)
{
    if (maxima.empty() || crest_index < 0 || crest_index >= static_cast<int>(maxima.size()) || width_ratio <= 0) {
        return 0;
    }
    int first_prev = -1;
    int first_next = -1;
    for (int i = 0; i != static_cast<int>(minima.size()); ++i) {
        if (minima.at(i) < maxima.at(0)) {
            first_prev = minima.at(i);
        } else if (first_next < 0) {
            first_next = minima.at(i);
        }
    }
    if (first_prev < 0 || first_next < 0) {
        return 0;
    }
    int width = (first_next - first_prev) / width_ratio;
    if (width <= 0) {
        return 0;
    }

    int center = maxima.at(crest_index) + index_offset;
    int prev = -1;
    for (int i = 0; i != static_cast<int>(minima.size()); ++i) {
        if (minima.at(i) < maxima.at(crest_index)) {
            prev = minima.at(i);
        }
    }
    if (prev < 0) {
        return 0;
    }
    double bottom = data.at(prev + index_offset);
    int begin = center - width / 2;
    if (begin < 0) {
        begin = 0;
    }
    int end = begin + width;
    if (end > static_cast<int>(data.size())) {
        end = static_cast<int>(data.size());
    }
    double area = 0;
    for (int i = begin; i != end; ++i) {
        area += data.at(i) - 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 maxima = extremum(dv_blur_dst, 1);
    auto minima = extremum(dv_blur_dst, -1);
    int index_offset = 2 * arg.dv_h_value;

    std::vector<double> result;
    if (maxima.empty()) {
        return result;
    }
    double first_crest_area = crest_area(blur_dst, minima, maxima, 0, index_offset, arg.integral_width_ratio);
    if (first_crest_area == 0) {
        return result;
    }
    for (int i = 1; i < static_cast<int>(maxima.size()); ++i) {
        double current_area = crest_area(blur_dst, minima, maxima, i, index_offset, 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;
};