How to compute the HRV and peak extraction in ECG?

Can you help me to compute the HRV and detection of the peaks from this code:

load sig_ecg.mat

Fs = 300;
L = numel(sig);
t = linspace(0, L-1, L)/Fs;
% t=linspace(0,length(sig)/fs,length(sig));
figure()
plot(t,sig/max(sig));
grid
title('Original Signal in Time Domain');
xlabel('Time(s)');
ylabel('Amplitude');

figure()
sig2=sig-mean(sig); % Per visualizzare meglio lo spettro
f=linspace(-Fs/2,Fs/2,length(sig2));
plot(f,fftshift(abs(fft(sig2))));
grid
title('Signal in Frequency Domain'); 
xlabel('Frequency (Hz)');
ylabel('Amplitude');

Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
NFFT = 2^nextpow2(L)

sig = bandstop(sig, [48 52], Fs, 'ImpulseResponse','iir');  % Remove 50 Hz Mains Interference
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Ws = 0.5/Fn;                                                % Passband Frequency Vector (Normalised)
Wp = 2.5/Fn;                                                % Stopband Frequency Vector (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs = 150;                                                   % Stopband Attenuation (dB)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                             % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp,'high');                         % Default Here Is A Lowpass Filter
[sos,g] = zp2sos(z,p,k);                                    % Use Second-Order-Section Implementation For Stability
sig_filtered = filtfilt(sos,g,sig);                           % Filter Signal (Here: isrg)
figure
freqz(sos, 2^14, Fs)                                        % Bode Plot Of Filter
set(subplot(2,1,1), 'XLim',[0 50])                          % Optional, Change Limits As Necessary
set(subplot(2,1,2), 'XLim',[0 50])                          % Optional, Change Limits As Necessary

figure
plot(t, sig_filtered)
grid
title('Filtered Signal in Time Domain');
xlabel('Time(s)');
ylabel('Amplitude');

figure()
sig=sig-mean(sig); % Per visualizzare meglio lo spettro
f=linspace(-Fs/2,Fs/2,length(sig));
plot(f,fftshift(abs(fft(sig))));
grid
title('Filtered Signal in Frequency Domain'); 
xlabel('Frequency (Hz)');
ylabel('Amplitude');

[wel,f_wel] = pwelch(sig, hamming(1000),250,1000,Fs); % Finestra Hamming di 1000 campioni con overlap di 1/4 e no zero padding

figure()
plot(f_wel,wel)
grid
title('PSD Welch');

Calcolo dello spettro attraverso il metodo Periodogram:
[per_rect, f_per_rect] = periodogram(sig,[],1000,Fs);
plot(f_per_rect,per_rect);
grid
title('PSD Periodogram with Rectangular Window')
figure;

[per_ham, f_per_ham] = periodogram(sig,hamming(length(sig)),1000,Fs);
plot(f_per_ham,per_ham);
grid
title('PSD Periodogram with Hamming Window')
figure;


% Yule-Walker
order=50; 
ar_yule=zeros(order,order+1);
for i=1:order
    [ar_yule(i,1:i+1),e_ord(i), rc_ord]= aryule(sig,i);
end 


% Yule
figure()
subplot(1,2,1)
stem(rc_ord);
title('Reflection coefficients')
subplot(1,2,2)
plot(e_ord);
title('Variance');

ordine = 9; % The optimal order is chosen by the behavior of the reflection coefficients (they should go to zero) and of the variance.
[Pxx_y,f_y] = pyulear(sig,ordine,Fs);
figure()
plot(f_y*(Fs/2)/pi,Pxx_y); % This is done because the output frequencies from pyulear are normalized frequencies.
grid
title('Yule-Walker')
xlabel('Frequency [Hz]'); 
ylabel('dB/Hz');

% Burg
ord = 50;
ar_burg=zeros(ord,ord+1);
for i=1:ord
    [ar_burg(i,1:i+1),eb_ord(i),rcb_ord]=arburg(sig,i);
end
figure()
subplot(1,2,1)
stem(rcb_ord);
title('Reflection coefficients');
subplot(1,2,2)
plot(eb_ord);
grid
title('Variance');

ordine_b = 9;
figure()
[Pxx_b,f_b]=pburg(sig,ordine_b,Fs);
plot(f_b*(Fs/2)/pi,Pxx_b,'r');
grid
title('Burg')
xlabel('Frequency [Hz]'); 
ylabel('dB/Hz');
hold on;

Below is how I tried to compute the HRV:

th=0.8;
[RR,N] = peak_detection(sig,Fs,th);

taco=diff(RR)/Fs;
% inverso intervallo è la nuova freq di campionamento
figure()
plot(taco);
grid
title('Tachogram');
xlabel('Beats');
ylabel('RR(s)');

taco=taco(15:end);

fs1=1/(abs((min(taco)))); 
taco=taco-mean(taco);

% Yule-Walker
ordine=60;
a_ord=zeros(ordine,ordine+1);
for i=1:ordine    
[a_ord(i,1:i+1),e_ord(i),rc]=aryule(taco,i);
end

plot(rc);
title('A COEFFICIENT:Pick the one close to 0');

figure()
plot(e_ord);
title('VARIANCE:the value before and after has to be similar')

figure()
o=7;
[Pxx_t,f_t]=pyulear(taco,o,length(taco),fs1);
%plot(f_t,10*log10(Pxx_t));
plot(f_t,(Pxx_t));
title ('HRV using Yule-Walker');
xlabel('frequency (Hz)'); 
ylabel('dB/Hz'); 
legend(num2str(o));
title(legend,'order');

% Burg
figure()
ordine=50; 
for i=1:ordine
[a_ord,e_ord(i),rc]=arburg(taco,i); 
end 
plot(rc);
title('A COEFFICIENT:Pick the one close to 0');

figure()
o=7;
[Pxx_b,f_b]=pburg(taco,o,length(taco),fs1);
plot(f_b,Pxx_b,'r');
title('HRV using Burg')
xlabel('frequency (Hz)'); 
ylabel('dB/Hz'); 
legend(num2str(o));
title(legend,'order');


figure()
plot(f_b,Pxx_b,'r'); hold on;
plot(f_t,Pxx_t,'b');
title('HRV')
xlabel('frequency (Hz)'); 
ylabel('dB/Hz'); 
legend('Burg','Yule-Walker');
title(legend,'AR-Methods');

Answer 1

[Varun]

Hello!

I ran your code on my end but without the ‘peak_detection’ method defined in the provided code, I couldn’t figure out where the issues lies. I also saw that you are attempting to compute the HRV using ‘Yule-Walker’ and ‘Burg’ methods. I am assuming that the method ‘peak_detection’ is the main issue.

I recommend using the ‘findpeaks’ function in MATLAB for the same. The signal ‘sig’ and sample rate ‘Fs’ can be passed to the method as input arguments, while the threshold value ‘th’ can be passed as a Name-Value argument, for ‘MinPeakHeight’. So, the modified call for function ‘findpeaks’ looks something like this:

[ pks, locs ]=findpeaks(sig, Fs, 'MinPeakHeight', th); 

After this, the RR intervals can be calculated by applying the ‘diff’ function on ‘locs’.

You may refer to the following documentation for more information:

Findpeaks function: https://www.mathworks.com/help/signal/ref/findpeaks.html

Peak Analysis: https://www.mathworks.com/help/signal/ug/peak-analysis.html

Diff function: https://www.mathworks.com/help/matlab/ref/diff.html

Hope this helps!

Thanks,

Varun