Digital Filtering - A Fundamental Pre-Processing Step

Difficulty Level:

The acquired electrophysiological signals have always two intrinsic components. The signal we really want to acquire/study and noise, i.e. the acquisition component that is not relevant for the study purposes.

Noise can have different origins, such as in random events or due to voluntary/involuntary movements of the subject under analysis that affect the signal acquisition .

So, filtering is a fundamental step that needs to be applied to the signal, in order to ensure the maximisation of signal to noise ratio . Filtering can be achieved by hardware, having the analogical systems a great relevance, or by software using digital filters.

In this Jupyter Notebook it will be demonstrated how to digital filter the signal.

1 - Importation of the needed packages

# biosignalsnotebooks own package for loading and plotting the acquired data
import biosignalsnotebooks as bsnb

# Scientific package
from numpy import array, mean, average, linspace, where
from numpy.random import normal

# Hide warnings.
# https://groups.google.com/a/continuum.io/forum/#!topic/bokeh/h817iNS2twk
from IPython.display import HTML
HTML('''<script>
code_show_err=false; 
function code_toggle_err() {
 if (code_show_err){
 $('div.output_stderr').hide();
 } else {
 $('div.output_stderr').show();
 }
 code_show_err = !code_show_err
} 
$( document ).ready(code_toggle_err);
</script>
To toggle on/off output_stderr, click <a href="javascript:code_toggle_err()">here</a>.''')

2 - Load of acquired ECG data

# Load of data
data, header = bsnb.load_signal("ecg_4000_Hz", get_header=True)

In the following cell, some relevant information is stored inside variables. This relevant information includes the mac-address of the device, channel number and signal acquisition parameters such as resolution and sampling rate.

For a detailed explanation of how to access this info, the "Signal Loading - Working with File Header" Notebook should be consulted.

ch = "CH1" # Channel
sr = header["sampling rate"] # Sampling rate
resolution = header["resolution"] # Resolution (number of available bits)

3 - Generation of signal power spectrum by Fast Fourier Transform (FFT)

With this step is possible to observe the frequency composition of ECG signal.

3.1 - Store the desired physiological data (channel 1) int an individual variable

# Acquired data
signal = data[ch]

3.2 - Removal of continuous component from our signal (baseline shift through the subtraction of the average value)

This task ensures more stability of our filtering system.

# Baseline shift.
signal = array(signal) - mean(signal)

3.3 - Generation of the power spectrum

# Power spectrum
freq_axis_1, power_spect_1 = bsnb.plotfft(signal, sr)

4 - The informational content of ECG signal is typically contained below the 40 Hz frequency component

With the next representation we can conclude that exist some unwanted information out of this frequency band.

fig_1 = bsnb.plot_informational_band(freq_axis_1, power_spect_1, signal, sr, band_begin=0.5, band_end=40, legend="ECG Power Spectrum", 
                                     x_lim=[0, 100], y_lim=[0, 5e6], show_plot=True)

5 - Application of a low-pass filter in order to be excluded the unwanted information above the 40 Hz frequency component

Some low-frequency noise can be present at [0, 0.5] Hz frequency band. To exclude it we can follow an identical procedure, but, instead of applying a low-pass filter, it should be used a band-pass filter for the frequencies inside [0.5, 40] Hz.
For now, we focused on the more problematic type of noise, i.e., the high frequency noise

# Digital lowpass filtering with a cutoff frequency f of 40 Hz
filter_signal_1 = bsnb.lowpass(signal, f=40, order=1, fs=sr)

# Power spectrum
freq_axis_2, power_spect_2 = bsnb.plotfft(filter_signal_1, sr)

6 - Comparison of the power spectrum of original and filtered signal

bsnb.plot_before_after_filter(signal, sr, band_begin=0.5, band_end=40, x_lim=[0, 100], y_lim=[0, 5e6], show_plot=True)

[[Figure(id='1180', ...), Figure(id='1282', ...)]]

In this first filtering attempt we used a first order filter (input argument order=1). It can be seen, in the previous figure, that some unwanted information have been removed, unfortunately no filter has an ideal behavior, so despite we specify a high cutoff frequency of 40 Hz, some information above this threshold is maintained after filtering.

The good news are that components greater than 80 Hz are almost completely removed.

The filter performance can be improved by increasing the filter order, because the higher the filter order is, more quickly the transition between the pass and stop band will be. The transition band will be smaller because of a higher attenuation rate (-20 x order dB/decade).

However, the filter order must be chosen with precaution in order to avoid system instability. Magnitude Bode plots are very useful to check the filter response, as can be seen in the figure below, taking into consideration the following Mathematical formulation.

\begin{equation} G = -20\times\log_{10}\Bigg(\sqrt{1 + \bigg(\frac{f}{f_c}\bigg)^{2.n}}\Bigg) \end{equation}

	$G$ - Gain factor (negative values reveal an attenuation)	$n$ - Filter order (integer)
	$f_{c}$ - Cutoff frequency of the filter (40 Hz, for the current implementation)	$f$ - Independent variable (input frequency to be filtered)

bsnb.plot_low_pass_filter_response()