2
$\begingroup$

I'm working on a little project involving digital filters and DFT, what I want to do is, given a sampled ECG, get the heart rate using power spectral density. I got the samples from this database. Here's the MATLAB code I wrote:

clear;
clc;

fc = 360; % sampling freq [Hz]
T = 1/fc;


Nfile = 10;     % number of files to join
filename = '103m';
ecg_samples = load([filename ' (0)']).val;

for i = 1:Nfile-1
    ecg_samples = [ecg_samples zeros(1,3600)];
    ecg_samples((3600)*i+1:end) = load([filename ' (',num2str(i),').mat']).val;
end

N = length(ecg_samples);    % number of ECG samples
times = (0:N-1)*T;
ecg_samples = detrend(ecg_samples);     % zero mean samples

%{
Creating some filters to clean the noise out of the data:
- high-pass to remove low frequency noise, up to 0.5 Hz
- band-stop to remove A/C power noise, centered on 60 Hz
- low-pass to remove other types of noise
%}
fstopband_hhp = 0.5;  % [Hz]
fpass_hhp = 0.6;    % [Hz]

fstopband_hlp = 155;    % [Hz]
fpass_hlp = 150;    % [Hz]

fstopband1_hbs = 59;    % [Hz]
fstopband2_hbs = 61;    % [Hz]
fpass1_hbs = 58;    % [Hz]
fpass2_hbs = 62;    % [Hz]

hlp_equiripple = designfilt("lowpassfir", ...
    PassbandFrequency = fpass_hlp, ...
    StopbandFrequency = fstopband_hlp, ...
    PassbandRipple = 1, ...
    StopbandAttenuation = 60, ...
    SampleRate = fc);
%fvtool(hlp_equiripple);
delay = mean(grpdelay(hlp_equiripple));

hhp = designfilt("highpassiir", ...
    PassbandFrequency = fpass_hhp, ...
    StopbandFrequency = fstopband_hhp, ...
    PassbandRipple = 1, ...
    StopbandAttenuation = 60, ...
    SampleRate = fc, ...
    DesignMethod = "butter");
%fvtool(hhp);
% Artefacts are present with Butterworth, but not with an elliptic 
% high-pass filter

hbs = designfilt("bandstopiir", ...
    PassbandFrequency1=fpass1_hbs, ...
    StopbandFrequency1=fstopband1_hbs, ...
    StopbandFrequency2=fstopband2_hbs, ...
    PassbandFrequency2=fpass2_hbs, ...
    PassbandRipple1=0.5, ...
    StopbandAttenuation=70, ...
    PassbandRipple2=0.5, ...
    DesignMethod="ellip", ...
    SampleRate=fc);
 %fvtool(hbs);

% Now let's filter the signal
ecg_samples_filtered = filter(hlp_equiripple, [ecg_samples zeros(1,delay)]);
ecg_samples_filtered = ecg_samples_filtered(delay+1:end);   % compensate the linear phase delay

ecg_samples_filtered = filtfilt(hhp,ecg_samples_filtered);

ecg_samples_filtered = filtfilt(hbs,ecg_samples_filtered);



% Spectrum plot
subplot(2,2,1);
plot((0:N-1)/(N*T),abs(fft(ecg_samples)),Color='#D95319');
title("spettro ecg");
xlabel("frequenza [Hz]");
ylabel("modulo |·|");
grid on

subplot(2,2,2);
plot((0:N-1)/(N*T),abs(fft(ecg_samples_filtered)), Color='#0072BD');
title("spettro ecg filtrato");
xlabel("frequenza [Hz]");
ylabel("modulo |·|");
grid on

subplot(2,2,[3,4]);
hold on
plot((0:N-1)/(N*T),abs(fft(ecg_samples)),Color='#D95319');
plot((0:N-1)/(N*T),abs(fft(ecg_samples_filtered)), Color='#0072BD');
title("confrontro tra gli spettri");
xlabel("frequenza [Hz]");
ylabel("modulo |·|");
grid on


% PSD calculation, let's divide the N samples into L sequences of length M
L = 10;
M = N/L;
freq = (0:M-1)/(M*T);
temp = zeros(1,M);

k = -M+1;
for i=1:L
    transf_values = fft(ecg_samples_filtered(k+i*M:i*M));
    mod_values = transf_values.*conj(transf_values);
    temp = temp + mod_values/L;
end

figure;
stem(freq(1:M/2)*60, temp(1:M/2));      % multiplying the x-axis by 60 to convert from Hz to bpm
title("PSD");
xlabel("bpm");
ylabel("power");

[max_value,index] = max(temp(1:M/2));
disp(['Stimated hear rate: ' num2str(freq(index)*60) ' bpm']);

I have some questions: How do I choose which type of filters to use? How can I choose between FIR filters and IIR filters, and which type should I use? (e.g. choosing between cheby2, cheby2, elliptic, butterworth ... ) I know the difference between FIR and IIR, so I also understand that theoretically, using a FIR filter is better due to the linear phase. However, I'm confused about why people on the internet use IIR filters to filter ECG data. In this case, the waveshape is important, so why use a filter that can alter the waveshape? In MATLAB there is the filtfilt() function, but the waveshape is still altered a bit. What is the preferred order for filtering? In which sequence should filtering be performed and why?

I can't figure out why the program only works (properly) with certain sets of ECG samples and settings. For example with the '103m' file, which has 36000 samples with an heart rate of about 66bpm using an elliptic filter as high-pass filter the program says that the heart rate is about 210 bpm, following the plot of both filtered and unfiltered ecg signals: enter image description here

But changing the filter design method to Butterworth gives me the correct heart rate but a big time domain distortion, and I'm not an expert but I don't think a filtered signal should appear like this:

enter image description here

If I use the cheby2 filter design method, the program still shows 210 bpm, but the signal distortions are not too bad: enter image description here

Why do Butterworth filters produce more distortions compared to cheby2 or elliptic filters, even when the pass-band and stop-band frequencies specified in the designfilt() function are the same? Why does the program work with that amount of distortion, but when the filtered data doesn't have all that garbage, it doesn't work anymore? pretty strange. All these graphs were referred to fpass_hhp = 0.6;, but by changing it to fpass_hhp = 0.7; even with a Butterworth filter there are not many distortion but the program still says 210 bpm enter image description here

This is just an example. I encounter this issue with all the sample sets I have. The program works fine with some sample sets, but with others, I need to adjust some parameters to make it work well (with some others doesn't work at all). I've been working on this project for a week, and I still haven't been able to find something that works well with all the data without making any changes in the code. Thanks if someone would like to help me

EDIT 1
To remove the bias I modified the code by adding this line and removing the detrend() function:

ecg_samples = (ecg_samples-ecg_samples(1))./(max(ecg_samples)-ecg_samples(1));

now the first part of the code looks like:

clear;
clc;

fc = 360;
T = 1/fc;


Nfile = 10;
nome_file = '103m';
ecg_samples = load([nome_file ' (0)']).val;

for i = 1:Nfile-1
    ecg_samples = [ecg_samples zeros(1,3600)];
    ecg_samples((3600)*i+1:end) = load([nome_file ' (',num2str(i),').mat']).val;
end

N = length(ecg_samples);
times = (0:N-1)*T;
ecg_samples = (ecg_samples-ecg_samples(1))./(max(ecg_samples)-ecg_samples(1));

But as you can see the problem remains, even without bias: enter image description here enter image description here enter image description here enter image description here

The filter transient seems to not start at the (0,0) point. Now let's examine the signal in the filtering chain. The initial low-pass FIR filtering yields:

enter image description here enter image description here

the low-pass introduces a little offset of about 0.005, now let's examine the high-pass IIR filtering:

enter image description here enter image description here

and finally the band-stop filtering:

enter image description here enter image description here

I believe the transient is present due to the filter chain, where the output signal of one filter is used as the input for another filter. The low-pass filter introduces an offset, affecting not only the low-pass, but also the high-pass and band-stop filters, causing the transient. Do I need to remove the offset from the output signal for each filter before the next filtering operation? I believe that the harmonics, which are multiples of the fundamental frequency and contain higher power, are just noise. Do you think there is a way to estimate the frequency of this noise just by looking at the ECG samples? This would help us get rid of the noise by applying a precise filter, By doing this, we should be able to give the fundamental frequency the highest power peak. Otherwise, I think I'll give up on getting heart rate from PSD and I'll work in the time domain.

EDIT 2 I tried to remove the little 0.005 offset introduced by the lowpass FIR filter but the artefact is still there

EDIT 3 Here's what happens if I run the code on just one signal:

clear;
clc;

fc = 360;
T = 1/fc; 

ecg_samples = load('105m (0).mat').val;
N = length(ecg_samples);
times = (0:N-1)*T;
ecg_samples = (ecg_samples-ecg_samples(1))./(max(ecg_samples)-ecg_samples(1));


fstopband_hhp = 0.5;
fpass_hhp = 0.6;   

fstopband_hlp = 155;
fpass_hlp = 150;  

fstopband1_hbs = 59;
fstopband2_hbs = 61;
fpass1_hbs = 58;
fpass2_hbs = 62;

hlp_equiripple = designfilt("lowpassfir", ...
    PassbandFrequency = fpass_hlp, ...
    StopbandFrequency = fstopband_hlp, ...
    PassbandRipple = 1, ...
    StopbandAttenuation = 60, ...
    SampleRate = fc);
%fvtool(hlp_equiripple);
delay_hlp = mean(grpdelay(hlp_equiripple));


hhp = designfilt("highpassiir", ...
    PassbandFrequency = fpass_hhp, ...
    StopbandFrequency = fstopband_hhp, ...
    PassbandRipple = 1, ...
    StopbandAttenuation = 60, ...
    SampleRate = fc, ...
    DesignMethod = "butter");
%fvtool(hhp);

hbs = designfilt("bandstopiir", ...
    PassbandFrequency1=fpass1_hbs, ...
    StopbandFrequency1=fstopband1_hbs, ...
    StopbandFrequency2=fstopband2_hbs, ...
    PassbandFrequency2=fpass2_hbs, ...
    PassbandRipple1=0.5, ...
    StopbandAttenuation=70, ...
    PassbandRipple2=0.5, ...
    DesignMethod="ellip", ...
    SampleRate=fc);
 %fvtool(hbs);

ecg_samples_filtered = filter(hlp_equiripple, [ecg_samples zeros(1,delay_hlp)]);
ecg_samples_filtered = ecg_samples_filtered(delay_hlp+1:end);
ecg_samples_filtered = filtfilt(hhp,ecg_samples_filtered);
ecg_samples_filtered = filtfilt(hbs,ecg_samples_filtered);


L = 10;
M = N/L;
freq = (0:M-1)/(M*T);
temp = zeros(1,M);

k = -M+1;
for i=1:L
    transf_values = fft(ecg_samples_filtered(k+i*M:i*M));
    mod_values = transf_values.*conj(transf_values);
    temp = temp + mod_values/L;
end

figure;
stem(freq(1:M/2)*60, temp(1:M/2));
title("PSD");
xlabel("bpm");
ylabel("power");

[max_value,index] = max(temp(1:M/2));
disp(['Heart rate: ' num2str(freq(index)*60) ' bpm']);

As you can see by the following pictures the artefacts are still there:

enter image description here

If I use an elliptic high-pass filter instead of a Butterworth high-pass filter there are no artefacts present: enter image description here

$\endgroup$
16
  • $\begingroup$ Are you guys the same person or two people in the same class? Two good questions! $\endgroup$
    – Jdip
    Commented May 19 at 16:51
  • $\begingroup$ None of the two $\endgroup$
    – minghierid
    Commented May 19 at 17:10
  • $\begingroup$ @Jdip No, I'm using Python while he is using MATLAB. Also he is using classical filter while I'm using machine learning to filter. $\endgroup$ Commented May 19 at 17:16
  • $\begingroup$ @minghierid "I need to adjust some parameters to make it work well (with some others doesn't work at all)." You might be interest with deep learning, where the parameters are setting up by model itself. All you have to do is just setting up the hyperparameter (how long the model is training, etc), and sit waiting with coffee. That's entirely black box, all you have to do is telling the model what are the clean signals and what are the noised signals. Here I give you reference to apply the model: github.com/fperdigon/DeepFilter $\endgroup$ Commented May 19 at 17:21
  • 1
    $\begingroup$ First of all: you are concatenating different signals which are clearly not continuous. You end up with gaps or double beats and the bpm varies between 67 and 72 over the set. Second: I tried your code to filter a single signal and it looks just fine. I don't see the artifacts that you are posting. Again: please trim this down to one signal and complete code to illustrate one problem. We can't help you if you make it difficult to reproduce the problem. $\endgroup$
    – Hilmar
    Commented May 21 at 18:09

1 Answer 1

2
$\begingroup$

A few points to discuss

IIR filter are typically minimum phase and FIR are typically linear phase. That's not always the case, but good enough for the current discussion. Linear phase filters are good at maintaining the overall time domain shape. On the downside they delay the signal have "pre-ringing", i.e. an impulse is smeared out both forward and backward in times.

FIR filter often add a lot of bulk delay (latency) and its difficult to do anything meaningful at low frequencies. In your example: having a steep highpass at 0.6Hz with a sample rate of 360 Hz will probably require several 1000 filter taps.

IIR filter maintain causality and have no bulk delay, i.e. impulses stay where that are in time but are smeared out only towards future time. On the downside, they can distort the waveform considerably since the delay is different at all frequencies.

I'm not sure why you get this massive transient with the Butterworth filter. This looks like a state discontinuity and could be a problem in how the filter state was initialized.

Normally Butterworth, cheby1+2 and elliptic filter behave similar. The main difference is that you get steeper transitions with an elliptic than a Butterworth but you need to "pay" for this with limited stopband attenuation and passband ripple. You can get a steeper response with the same order but the resulting post-ringing is pretty much the same. The two cheby flavors are in between Butterworth and Elliptic.

66bpm using an elliptic filter as high-pass filter the program says that the heart rate is about 210 bpm,

Overall your approach seems to be misguided. Your best shot at getting the heart rate is to identify the R-peaks and measure the spacing. You can do a running average of the spacing to get the hear rate as a function of time.

What you are seeing is the third harmonic of the hear rate. A heart rate is fairly spikey wave form with lots of harmonics. There is no reason to assume that the fundamental has he highest energy. In fact, it quite likely hasn't. If you really want to work in the frequency domain, you either need to look at the peak that's in frequency or identify a few prominent and map the frequency to a harmonic grid. Example: if you have prominent peaks at 300bpm, 500bpm & 700bpm, the heart rate would be 100 bpm

EDIT after questions.

What do you mean when you say "state discontinuity"

Let's take a look at the raw data (I'm using data set 4) enter image description here

The data has a large DC bias. The filter output is dependent on the past input samples including the one before your signal starts. If you don't give this information the function assumes the prior input is 0. So what the filter sees looks more like this

enter image description here

Hence you get a massive transient response at the beginning of the signal and you trigger a massive step response. That would be true for any type of filter (FIR or IIR). The workaround is to properly initialize the filter states or in this case, just remove the bias before filtering.

I'm asking you because the exercise is to find the heart rate by filtering and using PSD, no other methods.

That's unfortunate. It's quite easy to get the bpm just by looking at the R peaks.

enter image description here

We certainly can take the FFT of this and take a look at the spectrum. In this case it's easiest to take the FFT of the whole thing to get decent frequency resolution

enter image description here

I've scaled the frequency axis in bpm and penciled in the location of the fist few harmonics from the BPM estimate we got from analyzing the R peaks. And sure enough, these line up pretty well with the peaks we observe in the spectrum.

Note here that the highest energy is NOT at the fundamental (72Hz) but at the 3rd harmonic at about 217 Hz. That's perfectly normal and expected for a peaky waveform like a heartbeat. In order to determine the fundamental bpm you need to identify a good number of prominent peaks and identify the spacing between them. Looking for the maximum won't do.

So this is also peak finding exercise, as it would be in the time domain. But in the frequency domain its a lot harder, partially because the frequency resolution is terrible. In this case you get a point only every 6 bpm. This picture only looks so nice, since you got lucky and the fundamental is close to being a multiple of 6 bpm.

I don't think that filtering would help here at all. The filtering affects only the spectral area you are not looking at anyway. Filtering will improve the time domain signal but that looks pretty good already and you are not allowed to use it anyway, so why bother?

Here is the Matlab code that produced these graphs

%% read the files

% https://data.mendeley.com/datasets/7dybx7wyfn/3
fbase = "C:\Users\hilma\Downloads\7dybx7wyfn-3\ECG signals (1000 fragments)\MLII\1 NSR"
nf = 10;
nx = 3600;
x0 = zeros(nx,nf);
fs = 360;
for i = 1:nf
  p = load(sprintf("%s\\103m (%d).mat",fbase,i-1));
  x0(:,i) = p.val';
end

%% first, let's take out the mean to remove DC bias
x1 = x0 - mean(x0);

%% find the R peaks to determie the "correct" bpm
% this one is a bit sloppy and often misses the last peak and also isn't
% protected against peaks all the way at the edge of the file

timex = (0:nx-1)'/fs;  % time axis
bpm = zeros(nf,1);   % store results
for i = 1:nf
  y1 = x1(:,i);
  % find everything above a certain threshold, 80% seems fine
  ip = find(y1 > 0.8*max(abs(y1)));
  ip = ip(diff(ip) ~= 1);  % eliminate consecutives to islaote peaks
  % zero in in the actual peak
  for k = 1:length(ip)
    [~,imax] = max(y1(ip(k) + (-30:30)));
    ip(k) = ip(k) + imax - 31;
  end
  % find bpm as as mean of the differences
  bpm(i) = fs/mean(diff(ip))*60;
  % plot it
  figure(i); clf;
  plot(timex,y1);
  yy = get(gca,'ylim');
  hold on
  plot(timex(ip),y1(ip),'ro');
  title(sprintf('Data set %d, %6.2f pbm',i-1,bpm(i)));
end
  

%% Try the spectrum

nfft = nx; % FFT length
frbpm = (1:nfft/2)'/nfft*fs*60; % frequency axis in bpm
fx1 = fft(x1);
fx1 = fx1(2:nfft/2+1,:);  % remove DC and conjugate spectrum
for i = 1:nf
  figure(10+i); clf
  plot(frbpm, abs(fx1(:,i)),'LineWidth',2);
  set(gca,'xlim',[0 500]);
  xlabel('bpm ->');
  ylabel('FFR Magnitude');
  yy = get(gca,'ylim');
  % add the first harmonics
  hold on
  for k = 1:8
    plot(bpm(i)*k*[1 1],yy,'k');
  end
  title(sprintf('Data set %d, %6.2f pbm',i-1,bpm(i)));
end
$\endgroup$
4
  • $\begingroup$ What do you mean when you say "state discontinuity" ? $\endgroup$
    – minghierid
    Commented May 20 at 5:33
  • $\begingroup$ Could you explain me this please "If you really want to work in the frequency domain, you either need to look at the peak that's in frequency or identify a few prominent and map the frequency to a harmonic grid. Example: if you have prominent peaks at 300bpm, 500bpm & 700bpm, the heart rate would be 100 bpm" ? $\endgroup$
    – minghierid
    Commented May 20 at 5:35
  • $\begingroup$ "There is no reason to assume that the fundamental has he highest energy" is there a way to obtain the heart rate frequency just with PSD? Maybe by filtering out all the frequency an heart rate can't get, e.g. using a band-pass with band-pass edges at 0.7 Hz and 5 Hz ? I'm asking you because the exercise is to find the heart rate by filtering and using PSD, no other methods. Or maybe by reducing the amplitude of the frequency content outside this range 0,7 - 5 Hz ? I don't know if can work, I'm just asking some random things that popped up in my mind $\endgroup$
    – minghierid
    Commented May 20 at 5:49
  • $\begingroup$ Thank you for your response. I have just edited the answer in order to reply to you. $\endgroup$
    – minghierid
    Commented May 21 at 12:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.