# How to compute the HRV and peak extraction in ECG?

21 ビュー (過去 30 日間)
HelpAStudent 2023 年 6 月 17 日

Can you help me to compute the HRV and detection of the peaks from this code:
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');

サインインしてコメントする。

### 回答 (1 件)

Varun 2023 年 8 月 30 日
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 findpeakslooks 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.
Hope this helps!
Thanks,
Varun

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeScopes and Data Logging についてさらに検索

R2021a

### Community Treasure Hunt

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

Start Hunting!

Translated by