62 views (last 30 days)

Show older comments

Rp = 1; % passband ripple

Rs = 40; % stopband ripple

Fs = 5000; % sampling frequency

Fn = Fs/2; % nyquist frequency

Fpass = [0.05 150]; % passband

Fstop = [0.01 200]; % stopband

Wp = Fpass/Fn; % normalized

Ws = Fstop/Fn;

% Create a butterworth filter

[N,Wn] = buttord(Wp,Ws,Rp,Rs); % Returns the lowest order

[B,A] = butter(N,Wn); % Butterworth filter design

[sos1,g1] = tf2sos(B,A); % Convert digital filter transfer to second-order section form

ecg = filtfilt(sos1, g1, ecg_raw);

time = 1:length(ecg);

subplot(2,1,1)

plot(time,ecg_raw);

title('Raw ecg')

subplot(2,1,2)

plot(time,ecg)

title('Ecg filtered')

Hi.

Im trying to filter a ecg signal using a butterworth filter. Been looking a lot of postes trying to figure it out. Im new in signal in signal processing so it's possible that i dont fully understand the processing and the parametres which is required doing it.

can you point me the right direction?

Star Strider
on 27 Apr 2021

The filter may be unstable, and if so, may not produce any output. Another possibility is that one or more of the elements of the data are NaN (although unlikely).

Fs = 5000;

Fpass = [1 100];

ecg = bandpass(ecg_raw, Fpass, Fs);

The signal appears to have a significant noise component, however it is not possible to determine what that is. The way to determine that is to calculate and plot the Fourier transform of the signal —

ecg_raw = data(:,3);

L = numel(ecg_raw);

time = linspace(0, L, L)/Fs;

N = 2^nextpow2(L);

FTs = fft(ecg_raw,N)/L;

Fv = linspace(0, 1, N/2+1)*Fn;

Iv = 1:numel(Fv);

figure

plot(Fv, abs(FTs(Iv))*2)

grid

If this is 50-60 Hz supply frequency noise, it can be eliminated with a notch filter (the bandstop function). If it is broadband noise, the only effective ways to deal with it are wavelet denoising or the Savitzky-Golay filter (the sgolayfilt function).

Star Strider
on 28 Apr 2021

I was in the middle of responding when Windows crashed. Again. (I really hate the current version of Windows!)

Anyway, since it is not possible to load .mat files in the online application that I now routinely use to run my code in my Answers and Comments, I downloaded the .mat file and ran it on this computer, offline. I initially forgot to copy the ‘Fn’ assignment, however it is now added where it should be.

The last ‘(Detail)’ figure is specifically for the provided data. Note that it has a xlim call at the end, and this is likely the reason a plot using a different data set is empty, since it may not contain the same x-axis values. (The last figure is not necessary for the code. It simply demonstrates that the FIR filter works even better than I though it would!)

The filter itself is straightforward. A detailed discussion of it is in the kaiserord documentation. Beyond that, it is simply a standard FIR filter with the advantage that it combines a highpass filter (eliminiating the baseline variation) and two stopband filters to eliminate the mains frequency interference and the first harmonic (the only one that appears to add significant noise). There are several excellent signal processing texts available that go into significant detail on all types of filters, of course including FIR filters. They discuss them in much more detail (and with significantly greater authority) than I could, so I refer you to them. Many are listed in the references section at the end of each documentation page.

Mathieu NOE
on 29 Apr 2021

hello

dear all, an alternative using a 5th order IIR filter to remove the usual 50 Hz and harmonics mains hum

LD = load('ECG_sample.mat');

time_index = LD.time_index;

ecg_raw = time_index(:,2);

Fs = 5000;

Fn = Fs/2;

time = time_index(:,1)/Fs;

L = numel(ecg_raw);

N = 2^nextpow2(L);

FTs = fft(ecg_raw,N)/L;

Fv = linspace(0, 1, N/2+1)*Fn;

Iv = 1:numel(Fv);

figure

plot(Fv, abs(FTs(Iv))*2)

grid

title('Fourier Transform')

xlim([0 170])

%% two notch filters in series to remove 50 and 100 Hz hum

fc = 50; % notch freq for 50 Hz tone

Q = 7; % adjust Q factor for wider (low Q) / narrower (high Q) notch ; f = fc the filter has gain = 0

[B1,A1] = mynotchfilter(fc,Q, Fs);

fc = 2*50; % notch freq for 100 Hz tone

Q = 7; % adjust Q factor for wider (low Q) / narrower (high Q) notch ; f = fc the filter has gain = 0

[B2,A2] = mynotchfilter(fc,Q, Fs);

% two filter is series

B = conv(B1,B2);

A = conv(A1,A2);

ecg_filt = filtfilt(B, A, ecg_raw);

figure

freqz(B, A, 2^16, Fs)

figure

subplot(2,1,1)

plot(time, ecg_raw)

grid

title('Original Signal')

subplot(2,1,2)

plot(time, ecg_filt)

grid

title('Bandstop-Filtered Signal')

figure

subplot(2,1,1)

plot(time, ecg_raw)

grid

xlim([20 22])

title('Original Signal (Detail)')

subplot(2,1,2)

plot(time, ecg_filt)

grid

xlim([20 22])

title('Bandstop-Filtered Signal (Detail)')

function [B,A] = mynotchfilter(fc,Q, Fs)

%% notch filter section %%%%%%

% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)

% fc = 50; % notch freq

wc = 2*pi*fc;

% Q = 7; % 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;

B=[b0 b1 b2];

A=[a0 a1 a2];

end

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

Start Hunting!