c++ - 为噪声信号寻找实时可靠和精确的峰值检测算法

Question

我现在正在做一个项目，基本上我需要通过 RPi Pico 精确实时测量霍尔传感器测量的峰值，通过 Arduino-Pico 库在 Arduino IDE 中编码，问题是，信号非常嘈杂而且不是每个峰是完美的，很多都被摧毁了，我需要有可靠和精确的算法。如果有人处理过类似的问题并且能够给我一些建议，我将不胜感激。信号看起来像这样:

这是来自霍尔传感器的原始信号:

这是 4 个先前值的平均信号(数据与前一个不同):

我试过两种方法:一种是设置一个highThreshold，当这个值超过它时，程序开始寻找当前区域的最高数字；这行得通，尽管在数据有些损坏且图形没有正确峰值(输入 curVal)的部分无效。

  HighThresCoeff = 0.85
  //code for highThreshold generation
  vals[i]=curVal;
  i++;
  if(i==arrSize){
    low=getLow(vals);
    high=getHigh(vals);
    highThreshold=((high-low)*HighThresCoeff+low);
    i=0;
  }
  //peak detection
  if (curVal > highThreshold) {
    activated = true;
    if(curVal > lastHigh){
      lastHigh = curVal;
      lastHighTime = micros();
    }
  } else if (activated == true) {
    lastHigh = 0;
    activated = false;
    t2 = t1;
    t1 = lastHighTime;
    // code for processing the time of the peak
  }

我尝试的另一种方法也是基于highThreshold，虽然我是在找时间，当图形值超过和低于阈值时，再做一个平均值；虽然这样更好，但由于噪音，我仍然没有得到我希望的那么好的数据。

  HighThresCoeff = 0.85
  //code for highThreshold generation
  vals[i]=curVal;
  i++;
  if(i==arrSize){
    low=getLow(vals);
    high=getHigh(vals);
    highThreshold=((high-low)*HighThresCoeff+low);
    i=0;
  }
  //peak detection
  if (curVal > highThreshold) {
    tss = micros();
    activated = true;
  } else if (activated == true) {
    activated = false;
    tse = micros();
    t2 = t1;
    t1 = tss + ((tse - tss) / 2);
    //code for processing the time further
  }

附加信息:

• 延迟时间:如果延迟低于峰峰值时间的 1/3，并且是可预测的或恒定的，没关系。
• 示例数据: https://github.com/Atores1/exampleData

Answer 1

注意到 OP 提供了指向原始 int 数据的链接，我通过移动平均过滤器运行它。移动平均滤波器的优点是不需要将缓冲区中的所有样本相加，而只需减去掉落的样本并将新样本添加到缓冲区内容的初始总和中。更少的计算工作和内存访问。

这是用原始数据覆盖的过滤结果:

下面是读取原始数据以及输出同步原始数据和过滤数据的代码。

#include <iostream>
#include <fstream>
#include <vector>
#include <array>
#include <numeric>
#include <algorithm>
#include <cstdio>
#include <type_traits>

using std::array, std::vector, std::size_t;

using sample_type = int;            // data sample type, either int or double
constexpr int Global_Filter_N = 41; // filter length, must be odd

// moving average filter
template <typename T=sample_type, int N=Global_Filter_N>
class Filter_MA
{
public:
    T clk(T in)
    {
        sum += in - buf[index];
        buf[index] = in;
        index = (index + 1) % N;
        if constexpr (std::is_floating_point_v<T>)
            return sum / N;
        else
            return (sum + (N / 2)) / N;
    }
    bool update_vectors(const vector<T>& vin, vector<T>* pvout, vector<T>* prawout = nullptr)
    {
        if (vin.size() <= N || pvout == nullptr)
            return false;
        pvout->reserve(vin.size() - N);
        if (prawout != nullptr)
            pvout->reserve(vin.size() - N);
        for (size_t i = 0; i < N; i++)
            clk(vin[i]);
        for (size_t i = N; i < vin.size(); i++)
        {
            pvout->push_back(clk(vin[i]));
            if (prawout != nullptr)
                prawout->push_back(vin[i - N / 2]);
        }
        return true;
    }
private:
    array<T, N> buf{};  // moving average buffer
    T sum{};            // running sum of buffer
    size_t index{};     // current loc remove output, add input 
};

template <typename T=sample_type>
std::pair<T, T> peak_detect(T y1, T y2, T y3)
{
    // scale pk location by 100 to work with int arith
    T pk = 100* (y1 - y3) / (2 * (y1 - 2 * y2 + y3));
    T mag =  2 * y2 - y1 - y3;
    return std::pair{ pk, mag };
}

struct WaveInfo {
    sample_type w_mean{};
    sample_type w_max{};
    sample_type w_min{};
    vector<sample_type> peaks;
    vector<sample_type> mags;
};

inline WaveInfo get_wave_info(std::vector<sample_type> v)
{
    constexpr int N = Global_Filter_N;
    static_assert(Global_Filter_N & 1, "filter must be odd number");
    WaveInfo w;
    w.w_max = *std::max_element(v.begin(), v.end());
    w.w_min = *std::min_element(v.begin(), v.end());
    // "0ll + sample_type{}" Produces either a double or long long int depending on sample_type to stop overflow if > 2M samples
    w.w_mean = static_cast<sample_type>(std::accumulate(v.begin(), v.end(), 0ll + sample_type{}) / std::size(v));
    sample_type pos_thresh = w.w_mean + (w.w_max - w.w_mean) / 10;  // 10% above ave.
    sample_type neg_thresh = w.w_mean + (w.w_min - w.w_mean) / 10;  // 10% below ave
    int search_polarity = 0;    // if 0 prior peak polarity not determined
    for (int i = 0; i < int(v.size()) - N; i++)
    {
        const int center = N/2;
        if (v[i] > pos_thresh && v[i] > v[i + N - 1] && v[i] < v[i + center] && search_polarity >= 0)
        {
            search_polarity = -1;
            auto results = peak_detect(v[i], v[i + center], v[i + N - 1]);
            w.peaks.push_back(results.first * center / 100 + i + center);
            w.mags.push_back(results.second);
        }
        if (v[i] < neg_thresh && v[i] < v[i + N - 1] && v[i] > v[i + center] && search_polarity <= 0)
        {
            search_polarity = 1;
            auto results = peak_detect(v[i], v[i + N / 2], v[i + N - 1]);
            w.peaks.push_back(results.first * center / 100 + i + center);
            w.mags.push_back(-results.second);
        }
    }
    return w;
}

// Used to get text file int samples
vector<sample_type> get_raw_data()
{
    std::ifstream in("raw_data.txt");
    vector<sample_type> v;
    int x;
    while(in >> x)
        v.push_back(x);
    return v;
}

int main()
{
    Filter_MA filter;
    vector<sample_type> vn = get_raw_data();
    vector<sample_type> vfiltered;
    vector<sample_type> vraw;
    if (!filter.update_vectors(vn, &vfiltered, &vraw))
        return 1;   // exit if update failed

    // file with aligned raw and filtered data
    std::ofstream out("waves.txt");
    for (size_t i = 0; i < vfiltered.size(); i++)
        out << vraw[i] << " " << vfiltered[i] << '\n';

    // get filtered file metrics
    WaveInfo info = get_wave_info(vfiltered);
    out.close();

    //  file with peak locs and magnitudes
    out.open("peaks.txt");
    for (size_t i = 0; i < info.peaks.size(); i++)
        out << info.peaks[i] << " " << info.mags[i] << '\n';
}

这是前 4 个峰值的峰值信息输出。第一列是位置，第二列是峰值的相对大小，

116 43
344 32
577 44
812 37