1
$\begingroup$

Currently, I'm working on removing 50 Hz power line interference from an ECG signal. Before trying out notch filters, I decided to try out a simple lowpass filter with a cutoff less than 50 Hz. Here's the MATLAB code I used to make the filter:

Fs=500; %Sampling rate in Hz
Ast = 120; %Stopband attenuation
Ap = 1; %Passband ripple
Fp = 30; %Passband end
Fst = 45; %Stopband beginning
Hd=fdesign.lowpass('Fp,Fst,Ap,Ast',Fp,Fst,Ap,Ast,Fs);
d=design(Hd,'butter'); % Design a butterworth filter with the given characteristics
fvtool(d);

The magnitude response of the filter is: enter image description here

As you can see, the filter has around a 180 dB attenuation at 50 Hz and even more at higher frequencies.

Now, I run this filter on my data. Here's the original data, in the time and frequency domain:

enter image description here

And here's the data after applying the filter:

enter image description here

As you can see, the attenuation at 50 Hz is nowhere near the 180 dB that the filter advertised. For a better look, here's the plot of the magnitude response computed by taking the ratio of the FFTs of the signal before and after filtering:

The FFT ratio is in image 4.

Clearly, the attenuation is nowhere near the level it is supposed to be at. Is this the way it's supposed to be or am I doing something wrong in the implementation?

Thanks in advance for all the help!

Edit : Here's an enlarged view of a single beat of the ECG after applying the filter. Note the 50Hz noise right after the tall peak.

enter image description here

$\endgroup$
5
  • 1
    $\begingroup$ What command or commands are you using to plot the frequencies? I usually use something along the lines of "plot(20*log10(abs(fft(data))))". If you are doing 10*log10 the dB scale will be off because usually you want to get the power measurement, not the voltage measurement. $\endgroup$
    – Jim Clay
    Jun 7, 2014 at 15:23
  • $\begingroup$ For plotting the magnitude response of the filter, I'm using fvtool. I assume that this is equivalent to "plot(20*log10(abs(fft(data))))". For plotting the spectra of the signal before and after filtering, I'm using "semilogy(abs(fftshift(fft(data))))". $\endgroup$
    – Azura
    Jun 9, 2014 at 4:45
  • $\begingroup$ "semilogy" is doing 10*log10, so the frequency response of the before/after is off. Change it to "plot(20*log10" as shown in my previous comment. $\endgroup$
    – Jim Clay
    Jun 9, 2014 at 13:49
  • $\begingroup$ Assuming the x axis is in units of seconds then the ringing after the peak does not look like 50 Hz to me. I count 4 full cycles from about 3.82 s to 3.97 s. $\frac{1}{(3.97 - 3.82)/4} = 26.7 Hz$. $\endgroup$
    – Jim Clay
    Jun 9, 2014 at 15:45
  • $\begingroup$ Okay, it looks like the filter is working as it is supposed to. As for the ringing, I'd need to do some further research to figure out how to stop it. Thanks! $\endgroup$
    – Azura
    Jun 11, 2014 at 4:24

1 Answer 1

1
$\begingroup$

(Disclaimer: I don't have a Matlab license for trying your example.)

First, your ECG signal is quite short, and you don't seem to be applying any windowing before the FFT (which is equivalent to applying a rectangular window). So your final spectrum will be the circular convolution of your signal with the window transform, and hence the tails above 40Hz.

Second, you are designing for hundreds of dB attenuation. This yields filters of very high order, which are very sensitive to numerical precision errors during calculations. In real life you think in tens of dB, not hundreds.

In general, the best approach for cleaning ECG signals is:

  • Notch filter at 50 Hz (a band-pass of order 6 to 8 should suffice).
  • Low-pass filter with cutoff around 100 Hz (notice where you stop having ECG harmonics and start having only acquisition noise).

Edit for the downvoter: The filtered signal is OK, but using FFT without care gives a misleading result. At least you should apply a smooth window function to the signal (in this case, the already filtered signal) in order to get a more faithful spectrum representation. Try with a Hamming window: plot(abs(fft(y .* hamming(length(y))))), where y is the filtered signal.

$\endgroup$
5
  • $\begingroup$ I didn't -1 you but I don't see any indication that he is doing the filtering in the frequency domain. $\endgroup$
    – Jim Clay
    Jun 7, 2014 at 15:21
  • $\begingroup$ The main problem is that the OP is using the FFT of the resulting signal as the spectrum. FFT is not spectrum; FFT is spectrum with artifacts (as I explained in my answer: mainly circular convolution with window function). OP believes the artifacts are bad from bad filtering; filtering is OK, using FFT carelessly is the problem. $\endgroup$
    – Juancho
    Jun 7, 2014 at 15:33
  • $\begingroup$ I see now. Yes, that's a good point. In other words, the measurement is the problem, not the filtering. $\endgroup$
    – Jim Clay
    Jun 7, 2014 at 17:06
  • $\begingroup$ OP here. If I understand your response correctly, you're saying that the filter is working as intended but the problem is with using the FFT to visualize the spectrum. However, how does this explain my ECG signal still having noise in the time domain? Notice the oscillations in the ST segment (the part right after the sharp peak). I tried manually measuring the frequency of the oscillations by counting the number of peaks present per unit time, and it seems that the frequency is 50 Hz. If the filter is working, shouldn't this be gone? $\endgroup$
    – Azura
    Jun 9, 2014 at 4:51
  • $\begingroup$ Thanks for the zoomed-in plot. You've sure made a very high order filter. This results in what's called ringing, and is due to the sharp transition in the frequency domain; it's not noise. Using more reasonable specs for the filter should fix this. $\endgroup$
    – Juancho
    Jun 10, 2014 at 11:55

Your Answer

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

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