How to detect peaks in EKG signal

Hello everyone,

i have this project with a group of people and i am in charge of writing code to detect peaks of an EKG signal.

I read the PulseSensorPlayground github and tried to understand the code for detecting the peaks. I know this is for a Pulse signal and not for an actual EKG signal.

Would it still be possible to detect peaks of an EKG signal with this? I have some EKG-data but i dont know how to adjust the code so i can test if it would actually work.

Should i just write my own code or try to make it work?

Thanks in advance,

Xenoshell

Please do not post ."zip" files, since downloading such files is a potential vector for self-activating viruses. :astonished:

IOW, competent people will simply not download them.

ok i deleted the file. The EKG data isnt even relevant for the question.

Use the forum Google search function in the upper right of this page to search for the key words of your project. You will probably find many similar previous project discussions and code to get you started.

Google "peak detection algorithm" for several general approaches.

Ok i already searched for peak detection but just found some scientific papers with math functions and no code, but i searched again.

Anyway i found something on stack overflow that looked promising. The only problem is that its written in C++ and i have kind of a hard time differentiating what i can use for Arduino and what not.

Here is the code:

#include <iostream>
#include <vector>
#include <algorithm>
#include <unordered_map>
#include <cmath>
#include <iterator>
#include <numeric>

using namespace std;

typedef long double ld;
typedef unsigned int uint;
typedef std::vector<ld>::iterator vec_iter_ld;

/**
 * Overriding the ostream operator for pretty printing vectors.
 */
template<typename T>
std::ostream &operator<<(std::ostream &os, std::vector<T> vec) {
    os << "[";
    if (vec.size() != 0) {
        std::copy(vec.begin(), vec.end() - 1, std::ostream_iterator<T>(os, " "));
        os << vec.back();
    }
    os << "]";
    return os;
}

/**
 * This class calculates mean and standard deviation of a subvector.
 * This is basically stats computation of a subvector of a window size qual to "lag".
 */
class VectorStats {
public:
    /**
     * Constructor for VectorStats class.
     *
     * @param start - This is the iterator position of the start of the window,
     * @param end   - This is the iterator position of the end of the window,
     */
    VectorStats(vec_iter_ld start, vec_iter_ld end) {
        this->start = start;
        this->end = end;
        this->compute();
    }

    /**
     * This method calculates the mean and standard deviation using STL function.
     * This is the Two-Pass implementation of the Mean & Variance calculation.
     */
    void compute() {
        ld sum = std::accumulate(start, end, 0.0);
        uint slice_size = std::distance(start, end);
        ld mean = sum / slice_size;
        std::vector<ld> diff(slice_size);
        std::transform(start, end, diff.begin(), [mean](ld x) { return x - mean; });
        ld sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0);
        ld std_dev = std::sqrt(sq_sum / slice_size);

        this->m1 = mean;
        this->m2 = std_dev;
    }

    ld mean() {
        return m1;
    }

    ld standard_deviation() {
        return m2;
    }

private:
    vec_iter_ld start;
    vec_iter_ld end;
    ld m1;
    ld m2;
};

/**
 * This is the implementation of the Smoothed Z-Score Algorithm.
 * This is direction translation of https://stackoverflow.com/a/22640362/1461896.
 *
 * @param input - input signal
 * @param lag - the lag of the moving window
 * @param threshold - the z-score at which the algorithm signals
 * @param influence - the influence (between 0 and 1) of new signals on the mean and standard deviation
 * @return a hashmap containing the filtered signal and corresponding mean and standard deviation.
 */
unordered_map<string, vector<ld>> z_score_thresholding(vector<ld> input, int lag, ld threshold, ld influence) {
    unordered_map<string, vector<ld>> output;

    uint n = (uint) input.size();
    vector<ld> signals(input.size());
    vector<ld> filtered_input(input.begin(), input.end());
    vector<ld> filtered_mean(input.size());
    vector<ld> filtered_stddev(input.size());

    VectorStats lag_subvector_stats(input.begin(), input.begin() + lag);
    filtered_mean[lag - 1] = lag_subvector_stats.mean();
    filtered_stddev[lag - 1] = lag_subvector_stats.standard_deviation();

    for (int i = lag; i < n; i++) {
        if (abs(input[i] - filtered_mean[i - 1]) > threshold * filtered_stddev[i - 1]) {
            signals[i] = (input[i] > filtered_mean[i - 1]) ? 1.0 : -1.0;
            filtered_input[i] = influence * input[i] + (1 - influence) * filtered_input[i - 1];
        } else {
            signals[i] = 0.0;
            filtered_input[i] = input[i];
        }
        VectorStats lag_subvector_stats(filtered_input.begin() + (i - lag), filtered_input.begin() + i);
        filtered_mean[i] = lag_subvector_stats.mean();
        filtered_stddev[i] = lag_subvector_stats.standard_deviation();
    }

    output["signals"] = signals;
    output["filtered_mean"] = filtered_mean;
    output["filtered_stddev"] = filtered_stddev;

    return output;
};

int main() {
    vector<ld> input = {1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0,
                        1.0, 1.0, 1.0, 1.1, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 1.1, 1.0, 0.8, 0.9, 1.0,
                        1.2, 0.9, 1.0, 1.0, 1.1, 1.2, 1.0, 1.5, 1.0, 3.0, 2.0, 5.0, 3.0, 2.0, 1.0, 1.0, 1.0, 0.9, 1.0,
                        1.0, 3.0, 2.6, 4.0, 3.0, 3.2, 2.0, 1.0, 1.0, 0.8, 4.0, 4.0, 2.0, 2.5, 1.0, 1.0, 1.0};

    int lag = 30;
    ld threshold = 5.0;
    ld influence = 0.0;
    unordered_map<string, vector<ld>> output = z_score_thresholding(input, lag, threshold, influence);
    cout << output["signals"] << endl;
}

In the end i would like to have a function which i can call in the loop() to give a new AnalogRead value

The Arduino is programmed in C++

In the program you posted, main() does what setup() and loop() do in an Arduino program, by calling other functions (most are in included libraries).

I just googled "peak detection C code", which turned up many pages with example code. Here is one of them: Peak detection of a time series « Stack Overflow

Yeah true Arduino code is written in C++. Its still much simpler.

I just dont understand the code that well.

Even the last link jremington send me, i dont really understand what is going on with this code. It seems to detect only one peak for a whole array. I need to detect many peaks for a ECG. Then time those peaks to calculate the BPM.

OK.

stored_reading = -100000000000000000000000000000000000000.00

My first reading could be a +.01, so if first reading is greater then stored_reading then stored readings becomes first reading.

If 2nd reading is greater then stored_reading then 2nd reading become stored reading.

If 3rd reading is greater then stored_reading then 3rd reading become stored reading.

If 4th reading is greater then stored_reading then 4th reading become stored reading.

A pattern begins to emerge, no?

We strongly recommend that beginners start with the examples that come with the Arduino development software, in order to learn the programming language and the special features of the Arduino. If you skip that step, expect endless frustration.

That said, you could do a bit more research. Googling "arduino heart rate" turns up several working projects, so the problem has been solved many times.

@Idahowalker
Yeah i get that part, it makes sense to always store the highest value to get the peak.

What i struggle with is, if you look at the peak detection algorithm i posted earlier, this part:

 if (abs(input[i] - filtered_mean[i - 1]) > threshold * filtered_stddev[i - 1]) {
            signals[i] = (input[i] > filtered_mean[i - 1]) ? 1.0 : -1.0;
            filtered_input[i] = influence * input[i] + (1 - influence) * filtered_input[i - 1];
        } else {
            signals[i] = 0.0;
            filtered_input[i] = input[i];
        }
        VectorStats lag_subvector_stats(filtered_input.begin() + (i - lag), filtered_input.begin() + i);
        filtered_mean[i] = lag_subvector_stats.mean();
        filtered_stddev[i] = lag_subvector_stats.standard_deviation();

Ok maybe this is just the calculations how to get the peak, i just need to change the code a bit. For every peak i need a timestamp (we have a timer class) to calculate the BPM.
This code returns a vector, but of what? Is it really just a 1 or -1? Wouldnt i then just need to write something like this:

signals[i] = (input[i] > filtered_mean[i - 1]) ? 1.0 : -1.0;
timestamps[i] = timer.current_time();

@jremington
i googled something similar. I got to a page with a lot of Pulse detection projects. Im just not sure if Pulse Peak detection works as well for ECG signals since its not the same curve.

You sent me a link to an algorithm too. It seems this one just records the highest value given, which is obviously fine. For a ECG signal it could be that one peak, that needs to be detected, is lower than another peak.
I cant see if this would detect every peak in a series (i attached a ECG curve to better explain what i mean)

This line just sets the array element signals[ i] to be 1 or -1, depending on whether the current value of input is greater than the most recent filtered value. It does not "detect a peak".

signals[i] = (input[i] > filtered_mean[i - 1]) ? 1.0 : -1.0;

You don't need to detect peaks in order to measure heart rate. You need only measure the time difference between equivalent points in successive heart beat cycles, which can take advantage of any feature of the signal.

The signal is very clean and it would be trivial to extract heart rate from it. One idea is to measure the time period between two points in successive cycles, where the negative excursion, from one sample to the next, goes below a certain value, say 200.

OP's image:

As an aside, the graph shows all peaks going to zero (or lower). That suggests you are exposing the analog input to negative voltages, which can destroy it.

@jremington,

i just wrote this and it works pretty well,

how i call it in the loop:

 if (x < 980) {
            BPM = EKG_calculate_1.detect_peaks(ekgvalues[x]);
            delay(100);
            Serial.println("ekgvalues:");
            Serial.println(ekgvalues[x]);
            Serial.println("BPM:");
            Serial.println(BPM);
          }
          x++;

There is a boolean value to check if a peak was detected and a boolean first_peak value to see if the first peak was already detected.
Can you maybe give some feedback what could be improved?

How would i implement the variable middle instead of using 1500? I dont want to create an overflow with all_values but i probably would have to reset it at some point. How would i do that?

unsigned short EKG_calculate::detect_peaks(unsigned short value) {

  //counter++;
  //all_values += value;
  //middle = all_values / counter;

  if (time_on == false) {
    timer1.start();
    time_on = true;
  }

  if (value > (1500 - THRESHOLD)) {
    peak = false;
  }

  if (value < (1500 - THRESHOLD) && (peak == false)) { // Detection of Peak

    if (peak_one == true) { //This detects if the first peak was already detected
      peak_time_2 = timer1.elapsedMilliseconds();
      peak_one = false;
      peak = true;
      Serial.println("Peak time 2:");
      Serial.println(peak_time_2);

      BPM = (peak_time_1 - peak_time_2) * 60;

      return BPM;
    }

    if (peak_one == false) { //no first peak
      peak_time_1 = timer1.elapsedMilliseconds(); //first peak time
      peak_one = true;
      peak = true;
      Serial.println("Peak time 1:");
      Serial.println(peak_time_1);
    }
  }

  return BPM;
}

Thanks for your help btw :slight_smile:

How would i implement the variable middle instead of using 1500?

Use a running average, like the code posted in reply #10.

I dont want to create an overflow with all_values but i probably would have to reset it at some point.

You aren't using all_values right now, but the technique of using a running average takes care of the problem.

Note that there is a stupid bug in the code linked above. It won't work correctly until "numReadings" points have been collected.

What do you mean with "until "numReadings" points have been collected"?

Look at the Smoothing example code, and decide for yourself what happens after the first data point has been collected. Or run the example and study the output.