Calculating BPM - from PPG signal

Hi everyone,

I've been able to create a script that reads heart rate - displaying the value in real time on a graph.

Now i need to calculate the number of peaks in the graph and work out the BPM

Any simple script solution?

Any help will be most appreciated! Thanks!

[Ahmet YALCIN]

please share if you got a solution to this problem. thx

[willing2learn]

Hey, use try this code out:

 ydpos = diff(ydata>1000)>0;     % these data points are shown in the
                                 % attached figure below
 dxmax = max(xdata(ydpos==1));
 dxmin = min(xdata(ydpos==1));
 dx    = (dxmax - dxmin)/(sum(ydpos)-1)
 dx =

written by Mischa Kim

[Ahmet YALCIN]

Thank you Mischa,

I will try your suggestion… I developed the following algorithm (similar to yours) to find the peaks, and it seems also to work:

- Smooth by using 1D Gaussian filter and Convolution g = gausswin (5); g = g / sum(g); and by convolution with 'conv' --> conv_sig = conv(sig, g, 'same');

- Calculate and plot the 1st-derivative of the convolved signal

- Find the peaks of 1st-derivative signal [pks_1stderiv, loc_1stderiv_peaks]= findpeaks(conv_sig(d), 'MINPEAKHEIGHT',1,'MINPEAKDISTANCE' ,18)

[Federation]

Hi everyone, I need to calculate Heat Rate variability (HRV) by using Frequency domain measurement and find very low frequency(VL), low frequency(LF) , high frequency(HF) and total power. For example If X is RR interval X=[0.456, 0.567, 0.457, ..........,0.770] where length X = 600 and in second. How can I find (VL for X from 0.0-0.04), (LF for X from 0.04 - 0.15),(HF for X from 0.15 - 0.40) and total power in millisecond square.

[sarang bagade]

please tell me how you displayed the value in real time on a graph

[Supriya V]

Sir/Madam, could you share source code for ppg extraction from face video that is real time and along with I want to display a graph for that.

Answer 1

[Mischa Kim]

Check out the findpeaks command. Alternatively (if you do not have access to the Signal Processing Toolbox) you could do a Fourier Transform and look at the dominant frequency components.

[willing2learn]

Hey, thanks for the fast response!

I'm not great with matlab

Using the find peaks, how am i to calculate the frequency?

Using the fft function, i have got the attached figure - HR is less than 10Hz which is reasonable.

But i need to detect and display actual pulse count - e.g '70 beats over a 60second sample'

Thanks

[Mischa Kim]

I think you forgot to attach. If you have frequency in Hz (beats per sec) then you just need to multiply by 60 to get beats per minute. Does that help?

[willing2learn]

Hi, yes that helps

my problem is that i want to collect the data (several figures) and have a script that analyses the peaks for each figure. Then give an actual value. With my fft graph how will i do it? Or with the peakfind?

Thanks

[Mischa Kim]

Looking at the first figure, what resembles one beat? Signal values greater than 1000, 500? Does the peak lasting from x = 56 to x = 71 count as one single peak?

[willing2learn]

Hi, yes signal values > 1000 count as one peak - and x=56 - x=71 is a single pulse / peak.

The square wave is due to the saturation of my current circuit. When i turn down the gain the graph should look close to an ecg curve.

Thanks

[Mischa Kim]

OK. In this case you could simply use

 ydpos = diff(ydata>1000)>0;     % these data points are shown in the
                                 % attached figure below
 dxmax = max(xdata(ydpos==1));
 dxmin = min(xdata(ydpos==1));
 dx    = (dxmax - dxmin)/(sum(ydpos)-1)
 dx =
   37.500000000000000

xdata and ydata are the data from your HR1.fig, dx is the average duration between individual beats. I don't know what your units are on the x-axis, so you need to figure out how to convert dx into beats per minute.

Does this get you what you need and answer your question?

[willing2learn]

Cool! This is almost what i want, just that i don't know how convert the x-axis into time, it is currently number of samples i think - over a 10second period. Below is my code so far:

    % HR_real_time
  clear all;
  a = arduino('COM24')
  tic
  toc

time_start = tic; %start taking samples
time_stop = toc; %stop taking samples
sample_time = 60; %1k samples
xstop = 1000;
x = [0,0];

while(toc < sample_time)

      b = a.analogRead(0); %storing the analog values in b
      x = [x,b]; %x as a matrix?, updating x - causing the shift
      figure (1)
      plot(x);
      axis([0,xstop,0,1500]);
      grid
      %t = tic; % increment t by step size
      drawnow;    
      %pks = findpeaks(x)

end

      xfft = fft(x);
      figure (2)
      freq=0:(1000)/length(x):(1000)/2;
      plot(freq,abs(xfft(1:length(x)/2+1)))

      % code
    end

Thanks

[willing2learn]

Hi, i'm still unable to change the x-axis to time and use the code you provided. Could you further assist please?

Thanks

[willing2learn]

Hi Mischa, the time axis has been calculated and your solution added to the code. I am very close! The only issue is there seems to be a loss of counts somewhere - e.g when my heart rate is around 63BPM, Matlab returns around 55...

Could you please look at the updated code to see if i have missed something or misunderstood?

Thanks a lot

% HR_real_time
clear all;

a = arduino('COM24')
tic
toc

time_start = tic; %start taking samples
time_stop = toc; %stop taking samples
sample_time = 20; %1k samples
xstop = 20;
x = [];
time_ind=[];%zeros(600,1);
i=1;
BPM = 0;
while(toc < sample_time)

      b = a.analogRead(0); %storing the analog values in b
      time_ind=[time_ind, toc];%New array that stores time index
      i=i+1; %everytime loop runs, increment
      x = [x,b]; %x as a matrix?, updating x - causing the shift
      figure (1)
      plot(time_ind,x);%Plot the time index
      axis([0,xstop,0,1500]);
      ylabel('Magnitude');
      xlabel('Time');
      title('Heart Beat in Time Domain');
      grid
      drawnow;    

end

Plots with finished loop

      xfft = fft(x);
      figure (2)
      freq=0:(1000)/length(x):(1000)/2;
      plot(freq,abs(xfft(1:length(x)/2+1)));
      ylabel('Magnitude');
      xlabel('Frequency');
      title('Heart Beat in Frequency Domain');

      figure(3)
      plot(time_ind);
      ylabel('Magnitude');
      xlabel('Samples');
      title('Heart Beat');

      %%Calculating Frequency

      ydpos = diff(x>1000)>0;     
      %at pulse instant, current value is larger than previous therfore diff>0
      %at every other point, diff<0 or diff=0;
      dxmax = max(time_ind(ydpos==1));% dxmax assigned max time value when ydpos>0 
      dxmin = min(time_ind(ydpos==1));% dxmin assigned min time value when ydpos>0 
      dx    = (dxmax - dxmin)/(sum(ydpos)-1); % calculate frequency in 1 second???
      BPM   =  (dx*60) % Get BPM in a minute

[Mischa Kim]

Hard to say. In HR1 update.fig I am counting 21 beats in less than 19.5 sec

 21/(19.89-0.57)*60
 ans =
   65.217391304347828

Is this about what you'd expect?

[willing2learn]

Hi, yes that's what i expected around 63 BPM, but i ended up with around 55. I just ran another test getting 25 pulses, but the code only returns 48 BPM - should be 75 (25 beats in 20 seconds).

There must be a small mistake somewhere

[willing2learn]

Hi Mischa thanks a lot for your help so far.

I have increased the sampling period to 60 seconds, ran several tests and the code is 100% accurate to within 10 counts (when the pulses are manually counted).

Most times it is accurate to with 5 beats.

Any idea how to fine tune it a bit more?

Thanks

[Supriya V]

Sir/Madam, could you share source code for ppg extraction from face video that is real time and along with I want to display a graph for that.

Answer 2

[T. Thinh Nguyen]

Check out the code I write in here.

https://uk.mathworks.com/matlabcentral/fileexchange/53364-heart-rate--spo2-using-ppg

The heart rate estimation is done both by FFT or Peak detection.

[Supriya V]

Sir/Madam, could you share source code for ppg extraction from face video that is real time and along with I want to display a graph for that.