Calculate LF HF for ECG using Welchs PSD

I have an ECG timeseries that I am extracting time and frequency metrics.

For the time metrics (SDNN, RMSSD, etc) I am using a window of 5 minutes and a step of 30 sec, so there is an overlap of 90% in every sequential window. I want to calculate the high and low frequency power in each of these windows as I want to track the changes of the different frequencies with time. I can not upload the data but you can recreate it with the following code.

However I am not sure how to procced with the Welch's PSD. In the help page there is [pxx,f]=pwelch(x,window,nooverlap,f,fs).

What is 'nooverlap' if I just want to calculate for each window of 5 minutes as I have defined below (bwin). Also is 'window' the window size I am using in this case 5 minutes or that is something else? The ultimate goal is to calculate LF and HF but I am not sure how to procceed. I was wondering if anyone can help.

fs=1000;
t=[0:1/fs:3000];
x=zeros(length(t),1);
for i=1:5
x=x+i*sin(2*pi*2*i*t');
end

plot(t,x)
[peak loc]=findpeaks(x,'MinPeakProminence',12);
hold on
plot(t(loc), x(loc),'*r')
xlim([0 20])
hold off

wind=300*fs;
step=30*fs;
p=[1:step:length(t)-wind];
for i=1:length(p)
    bwin{i}=[p(i):p(i)+wind-1];
  ecg=x(bwin{i});
  [peak loc]=findpeaks(ecg,'MinPeakProminence',12);
  RR{i}=diff(loc);
end

Answer 1

[Mathieu NOE]

hello

this would be my suggestion

the nfft is choosen based on what frequency resolution (df) you need

remember that df = fs/nfft , so the higher nfft the finer the resolution

asyou generate a signal with 5 tones with a spacing of 2 Hz, a df below 0.5 hz should suffice to distinguish the frequencies in the pwelch spectrum. see my options in the code below.

,ow you can go up to the limit and use nfft = number of samples (step) and you get the finest resolution, but then there is no overlap possible.

overlap can get you a smoother output when your data are a bt corrupted with noise, at the price of a lower frequency resolution (there is no free lunch here either !) so it's up to you to play with the parameters (on real signals) and see if you have to improve the frequency resolution (increase nfft) or improve the noise rejection (increased overlap)

I have a bit modified the signal itself so its frequency is slightly going up with time (chirp like) , so you would actually see in the last plot that the LF and HF components are following that trend.

hope it helps !

fs=1000;
t=[0:1/fs:3000]';
x = 0;

for k=1:5
    x=x+k*sin(2*pi*2*k*t + 0.01*t.^2); % nb the signal frequency increases with time
end

figure(1);
plot(t,x)

[peak loc]=findpeaks(x,'MinPeakProminence',12);
hold on
plot(t(loc), x(loc),'*r')
xlim([0 20])
hold off

wind=300*fs;
step=30*fs;
p=[1:step:length(t)-wind];

% fft (pwelch) parameters
% "good enough" resolution (allows some overlapping for noise reduction - beneficial on real signals)
NFFT = 5000; % then the fft spectrum will have freq resolution df = 0.2 Hz (because df = fs/NFFT)
NOVERLAP = round(0.5*NFFT); % typical values are 0.5, 0.75 of NFFT
WINDOW = hanning(NFFT);

% % maximal resolution option  
% NFFT = step; % then the fft spectrum will have freq resolution df = 0.0333 Hz (because df = fs/NFFT)
% NOVERLAP = 0; % typical values are 0.5, 0.75 of NFFT
% WINDOW = hanning(NFFT);


figure(2);hold on
xlim([0 30])
title('Overlaid FFT Pwelch Spectra');

for i=1:length(p)
    bwin=[p(i):p(i)+wind-1];
    ecg=x(bwin);
    % fft (pwelch)
    [Pxx,F] = pwelch(ecg,WINDOW,NOVERLAP,NFFT,fs);
    plot(F,Pxx); 
  
    [peak loc]=findpeaks(Pxx,'MinPeakProminence',0.25);
    LF(i)=min(F(loc));  % LF is first (lowest) identified frequency 
    HF(i)=max(F(loc));  % HF is last (highest) identified frequency     
    tc(i) = mean(t(bwin)); % time index centered in this data (buffer)
end
hold off

% plot LF, HF vs time
figure(3);
plot(tc,LF,tc,HF)
title('LF, HF vs time');
xlabel('Time (s)');
ylabel('Freq (Hz)');