obtaining heart rate from ECG signal

39 次查看(过去 30 天)
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

采纳的回答

Diego Caro
Diego Caro 2024-5-12
移动:Image Analyst 2024-5-12
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
  3 个评论
Diego Caro
Diego Caro 2024-5-13
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
Daniel 2024-5-14
I will try the movmean function to get better results.
Once again, thanks!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 ECG / EKG 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by