Heart Rate Variability

I would like to plot heart rate variability from an ECG data signal using RR peak detection.

[Ashish Uthama]

Adding a specific MATLAB question with more details might help you get an answer.

[NICOLE MIN]

hi, ive plotted an ECG and etxracted QRS complex using pan tompkin algorithm. how can i plot the HRV signal from the QRS peak extracted?

my code is stated as below:

% pan tompkin algorithm

ecg=(val-0)/200;% extract signal

fs=1000;%sampled frequency

N=length(ecg);% length of signal to extract

t=[0:N-1]/fs; % time period(total sample/Fs )

figure, plot(t,ecg); title('Raw ECG Data plotting ')

xlabel('time')

ylabel('amplitude')

x1=ecg;

fs = 128; % Sampling rate

N = length (x1); % Signal length

t = [0:N-1]/fs; % time index

figure(1)

subplot(2,1,1)

plot(t,x1)

xlabel('second');ylabel('Volts');title('Input ECG Signal')

subplot(2,1,2)

plot(t(200:600),x1(200:600))

xlabel('second');ylabel('Volts');title('Input ECG Signal 1-3 second')

xlim([1 3])

x1 = x1 - mean (x1 ); % cancel DC conponents

x1 = x1/ max( abs(x1 )); % normalize to one

figure(2)

subplot(2,1,1)

plot(t,x1)

xlabel('second');ylabel('Volts');title(' ECG Signal after cancellation DC drift and normalization')

subplot(2,1,2)

plot(t(200:600),x1(200:600))

xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')

xlim([1 3])

% LPF (1-z^-6)^2/(1-z^-1)^2

b=[1 0 0 0 0 0 -2 0 0 0 0 0 1];

a=[1 -2 1];

%high pass filter

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 .

figure(3)

subplot(2,1,1)

plot([0:length(x2)-1]/fs,x2)

xlabel('second');ylabel('Volts');title(' ECG Signal after LPF')

xlim([0 max(t)])

subplot(2,1,2)

plot(t(200:600),x2(200:600))

xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')

xlim([1 3])

% HPF = Allpass-(Lowpass) = z^-16-[(1-z^-32)/(1-z^-1)]

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 (16+[1: N]); %cancle delay

x3 = x3/ max( abs(x3 ));

figure(4)

subplot(2,1,1)

plot([0:length(x3)-1]/fs,x3)

xlabel('second');ylabel('Volts');title(' ECG Signal after HPF')

xlim([0 max(t)])

subplot(2,1,2)

plot(t(200:600),x3(200:600))

xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')

xlim([1 3])

% 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(5)

subplot(2,1,1)

plot([0:length(x4)-1]/fs,x4)

xlabel('second');ylabel('Volts');title(' ECG Signal after Derivative')

subplot(2,1,2)

plot(t(200:600),x4(200:600))

xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')

xlim([1 3])

x5 = x4 .^2;

x5 = x5/ max( abs(x5 ));

figure(6)

subplot(2,1,1)

plot([0:length(x5)-1]/fs,x5)

xlabel('second');ylabel('Volts');title(' ECG Signal Squarting')

subplot(2,1,2)

plot(t(200:600),x5(200:600))

xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')

xlim([1 3])

% 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 ));

figure(7)

subplot(2,1,1)

plot([0:length(x6)-1]/fs,x6)

xlabel('second');ylabel('Volts');title(' ECG Signal after Averaging')

subplot(2,1,2)

plot(t(200:600),x6(200:600))

xlabel('second');ylabel('Volts');title(' ECG Signal 1-3 second')

xlim([1 3])

figure(7)

subplot(2,1,1)

max_h = max(x6);

thresh = mean (x6 );

P_G= (x6>0.01);

difsig=diff(P_G);

figure (8)

subplot(2,1,1)

hold on

plot (t(200:600),x1(200:600)/max(x1))

box on

xlabel('second');ylabel('Integrated')

xlim([1 3])

subplot(2,1,2)

plot (t(200:600),x6(200:600)/max(x6))

xlabel('second');ylabel('Integrated')

xlim([1 3])

left=find(difsig==1);

raight=find(difsig==-1);

left=left-(6+16);

raight=raight-(6+16);

for i=1:length(left)-1

[R_value(i) R_loc(i)] = max( x1(left(i):raight(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):raight(i)) );

S_loc(i) = S_loc(i)-1+left(i); % add offset

[P_value(i) P_loc(i)] = min( x1(left(i):raight(i)) );

P_loc(i) = P_loc(i)-1+left(i); % add offset

end

Q_loc=Q_loc(find(Q_loc~=0));

R_loc=R_loc(find(R_loc~=0));

S_loc=S_loc(find(S_loc~=0));

P_loc=P_loc(find(P_loc~=0));

figure

subplot(2,1,1)

title('ECG Signal with R points');

plot (t,x1/max(x1) , t(R_loc) ,R_value , 'r^', t(S_loc) ,S_value, '*',t(Q_loc) , Q_value, 'o');

legend('ECG','R','S','Q');

subplot(2,1,2)

plot (t,x1/max(x1) , t(R_loc) ,R_value , 'r^', t(S_loc) ,S_value, '*',t(Q_loc) , Q_value, 'o');

xlim([1 3])

thanks