Sparsity-based correction of exponential artifacts

Highlights

• We proposed a method to suppress piecewise exponential transient.

• We proposed a method to suppress smooth protuberance transient.

• We give mathematical proof of quadratic majorizer of smoothed penalty functions.

Abstract

This paper describes an exponential transient excision algorithm (ETEA). In biomedical time series analysis, e.g., in vivo neural recording and electrocorticography (ECoG), some measurement artifacts take the form of piecewise exponential transients. The proposed method is formulated as an unconstrained convex optimization problem, regularized by smoothed ℓ₁-norm penalty function, which can be solved by majorization–minimization (MM) method. With a slight modification of the regularizer, ETEA can also suppress more irregular piecewise smooth artifacts, especially, ocular artifacts (OA) in electroencephalography (EEG) data. Examples of synthetic signal, EEG data, and ECoG data are presented to illustrate the proposed algorithms.

Introduction

This work is motivated by the problem of suppressing various types of artifacts in recordings of neural activity. In a recent study [25], typical artifacts in in vivo neural recordings are classified into four types (Type 0–3, see Section 2.2 and Figure 3 in [25]). This classification covers many artifacts in the scope of human brain activity recordings, e.g., electroencephalography (EEG) and electrocorticography (ECoG). In this paper, we consider the suppression of Type 0 and Type 1 artifacts. For the purpose of flexibility and generality, we redefine them in terms of morphological characteristics:

• Type 0: a smooth protuberance that can be modeled as $\widehat{x}\left(t\right)={\mathit{te}}^{-\alpha t}$, when $t\ge t_0$.

• Type 1: an abrupt jump followed by an exponential decay that can be modeled as $\widehat{x}\left(t\right)=e^{-\alpha t}$, when $t\ge t_0$.

Fig. 1 shows examples of the two types of artifacts.We do not consider the other two types in this work because our previous works have addressed efficient algorithms to remove such artifacts. For instance, low-pass filtering/total variation denoising ($\mathrm{LPF}/\mathrm{TVD}$) [61] suppresses Type 2 artifacts (Fig. 1(c)), and lowpass filtering/compound sparse denoising (LPF/CSD) [60], [61] can remove sparse and blocky spikes (Type 3 shown in Fig. 1(d)).

The approach proposed in this paper is based on an optimization problem intended to capture the primary morphological characteristics of the artifacts using sparsity-inducing regularization. To formulate the problem, we model the observed time series as

$$y(t)=f(t)+x(t)+w(t),$$

where f is a lowpass signal, x is a piecewise smooth transient signal (i.e., Type 0 or Type 1 artifacts), and w is stationary white Gaussian noise. More specifically, f is assumed to be restricted to a certain range of low frequencies. In other words, $\mathbf{H}\left(f\right)\approx \mathbf{0}$, where H is a high-pass filter. Note that in the signal model (1), conventional LTI filtering is not suitable to estimate either f or x from y, because component x, as a piecewise smooth signal comprised of transients, is not band limited.

In order to estimate the components, we combine LTI filtering with sparsity-based techniques. We formulate an optimization problem for both decomposition and denoising. A computationally efficient algorithm is derived to solve the optimization problem, based on the theory of majorization–minimization (MM) [15], [24], [36].

In addition, this paper specifies how to generate a smoothed penalty function and its majorizer from a non-smooth one, in order to overcome a numerical issue that arises when the penalty function is not differentiable.

Some recent works recover signals with transients by various algorithms. In [19], a slowly varying signal is modeled as a local polynomial and an optimization problem using Tikhonov regularization is formulated to capture it. In [47], the slowly varying trend is modeled as a higher-order sparse-derivative signal (e.g., the third-order derivative is sparse).

Instead of estimating the slowly varying component via regularization, the $\mathrm{LPF}/\mathrm{TVD}$ method [61] estimates a lowpass component by LTI filtering and a piecewise constant component by optimization. In this case, an optimization problem is formulated to estimate the piecewise constant component. The approach proposed here uses a similar technique to recover the lowpass component, but in contrast to $\mathrm{LPF}/\mathrm{TVD}$, it is more general—the regularization is more flexible with a tunable parameter, so that $\mathrm{LPF}/\mathrm{TVD}$ can be considered as a special case.

Another algorithm related to the approach taken in this paper is the transient artifact reduction algorithm (TARA) [60] which is utilized to suppress additive piecewise constant artifacts and spikes (similar to a hybrid of Type 2 and Type 3 artifacts). The approaches proposed in this work target different types of artifacts (Type 0 and Type 1) and applied in different applications.

The irregularity of Type 0 transients leads to a more complicated artifact removal problem, where the artifact are irregular fluctuations. A typical example in EEG is ocular artifacts (OA) caused by the blink and/or movement of eyes. To suppress OA, there are approaches based on empirical mode decomposition (EMD) [41], [42], [43], [69], and on independent component analysis (ICA) methods [1], [12], [20], [26], [37], [48]. The concept of spatial-frequency in acoustic analysis is also used to remove OA from multichannel signals [44], [64]. In this work, we present a new method to suppress ocular artifacts by proposing a specific model and using sparse optimization.

This paper adopts a regularizer inspired by the generalized 1-D total variation [27], wherein the derivative operator in conventional total variation regularizer is generalized to a recursive filter. The regularizers adopted in $\mathrm{ETEA}$ and second-order $\mathrm{ETEA}$ coincide with first-order and second-order cases of generalized 1-D total variation, respectively. Some differences to the problem discussed in [27] are as follows. Firstly, the signal model (1) allows a lowpass baseline as a component, hence, $\mathrm{ETEA}$ can be seen as a combination of conventional LTI filtering and generalized 1-D total variation. Secondly, we consider a formulation in terms of banded matrices, for computational efficiency. Thirdly, we give optimality conditions of the proposed problems, and use these conditions as a guide to set the regularization parameters.