Posts

The Removal of 50/60Hz powerline interference from delicate information rich ECG biomedical waveforms is a challenging task! The challenge is further complicated by adjusting for the effects of EMG, such as a patient limb/torso movement or even breathing. A traditional approach adopted by many is to use a 2nd order IIR notch filter:

\(\displaystyle H(z)=\frac{1-2cosw_oz^{-1}+z^{-2}}{1-2rcosw_oz^{-1}+r^2z^{-2}}\)

where, \(w_o=\frac{2\pi f_o}{fs}\) controls the centre frequency, \(f_o\) of the notch, and \(r=1-\frac{\pi BW}{fs}\) controls the bandwidth (-3dB point) of the notch.

What’s the challenge?

As seen above, \(H(z) \) is simple to implement, but the difficulty lies in finding an optimal value of \(r\), as a desirable sharp notch means that the poles are close to unit circle (see right).

In the presence of stationary interference, e.g. the patient is absolutely still and effects of breathing on the sensor data are minimal this may not be a problem.

However, when considering the effects of EMG on the captured waveform (a much more realistic situation), the IIR filter’s feedback (poles) causes ringing on the filtered waveform, as illustrated below:

Contaminated ECG with non-stationary 50Hz powerline interference (IIR filtering)

As seen above, although a majority of the 50Hz powerline interference has been removed, there is still significant ringing around the main peaks (filtered output shown in red). This ringing is undesirable for many biomedical applications, as vital cardiac information such as the ST segment cannot be clearly analysed.

The frequency reponse of the IIR used to filter the above ECG data is shown below.

IIR notch filter frequency response

Analysing the plot it can be seen that the filter’s group delay (or average delay) is non-linear but almost zero in the passbands, which means no distortion. The group delay at 50Hz rises to 15 samples, which is the source of the ringing – where the closer to poles are to unit circle the greater the group delay.

ASN FilterScript offers designers the notch() function, which is a direct implemention of H(z), as shown below:

ClearH1;  // clear primary filter from cascade
ShowH2DM;   // show DM on chart

interface BW={0.1,10,.1,1};

Main()

F=50;
Hd=notch(F,BW,"symbolic");
Num = getnum(Hd); // define numerator coefficients
Den = getden(Hd); // define denominator coefficients
Gain = getgain(Hd); // define gain

Savitzky-Golay FIR filters

A solution to the aforementioned mentioned ringing as well as noise reduction can be achieved by virtue of a Savitzky-Golay lowpass smoothing filter. These filters are FIR filters, and thus have no feedback coefficients and no ringing!

Savitzky-Golay (polynomial) smoothing filters or least-squares smoothing filters are generalizations of the FIR average filter that can better preserve the high-frequency content of the desired signal, at the expense of not removing as much noise as an FIR average. The particular formulation of Savitzky-Golay filters preserves various moment orders better than other smoothing methods, which tend to preserve peak widths and heights better than Savitzky-Golay. As such, Savitzky-Golay filters are very suitable for biomedical data, such as ECG datasets.

Eliminating the 50Hz powerline component

Designing an 18th order Savitzky-Golay filter with a 4th order polynomial fit (see the example code below), we obtain an FIR filter with a zero distribution as shown on the right. However, as we wish to eliminate the 50Hz component completely, the tool’s P-Z editor can be used to nudge a zero pair (shown in green) to exactly 50Hz.

The resulting frequency response is shown below, where it can be seen that there is notch at exactly 50Hz, and the group delay of 9 samples (shown in purple) is constant across the frequency band.

FIR  Savitzky-Golay filter frequency response

Passing the tainted ECG dataset through our tweaked Savitzky-Golay filter, and adjusting for the group delay we obtain:

Contaminated ECG with non-stationary 50Hz powerline interference (FIR filtering)

As seen, there are no signs of ringing and the ST segments are now clearly visible for analysis. Notice also how the filter (shown in red) has reduced the measurement noise, emphasising the practicality of Savitzky-Golay filter’s for biomedical signal processing.

A Savitzky-Golay may be designed and optimised in ASN FilterScript via the savgolay() function, as follows:

ClearH1;  // clear primary filter from cascade

interface L = {2, 50,2,24};
interface P = {2, 10,1,4};

Main()

Hd=savgolay(L,P,"numeric");  // Design Savitzky-Golay lowpass
Num=getnum(Hd);
Den={1};
Gain=getgain(Hd);

Deployment

This filter may now be deployed to variety of domains via the tool’s automatic code generator, enabling rapid deployment in Matlab, Python and embedded Arm Cortex-M devices.

 

 

Download demo now

 

Licencing information

 

Comb filters have found use as powerline (50/60Hz) harmonic cancellation filters in audio applications, and form the basis of so called CIC (cascaded integrator–comb) filters used for anti-aliasing in decimation (sample rate reduction), and anti-imaging in interpolation (sample rate increase) applications.

The frequency response of a comb filter consists of a series of regularly-spaced troughs, giving the appearance of a comb. As seen in the plot below, the spacing of each trough appears at either odd or even harmonics of the desired fundamental frequency.

Frequency response of a typical FIR comb filter (odd harmonics cancellation):
\(f_s=500Hz\),  \(f_c=25Hz\), \(L=10\) and \(\alpha=1\)

An FIR comb filter can be described by the following transfer function:

\(H(z)=1+\alpha z^{-L}\)
\(\Rightarrow Y(z)=X(z)\left[1+\alpha z^{-L}\right]\)

Clearly, the comb filter is simply a weighted delayed replica of itself, specifiied by \(L\). Taking inverse z-transforms, we obtain the difference equation needed for implementation,

\(y(n)=x(n)+\alpha x(n-L)\)

where, \(\alpha\) is used to set the Q (bandwidth) of the notch and may be either positive or negative depending on what type of frequency response is required. In order to elaborate on this, negative values of \(\alpha\) have their first trough at DC and their second trough at the fundamental frequency. Clearly this type of comb filter can be used to remove any DC components from a measured waveform if so required. All subsequent troughs appear at even harmonics up to and including the Nyquist frequency.

Positive values of \(\alpha\) on the other hand, only have troughs at the fundamental and odd harmonic frequencies, and as such cannot be used to remove any DC components.

Application to powerline interference cancellation

The affectivity of the comb filter is dependent on the sampling frequency, \(f_s\), as \(L\) is limited to integer values only. Also, a relationship between \(f_s\), as \(L\) and will be dependent on the sign of \(\alpha\). Thus, for the purposes of the mains cancellation application considered in this discussion, only positive values of will be considered, as we need only cancel odd harmonics.

A simple relationship for determining  \(L\) can be summarized for positive values of \(\alpha\) as follows:

\(L=ceil\left( \large{\frac{f_s}{2f_c}}\right)\)

where, \(f_c\) is the desired centre point of the fundamental notch frequency. Based on this expression, we can re-calculate the sampling frequency, such that \(f_c\) is a true multiple of \(f_s\)

\(f_{snew}=2f_c L\)

Example

For the example considered herein, i.e. \(f_s=500Hz\) and \(f_c=25Hz\), we obtain \(L=10\). However, if \(f_c=60Hz\), we would need \(L=5\), and a new sampling rate of \(600Hz\) respectively, although it’s interesting to note that \(f_s=480Hz\) for \(L=4\) would also suffice.

Implementation

An FIR comb filter may be implemented in ASN FilterScript as follows:

ClearH1;  // clear primary filter from cascade
interface L = {4,20,1,5}; // delay
interface alpha = {-1,1,0.01,1};

Main()
Num = {1,zeros(L-1),alpha}; // numerator coefficients
Den = {1};  
Gain = 1/sum(abs(Num));

 

 

Download demo now

 

Licencing information