ECG signal filtering problem

29 views (last 30 days)
M
M on 29 Mar 2021
Commented: Star Strider on 29 Mar 2021
This is my ECG signal I wanna improve. Firstly i would like to cut and make high frequencies kind of equal:
Secondly I would like to smooth the noise below.
How to create a filter like this using filterdesign? I found some parameters from other questions concerning it but after I put it in filterdesign, it was destroying my signal like doubling its amplitude etc. I'd be grateful for any help

Accepted Answer

Mathieu NOE
Mathieu NOE on 29 Mar 2021
hello
a little tool to do the fft anamysis (and notch filter demo included)
clc
clearvars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% dummy data
Fs = 10000;
samples = 20000;
t = (0:samples-1)'*1/Fs;
signal = sin(2*pi*50*t)+2*sin(2*pi*440*t)+1e-2*randn(samples,1); % two sine + some noise
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 5;
if decim>1
signal = decimate(signal,decim);
Fs = Fs/decim;
end
samples = length(signal);
t = (0:samples-1)*1/Fs;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = Fs; % so delta f = 1 Hz
Overlap = 0.75;
w = hanning(NFFT); % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
fc = 50; % notch freq
wc = 2*pi*fc;
Q = 5; % adjust Q factor for wider (low Q) / narrower (high Q) notch
% at f = fc the filter has gain = 0
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
b0 = 1;
b1 = -2*cos(w0);
b2 = 1;
a0 = 1 + alpha;
a1 = -2*cos(w0);
a2 = 1 - alpha;
% analog notch (for info)
num1=[1/wc^2 0 1];
den1=[1/wc^2 1/(wc*Q) 1];
% digital notch (for info)
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = (fc-1:0.01:fc+1);
[g1,p1] = bode(num1,den1,2*pi*freq);
g1db = 20*log10(g1+1e-6);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1+1e-6);
figure(1);
plot(freq,g1db,'b',freq,gd1db,'+r');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
legend('analog','digital');
xlabel('Fréquence (Hz)');
ylabel(' dB')
% now let's filter the signal
signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap); %
sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)
% demo findpeaks
[pks,locs]= findpeaks(sensor_spectrum_dB,'SortStr','descend','NPeaks',2);
[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),plot(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')
text(freq(locs)+1,pks,num2str(freq(locs)))
function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).
% Linear averaging
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,channels);
s_tmp((1:samples)) = signal;
signal = s_tmp;
samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
% compute fft with overlap
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
% % for info is equivalent to :
% noverlap = Overlap*nfft;
% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows
% main loop
fft_spectrum = 0;
for i=1:spectnum
start = (i-1)*offset;
sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only
end
fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

More Answers (1)

Star Strider
Star Strider on 29 Mar 2021
The signal appears to have two problems, baseline variation and high-frequency noise.
Ideally, it would be best to first calculate the fft of the signal in order to determine where the high-frequency noise begins, and if it is a specific frequency or a band of frequencies.
Otherwise, use the bandpass function with a lower passband of 1 Hz and an upper passband of between 45 Hz and 100 Hz. If the high-frequency noise is at a specific frequency, experiment with a bandstop filter to eliminate it specifically, and a highpass filter to eliminate the baseline variation. Other filter designs and approaches are also possible, depending on the signal what you want to do.
  2 Comments
M
M on 29 Mar 2021
I understand that it should be IIR butterworth filter? What about such parameters as Fstop1/Fstop2 or Astop1/Astop2?
Star Strider
Star Strider on 29 Mar 2021
I prefer elliptic IIR filters, since they are more computationally efficient. The functions I linked to take care of all the design details, so it is only necessary to specify the stopband frequencies and the sampling frequency.
It is straightforward to design any number of filters with the Signal Processing Toolbox, so if you want to design your own filters from command-line functions, I can show you how to do that. It is more difficult to demonstrate interactive functions such as DesignFilt so I will leave that to you to explore.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!