Powerline noise elimination in biomedical signals via blind source separation and wavelet analysis

Empirical mode decomposition

Suppose we want to eliminate a chronic sinusoidal xHz noise from a signal y(t) in a single channel. We can employ empirical mode decomposition (EMD) to split the time series signal into narrow-band amplitude–frequency sections (Huang et al., 1998; Wang et al., 2014). These sections, called intrinsic mode functions (IMFs), are obtained by an iterative sifting procedure. Initially, the local maxima and minima of the signal are detected and are connected by cubic splines to form the upper and lower envelopes respectively. The mean of the envelopes is then subtracted from the original signal to give rise to the first IMF. Subsequently, the first IMF is subtracted from the original signal to result in the residue. This process is repeated k times with the residue obtained at the end of each iteration serving as the input for the next. Ultimately, k IMFs and one residue will be obtained. The signal y(t) can be reconstructed by a summation of the IMFs c_i(t) and the residue r(t): (1)${\displaystyle y\left(t\right)=\sum _{i=1}^k{c}_i\left(t\right)+r\left(t\right).}$

The method by which the local extrema are detected is described in the Appendix (Ensemble empirical mode decomposition: detecting the local extrema).

Ensemble empirical mode decomposition

The ensemble empirical mode decomposition algorithm (EEMD) is an extension of EMD which is aimed at making the extraction of IMFs a robust process (Wu & Huang, 2009). In this algorithm, Gaussian noise of zero mean and a specified standard deviation is added to the input signal before EMD is applied. For the extraction of each IMF, the addition of Gaussian noise to the input is done over a specified number of ensembles and the mean extract in all the ensembles of putative IMFs is selected as the true IMF. This is a truly noise-assisted method of blind source separation because Gaussian noise forces the EMD algorithm to consider all options when sifting. With a high enough number of ensembles for each EEMD iteration, it can be inferred by central limit theorem that the mean of the ensembles is representative of the most likely IMF; thus, the fittest survive. The input signal for each ensemble can be summarized as follows: (2)${\displaystyle G\sim \mathcal{N}\left(0,{\sigma }_e^2\right);{\sigma }_e^2={\left[nl\ast \sigma \left(y\left(t\right)\right)\right]}^2;e_{h_{in}}=y\left(t\right)+G}$ where nl is the desired inverse signal-to-noise ratio of the input signal and e_hin is the input signal for each ensemble h.

IMF selection via wavelet analysis

A wavelet φ_(a,b)(t) is square-integrable function that can be dilated (or constricted) along a and translated along b; it is written in the following form: (3)${\displaystyle {\phi }_{\left(a,b\right)}\left(t\right)=\frac{1}{\sqrt{a}}\left(\frac{t-b}{a}\right),\forall a\in {\Re }^{+}\wedge \forall b\in \Re .}$ For this framework, the real component of the consecutive Morlet wavelet was employed: (4)${\displaystyle {\phi }_{\left(f,j\right)}\left(t\right)=\left(\sqrt{1+e^{-f^2}-2e^{\frac{-3f^2}{4}}}\right)\left(\frac{1}{\sqrt{2\pi }}{e}^{2\pi if}{e}^{\frac{-{\left(t-j\right)}^2}{2}}\right),}$ where f and j represent the center frequency (scale) and translation respectively. In order to describe the frequency properties of each independent component, a projection of the Morlet wavelet unto c_i(t) = [c₁(t), …, c_k(t), r(t)]^T was used. This projection p(F, i) was by accomplished by finding the frequency f in the set of frequencies $F={\left\{f_i\right\}}_{f_i=20}^{80}$ that maximizes a pseudo-convolution between φ_(f,j)(t) and c_i(t) and a summation of the pseudo-convolution values at all temporal locations t: (5)${\displaystyle p\left(F,i\right)=\underset{f\in F}{\mathrm{argmax}}\sum _{\underset{k}{t}}\left[\left|\left|\sum _{\underset{k}{j}}{\phi }_{\left(f,j\right)}\left(t\right)\ast c_i\left(t\right)\right|-\frac{1}{t_k{f}_k}\sum _{\underset{k}{t}}\sum _{\underset{k}{f}}\left|\sum _{\underset{k}{j}}{\phi }_{\left(f,j\right)}\left(t\right)\ast c_i\left(t\right)\right|\right|\right],}$ where j_k, t_k and f_k are integers representing the k^th translation, time and frequency evaluated respectively. A summation of the pseudo-convolution values at all temporal locations exposes the overall frequency morphology that describes a particular IMF. As mentioned previously, although each consecutive IMF exist within a lower frequency range, there is a likelihood of obtaining IMFs with partially overlapping frequency spectra (Fig. 1B). Thus, by choosing IMFs with dominant frequencies within a specified range, there is an appreciable degree of assurance that the xHz signal is in either of the chosen IMFs. Essentially, with Eq. (5), if xHz − bHz ≤ p(F, i) ≤ xHz + bHz (with bHz serving as a threshold for IMF selection), then ∃i such that i is the index that defines c_i(t) as a putative xHz powerline noise. In the following subsections, any c_i(t) such that xHz − bHz ≤ p(F, i) ≤ xHz + bHz is denoted as n_q(t). This method is further elucidated in the Appendix (IMF selection).

Independent component analysis

Independent component analysis (ICA) is a data-driven approach to blind source separation of mixed input signals into their independent components. The premise of ICA is the assumption that each vectorial input signal n_q(t) in n = [n₁(t), …, n_k(t)]^T is an observed signal which is a linearly mixed representation of its intrinsic source component s_i(t) and components of statistically independent origin (Hyvärinen, Karhunen & Oja, 2004). This assumption infers the following: (6)${\displaystyle n=As,}$ where A is the mixing matrix and the source component s = [s₁(t), …, s_k(t)]. Inevitably, finding the un-mixing matrix R to obtain a matrix of output signals u which are estimates of s_i(t) is the essence of ICA: (7)${\displaystyle u=Rn=RAs.}$

In this implementation of ICA, R was determined by an approximate simultaneous diagonalization of the third and fourth order cumulant tensors. Independent component analysis: finding R of the Appendix provides an explanation of how R was computed.

Powerline noise recognition

The powerline noise in u is recognized by extending the method outlined in Eq. (5) to search for the u_i whose frequency properties resemble that of xHz powerline noise the most: (8)${\displaystyle d_{comp}\left(t\right)=\underset{u_i\left(t\right)}{\mathrm{argmin}}\left(\left|\left[\underset{f\in F}{\mathrm{argmax}}\sum _{\underset{k}{t}}\left(\left|\left|\sum _{\underset{k}{j}}{\phi }_{\left(f,j\right)}\left(t\right)\ast u_i\left(t\right)\right|-\frac{1}{t_k{f}_k}\sum _{\underset{k}{t}}\sum _{\underset{k}{f}}\left|\sum _{\underset{k}{j}}{\phi }_{\left(f,j\right)}\left(t\right)\ast u_i\left(t\right)\right|\right|\right)\right]-x\hspace{0.17em}\right|\right),}$ where x represents the expected frequency of the powerline noise (50 Hz or 60 Hz). In essence, the u_i(t) that satisfies Eq. (8) is the most appropriate approximation of powerline noise. The motivation behind this approach has been explicated in Powerline noise recognition of the Appendix.

Characteristics of the proposed approach

To demonstrate the viability of the process of powerline noise removal, artificial 60 Hz noise was added to local field potentials (SNR = 4.2816dB; amplitude = 700) and was subsequently reconstructed. A 200 ms time window (without overlap) was used to denoise 250 ms of the mentioned neural data. As shown in Figs. 2A and 2C, the procedure was able to recover a sufficient amount of the original signal with little Manhattan distances. The morphology of the Manhattan distance suggests that the procedure extracted the appropriate frequency, but obtained relatively erroneous amplitude information. As shown in Figs. 2B and 2D, the pseudo-convolution procedure for extraction of powerline noise chose the 60 Hz noise—which is the right frequency for the noise model in this context. Since the pseudo-convolution between the extracted alternating current noise peaks at about 60 Hz on all translations, it implies that the summation of the pseudo-convolution values on every translation should result in a function that also peaks at about 60 Hz. Thus, whichever independent component that has its summed pseudo-convolution peak closest to 60 Hz is indeed the desired powerline noise. Although the final round of blind-source separation in this algorithm requires prior knowledge about the statistical properties of the source data being approximated, it is plausible that this constraint is not a hindrance but, instead, facilitates the process. With the assumption of statistical independence, at least one of the resulting approximations of the source data is forced to look undeniably unique. Signals of this form are usually some fluctuations that are chronically present in the mixed data. For this reason, ICA effectively serves as a helper for the extraction and identification of the desired powerline noise to be removed. The Manhattan distance between the original signal and the denoised signal increased during the final 50 ms of the neural data. This suggests the approach does not provide good results with very small time windows. Essentially, the best approach will be to run the algorithm without a time window. In conclusion, it is most appropriate for offline signal analysis.

Learning the noise amplitude

In accordance with the fact that the problem of deconvolution is ill-posed, the solutions are not unique; this issue is tackled by learning the appropriate amplitude for the noise model. To exhibit the amplitude learning procedure, an EEG with natural 50 Hz noise was denoised using the proposed framework. As presented in Fig. 3C, the learning procedure aimed at minimizing the cost function has a unique solution. This indicates that there is some level of consistency in the process. In view of the fact that there is virtually no ground truth with natural AC noise, the power spectra of the original signal and the denoised signal were compared. The traces in Fig. 3A show that there is little difference between the original signal and the recovered signal. However, Fig. 3B reveals that this manifestation is primarily due to the low amplitude of the 50 Hz AC noise. Although the amplitude is very small, the algorithm was able to extinguish an appreciable amount of powerline noise. The same figure indicates that if the amplitude of the AC noise is small enough, it might go unnoticed after Fourier transformation. Further, in order to obtain an accurate view of how a biomedical signal behaves under various frequency specifications, it is better to employ the wavelet transform (which is a convolution between a wavelet and a signal) rather than the Fourier transform. This is due to the fact that the Fourier transform assumes the signal being transformed is of a sinusoidal origin (Bloomfield, 2004), which is not true in most cases. Although the short-term Fourier transform is used to avoid this pitfall, the wavelet transform (which handles frequencies in a logarithmic fashion) adapts better to highly variable signals (Daubechies, 1990). The suggested algorithm employed a transformation analytically similar to the wavelet transform. By the same token, the frequency properties revealed by a summation of the pseudo-convolution between the real component of the Morlet wavelet and the signals aided in reliably in extracting powerline noise; this, in turn, provided a solid model foundation for learning the most appropriate amplitude.

Neural signals

Sub-cortical local field potentials and extra-cranial EEGs were resampled to have equal sampling rate, concatenated and adulterated with artificial powerline noise. A band-stop 4th order Butterworth filter was used to filter out signals between 59.5 Hz and 60.5 Hz, and the proposed framework was also applied to the adulterated signal. The introduction of powerline noise introduced a slight amplitude decorrelation between the original signal and the adulterated signal. The aim of Fig. 4 was to characterize the level of re-correlation after removing the powerline noise using the proposed approach and the mentioned infinite impulse response filter (IIR filter). It was noted that although the frequency spectra indicated that both methods expelled powerline noise, their correlations differed minutely. The proposed approach recovered the signals with a slightly higher correlation than the named IIR filter (proposed = 1.0000, Butterworth filter = 0.9999).

In order to understand the subtle differences between the proposed method and the mentioned IIR filter, the Manhattan distance between the Fourier transformation of the original signal which had no powerline noise and the denoised signal using both approaches were evaluated (Fig. 5). It is shown that the difference between the power spectrum of the original signal and the one reconstructed using the proposed approach was smaller than that of the band-stop Butterworth filter. This was further validated via a Wilcoxon rank-sum test; the p-value obtained by comparing the original power spectrum with that of the one reconstructed using the proposed approach was 0.9903, while that of the band-stop Butterworth filter was 0.7627. In spite of the fact that there was no statistically significant difference between the reconstructed power spectra and that of the original, there was a statistically significant difference between the power spectrum of the signals reconstructed using the suggested algorithm and that of the band-stop IIR filter (p-value < 0.0001). On a grand scale, the power spectra of the denoised signals were highly correlated. This is an indication that there was a relatively similar response at the powerline noise frequency for neural signals.

Electrocardiogram signals

Electrocardiogram signals without powerline noise were corrupted with artificial 60 Hz powerline noise. Thereafter, they were reconstructed via band-stop Butterworth filtering and the proposed method. The results obtained were almost identical in the time series (Fig. 6A), however there was a substantial difference circa 60 Hz. Similar to the results obtained with neural signals, the power of the extract using the proposed method at 60 Hz was lower than that of the Butterworth filter. One of the many properties of ECG signals sought after is the waveform. In accordance with the knowledge that the ECG signal used was a concatenation of three data sets with two electrode recordings each, it is expected that at least six different clusters ECG waveforms can be detected. The waveforms were detected via simple thresholding and were partitioned via the k-means algorithm. Naturally, clustering such high dimensional data requires dimensionality reduction via principal component analysis or laplacian eigenmaps. This clustering process was done without dimensionality reduction in order to accurately compare the effect of the proposed approach and that of the Butterworth filter on the reconstructed signal. As shown in Fig. 6B, although the labels of the classes were different, their elements were not. This further proves that the relative temporal morphology of each reconstructed signal is preserved.