obtaining heart rate from ECG signal

Hello,

I need to obtain heart rate values from an ECG signal. While searching for information, I discovered that I need to (more or less) get the RR interval and then use the time difference between samples to calculate the instantaneous HR.

I have a database acquired in the local hospital, but the raw signal is full of noise and invalid samples but I dont know how to remove both.

Attached you can find of the raw ECG. It is sampled at 1kHz, stored in int16 data type and the measured units are milivolts.

fl = fopen("ficheroAnalogico.dat"); 
A = fread(fl, inf, 'int16');                
fclose(fl); 

figure(1)
subplot(3,1,1); 
plot(A);
title('raw');

subplot(3,1,2);
plot(A_4000);
title('zoom');

I use this short code to read the data and plot it, but, as you can see in the capute, the data is full of noise and invalid samples.

The second plot show the first 4000 samples, where you can see the R peaks, but I need to filter and remove invalid samples before calculating HR. Can you please help me? I know I'm asking for too much, but I haven't programmed in more than 10 years, so my knowlegde is rusted and I'm really stucked.

Thanks you in advance and sorry for my english.

Daniel

[Mario Malic]

Amazing, put this as an answer @Diego please.

[Diego Caro]

Done!

Answer 1

[Diego Caro]

If you only want to estimate the heart rate, my suggestion is that you take the derivative of the ecg signal, so that R peaks increase their amplitude considerably.

Getting the raw signal:

fl = fopen("ficheroAnalogico.dat");
A = fread(fl,inf,'int16');
fs = 1e3; % Setting Sample Frequency
fclose(fl);

A = A - mean(A); % Removing offset
n = length(A);
t = (1:n)/fs; 

figure
plot(t,A)
title('Raw data')
xlabel('Time (s)')
ylabel('Amplitude (mV)')

Now getting the derivative:

A_der = A;
for k = 2:length(A)-1         
    A_der(k) = (A(k+1)-A(k))./(1/fs);
end 
figure
plot(t,A_der)
title('Derivative of raw data')
xlabel('Time (s)')
ylabel('Amplitude (mV)')

If you zoom into the figure above (R peaks show consistently as negative spikes):

Then you could use findpeaks to get RR intervals:

[peaks,tpeaks] = findpeaks(-A_der,'MinPeakHeight',1.5,'MinPeakDistance',0.25*fs);
hold on
plot(t(tpeaks),-peaks,'xr')
hold off

Zooming in to verify proper detection of negative peaks:

xlim([0 60])

Finally, computing the heart rate from the average RR interval (AVNN)

NN = 0;
j = 1;
peak_count = length(peaks);
NNs = zeros(1,peak_count-1);

for i = 1:peak_count-1
    NNi = tpeaks(j+1) - tpeaks(i);
    NNs(i) = NNi;
    NN = NN + NNi;
    j = j + 1;
end

AVNN = NN/(fs*length(peaks));
BPM = 60/AVNN

BPM =

152.4790

[Daniel]

Thank you so much for your answer Diego. Is what I needed.

Now I want to do one more step and calculate the FHR each two R peaks, like the "instantaneous FHR" (I don't know if this term is correct). Using the peaks and tpeaks values, I can do this:

% Inicializar el vector vectorAnalogFHR
vectorAnalogFHR = zeros(size(tpeaks));
diferencia_absoluta = zeros(size(tpeaks));

% Calcular la diferencia absoluta entre muestras consecutivas del vector que 
% contiene las posiciones de las ondas R (distancia RR)y dividirlas por 
% 60000 para obtener la FHR en cada instante
for i = 2:length(tpeaks)
    diferencia_absoluta(i) = abs(tpeaks(i) - tpeaks(i-1));
    B(i) = 60000 / diferencia_absoluta(i);
end

Then, I can plot the all the signals. First the raw data and the derivated:

Second the signal with their peaks detected, the values of "instantaneous FHR" and HR2 (this signal is the FHR values obtained directly from the cardiotocographic monitor).

My objetive is that the FHR values obtained from raw ECG (valores HR, second plot) looks like HR2 (third plot). Is almost done, but as you can see in the graph, near the end there are some incorrect values, probably due to the gap between marked peaks in that part of the signal (I used the zoom on the first plot of the graph to show that). Why is this happening? Because of the offset removal? Or due to the derivation of the signal?

Once again, thank you so much for your first answer and I hope you can help my with this too.

Daniel

[Diego Caro]

In fact, the instantaneous BPM you obtained is correct for the signal provided. If you look at whats happening with the raw signal at the moment the derivative goes flat, you see the raw signal is also flat.

If you want to obtain a figure like HR2, I'd suggest using a moving average filter with movmean, to smooth out the spike at the end of the plot you obtained.

BPM2 = movmean(B,5); %try different values for k 
tt = max(t)*(0:length(B)-1)/length(B);
figure
plot(tt,BPM2)
title('Smooth HR')
xlabel('Time (s)')
ylabel('BPM')
ylim([0 200])

[Daniel]

I will try the movmean function to get better results.

Once again, thanks!