How can I calculate heart rates?

I have a ".txt" file obtained by an ECG, now I need calculate the cardiac frequencies from that ECG. I've built some code but I'm pretty confused, to be able to do the rest. Someone can help me?

clear all
close all
clc
 Fs=1000;
%t=1/Fs;
%N = length (x2)
file =importdata('ecg16.txt'); 
save('ecg16.mat','-struct','file'); 
load('ecg16.mat'); 
file1=importdata('ecg16.mat'); 
sinal=file1.data(:,4); 
ecg16=sinal(60001:70000);
t=1/Fs:1/Fs:length(ecg16)/Fs; 
figure;
subplot (2,1,1);
plot(t,ecg16); 
xlabel('t/s');
title ('ECG without filter')
y = sgolayfilt(ecg16,3,41);  
subplot (2,1,2)
%plot(60000:70000,y(60000:70000));
plot (1: 2000, y (1: 2000))
xlabel('t/s')
title ('ECG with filter')
%eixo ([2000 -4 0 4])
%grade
T=0.001:0.001:10;
figure,
plot(T,ecg16);
xlabel('T');
title ('interval ECG');
F=1./T;
figure,
plot(T,F);
xlabel('T');
ylabel('F');
title ('Frequency do ECG');
n=60;
l=length(ecg16);
if n>=l 
    r=data-mean(ecg16)*ones(l,1);
else
    r=zeros(l,1);
    for i=1:1:n+1; 
        r(i)=ecg16(i)-mean(ecg16(1:n+i));
    end
    for i=l-n+1:1:l; 
        r(i)=ecg16(i)-mean(ecg16(i-n:l));
    end
    for i=n+2:1:l-n; 
        r(i)=ecg16(i)-mean(ecg16(i-n:i+n));
    end
end
figure; plot(t,r); 
xlabel('t/s');
title('Alinhamento dos picos')
s=(r>0.7*max(r)); 
i=1;
picos=[];
while i<=length(s)
    if s(i)==0 %ignorar os 0s
        i=i+1;
    else
        v=[];
        while s(i)==1 
            v=[v i];
            i=i+1;
        end
        [y j]=max(r(v)); 
        picos=[picos [y;v(j)]]; 
    end
end
figure; hold on
plot(picos(2,:)/Fs,picos(1,:),'.r');
plot(t,r);
xlabel('t/s'); hold off;
title('')
freq = [];
for i=2:length(picos) 
    freq = [freq ((picos(2,i)-picos(2,i-1))/Fs)^-1]; %f=1/(T(i)-T(i-1))
end
x=picos(2,2:end);
y=freq;
xi = picos(2,2):picos(2,end);
yi = interp1(x,y,xi,'linear'); 
figure;
hold on;
plot(x/Fs,y,'or');
plot(xi/Fs,yi);
xlabel('t/s');
ylabel('f/Hz');
title('frequency/ time')
hold off;
soma=zeros(1,600); 
amostras=zeros(47,600); 
for i=1:length(picos)
    soma=soma+r(picos(2,i)-300:picos(2,i)+299)';
    amostras(i,:)=r(picos(2,i)-300:picos(2,i)+299)';
end
soma=soma/length(picos); 
figure;
plot(1/Fs:1/Fs:length(soma)/Fs,soma);
xlabel('t/s');
title ('Onda PQRST')

[Walter Roberson]

You need to indicate what you are confused about, or what in your code is not working the way you expect.

[Biza Ferreira]

I need calculate the frequency of this ECG acquisition, the file was obtained from a real ambulatory ECG

Answer 1

[Star Strider]

I will help you out a bit with the signal processing, and I will let you do the heart rate and all other necessary calculations:

fidi = fopen('Biza Ferreira ecg16.txt', 'rt');
Data = textscan(fidi, repmat('%f', 1, 11), 'Delimiter','\t', 'CollectOutput', 1, 'HeaderLines',8);
fclose(fidi);
Fs = 1000;                                                          % Sampling Frequency (Hz)
Fn = Fs/2;                                                          % Nyquist Frequency
Ts = 1/Fs;                                                          % Sampling Interval (sec)
EKG = Data{:}(:,4);                                                 % EKG Data
T = [0:size(EKG,1)-1]*Ts;                                           % Time Vector

Rp =  5;                                                            % PAssband Ripple (dB)
Rs = 25;                                                            % Stopband Ripple (dB)
Wp = [2  50]/Fn;                                                    % Normalised Passband Frequencies
Ws = Wp.*[0.2  1/0.2];                                              % Normalised Stopband Frequencies

[n,Wn]  = buttord(Wp, Ws, Rp, Rs);                                  % Calculate Butterworth Filter Order
[b,a]   = butter(n,Wn);                                             % Calculate Butterworth Filter Coefficients
[sos,g] = tf2sos(b,a);                                              % Convert To ‘SOS’ For Stability
EKGf = filtfilt(sos,g,EKG);                                         % Filter Signal

figure(1)                                                           % Filter Characteristic (Bode Plot)
freqz(sos, 1024, Fs);

figure(2)                                                           % Examine Fidelity Of Filtered EKG Here
plot(T, EKGf)
hold on
plot(T, EKG)
hold off
grid
axis([0  2.5    ylim])                                              % Comment This Out To See The Entire Trace
legend('Original', 'Bandpass Filtered')
title('Filtered and Unfiltered EKG')
xlabel('Time (sec)')
ylabel('Amplitude (mV)')

[R, Rt] = findpeaks(EKGf, Fs, 'MinPeakHeight',500);                 % Find R-Waves & Times

figure(3)                                                           % Plot Filtered EKG
plot(T, EKGf)
hold on
plot(Rt, R, '^r')                                                   % Plot Identified R-Waves
hold off
grid
axis([0  2.5    ylim])                                              % Comment This Out To See The Entire Trace
legend('EKG', 'R-Wave Peaks')
title('Filtered EKG With R-Wave Locations')
xlabel('Time (sec)')
ylabel('Amplitude (mV)')

It is generally necessary to filter an EKG signal to remove the wandering baseline (usually due to respiratory motion) and high-frequency noise. This EKG had a bit of high-frequency noise, so I used a Butterworth bandpass filter to filter the noise and the baseline drift to get a stable baseline. I then used the findpeaks function to locate the R-waves and their times, since in this EKG it will work for that purpose. The best way to understand this code is to read the documentation on the various functions to understand what they do and how I used them here.

ABNORMAL EKG I am not certain where you got this EKG, and I am not certain what lead it is, but it is distinctly abnormal. The ST segment should be isoelectric in a healthy EKG, and there are about 2 mV of ST segment depression in this one. There is also a degree of electrical alternans in the R-wave amplitudes. The intervals and amplitudes are otherwise within normal limits. The abnormalities could be normal variants, but both need to be investigated with a full workup.

[Biza Ferreira]

 Star Strider, it's a great help. Can you provide-me some websites or
 advise some documentation to read more about this. This data was
 acquired with a bioplux ambulatory ecg device

[Star Strider]

The best way to read about the functions is to read the documentation in the Signal Processing Toolbox.

You will need to review the documentation on the Bioplux recorder to understand the recording and filtering it uses.

With respect to reading about cardiovascular medicine, your best resource is Braunwald's Heart Disease: A Textbook of Cardiovascular Medicine. A somewhat more condensed discussion on the electrocardiogram itself is ‘Basic ECG Theory, 12-Lead Recordings and Their Interpretation’. The instrumentation involving EKG acquisition is discussed in Webster’s Medical Instrumentation: Application and Design, 4th Ed. and other biomedical engineering textbooks.

[Biza Ferreira]

Thanks for your concern. One last thing how can I count the peaks in your code?!

[Star Strider]

My pleasure.

The peak values are in the ‘R’ vector and their locations in time in ‘Rt’ (in seconds), so the number of R-waves in the recording are the lengths of either vector.