Calculate RR from ECG signals

I need to calculate RR interval from ECG signals using matlab. Surfing in this forum I find the follow code:

x1 = load('ecg.dat');
x2=x1; fs = 500; % Sampling rate
N = length (x2); % Silength
t = [0:N-1]/fs; % tiidx
figure(1)
subplot(2,1,1)
plot(t,x1)
xlabel('second');ylabel('Volts');title('Input ECG Signal')
% Cancellation DC drift and normalization
x1 = x1 - mean (x1 ); % cancel DC conponents
x1 = x1/ max( abs(x1 )); % normalize to one
subplot(4,1,2)
plot(t,x1)
xlabel('second');ylabel('Volts');title(' ECG Signal after cancellation DC drift and normalization')
% Low Pass Filtering
b=[1 0 0 0 0 0 -2 0 0 0 0 0 1]; a=[1 -2 1];
h_LP=filter(b,a,[1 zeros(1,12)]); % transfer function of LPF
x2 = conv (x1 ,h_LP);
x2 = x2 (6+[1: N]); %cancle delay
x2 = x2/ max( abs(x2 )); % normalize , for convenience .
subplot(4,1,3)
plot([0:length(x2)-1]/fs,x2)
xlabel('second');ylabel('Volts');title(' ECG Signal after LPF')
xlim([0 max(t)])
% High Pass Filtering
b = [-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 32 -32 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];
a = [1 -1];
h_HP=filter(b,a,[1 zeros(1,32)]); % impulse response iof HPF
x3 = conv (x2 ,h_HP);
x3 = x3/ max( abs(x3 ));
subplot(4,1,4)
plot([0:length(x3)-1]/fs,x3)
xlabel('second');ylabel('Volts');title(' ECG Signal after HPF')
xlim([0 max(t)])
% Derivative Filter
% Make impulse response
h = [-1 -2 0 2 1]/8;
% Apply filter
x4 = conv (x3 ,h);
x4 = x4 (2+[1: N]);
x4 = x4/ max( abs(x4 ));
figure(2)
subplot(4,1,1)
plot([0:length(x4)-1]/fs,x4)
xlabel('second');ylabel('Volts');title(' ECG Signal after Derivative')
% Squaring
x5 = x4 .^2;
x5 = x5/ max( abs(x5 ));
subplot(4,1,2)
plot([0:length(x5)-1]/fs,x5)
xlabel('second');ylabel('Volts');title(' ECG Signal Squarting')
% Moving Window Integration
% Make impulse response
h = ones (1 ,31)/31;
Delay = 15; % Delay in samples
% Apply filter
x6 = conv (x5 ,h);
x6 = x6 (15+[1: N]);
x6 = x6/ max( abs(x6 ));
subplot(4,1,3)
plot([0:length(x6)-1]/fs,x6)
xlabel('second');ylabel('Volts');title(' ECG Signal after Averaging')
% Find QRS Points Which it is different than Pan-Tompkins algorithm
max_h = max(x6);
thresh = mean (x6);
poss_reg =(x6>thresh*max_h);
left = find(diff([0 poss_reg])==1);
right = find(diff([poss_reg 0])==-1);
left=left-(6+15); % cancle delay because of LP and HP
right=right-(6+15);% cancle delay because of LP and HP
left=abs(left);
right=abs(right);
for i=1:length(left)
    [R_value(i) R_loc(i)] = max( x1(left(i):right(i)) );
    R_loc(i) = R_loc(i)-1+left(i); % add offset
    
    [Q_value(i) Q_loc(i)] = min( x1(left(i):R_loc(i)) );
    Q_loc(i) = Q_loc(i)-1+left(i); % add offset
    
    [S_value(i) S_loc(i)] = min( x1(left(i):right(i)) );
    S_loc(i) = S_loc(i)-1+left(i); % add offset
end
subplot(4,1,4)
plot (t,x1,t(R_loc) ,R_value , 'r^');
% there is no selective wave
Q_loc=Q_loc(find(Q_loc~=0));
R_loc=R_loc(find(R_loc~=0));
S_loc=S_loc(find(S_loc~=0));
values=x2(R_loc);
% 1 beat/sec x 60 sec/min = 60 beats/min.
%%R-R interval using diffrence operation method
tiscl=t(2)-t(1);
px=diff(R_loc);
d_tms=find(diff(R_loc)>200);
d_men=mean(px(d_tms));
beat=d_men*tiscl;
heart_rate=1/beat*60;
fprintf('%2.1f seconds/beat\n',beat);
rs=sprintf('Heart Rate= %2.1f beats per minit \n',heart_rate);

but if I try to run this code I m getting the follow error:

Error using horzcat
Dimensions of arrays being concatenated are not consistent.

This is the line got to error:

left = find(diff([0 poss_reg])==1);

I m not able to find a solution. I attach also the ecg input file.

Could you help me ?

Answer 1

[Star Strider]

A somewhat more direct approach —

EKG = readmatrix('ecg.txt');
L = numel(EKG)

L = 30000

Fs = 500;                                                       % Samplilng Frequency (Hz)

t = linspace(0, L-1, L)/Fs;                                     % Time Vector (s)

[pks,locs] = findpeaks(EKG, 'MinPeakProminence', 250);

Rlocs = t(locs);
dRlocs = diff(Rlocs);
Rlocs_stats = [mean(dRlocs); median(dRlocs); std(dRlocs); mean(1./dRlocs); median(1./dRlocs); std(1./dRlocs)];

fprintf('\nMean R-R Interval   = %8.3f\nMedian R-R Interval = %8.3f\nStDev R-R Interval  = %8.3f\nMean HR\t\t    = %8.3f\nMedian HR\t    = %8.3f\nStDev HR\t    = %8.3f\n', Rlocs_stats.*[1; 1; 1; 60; 60; 60])

Mean R-R Interval   =    0.604
Median R-R Interval =    0.604
StDev R-R Interval  =    0.090
Mean HR		    =  101.727
Median HR	    =   99.342
StDev HR	    =   16.656

figure
plot(t, EKG, 'DisplayName','EKG')
hold on
plot(t(locs), EKG(locs), '^r', 'DisplayName','R-Deflections')
hold off
grid
legend('Location','best')

Versions of this approach have already appeared elsewhere in MATLAB Answers.

.

[MICHELE CASTRIOTTA]

Thanks for your answer, but, with this code, how can I extract only RR interval RR from input ECG ?

[Star Strider]

It detects and calculates all of them. Choose the one you want!

[MICHELE CASTRIOTTA]

EKG(locs) contains RR interval or t(locs) contains RR interval ?

[Star Strider]

The ‘dRlocs’ vector contains the R-R intervals (in seconds). The rates are in beats/minute.