Multivariate denoising using wavelets and principal component analysis

Abstract

A multivariate extension of the well known wavelet denoising procedure widely examined for scalar valued signals, is proposed. It combines a straightforward multivariate generalization of a classical one and principal component analysis. This new procedure exhibits promising behavior on classical bench signals and the associated estimator is found to be near minimax in the one-dimensional sense, for Besov balls. The method is finally illustrated by an application to multichannel neural recordings.

Introduction

On one hand, denoising algorithms based on wavelet decompositions are a popular method for one-dimensional statistical signal extraction and filtering. On the other hand, principal component analysis (PCA) is among the most notorious data-analysis tools designed to simplify multidimensional data by tracking new factors supposed to capture the main features.

This paper proposes a multivariate extension of wavelet denoising procedures, combining a straightforward multivariate generalization of the classical one for scalar valued signals and principal component analysis. This proposal takes place among the various recent approaches combining wavelet strategies and data analytic tools to cope with the problem of feature extraction in regression models. Numerous applied situations strongly motivate this interest. Let us mention some of them together with some references focusing on a wavelet-data analysis approach: spectral calibration problems (Vannucci et al., 2003), multivariate statistical process control (see Bakshi, 1998; Teppola and Minkkinen, 2000), blind source separation (Roberts et al., 2004), functional magnetic resonance imaging (fMRI) analysis (Meyer and Chinrungrueng, 2003), spike detection and sorting (Oweiss and Anderson, 2001).

This paper focuses on multivariate wavelet denoising and deals with regression models of the form $X(t)=F(t)+\epsilon (t)$, where the observation X is p-dimensional, F is deterministic and is the signal to be recovered and $\epsilon$ is a spatially correlated noise. This kind of model is well suited for situations for which such an additive spatially correlated noise is realistic. For example, a longitudinal study on p subjects, the analysis of a part a fMRI region (involving p voxels) or the noise reduction in multichannel neural recordings (using p channels). Let us be a little bit more precise on this last example which is chosen to illustrate the behavior of the proposed procedure on a real world data set, at the end of the paper. Following Oweiss and Anderson (2001), extra-cellular neural recordings can be modeled as an invariant deterministic signal and an additive noise which obscures neural discharges from cells of interest. This noise contains a component exhibiting spatial correlation coming from background activity caused by neural cells.

To close this introduction, let us recall some facts about classical univariate wavelet denoising dealing both with signal processing and functional estimation in statistics and which is of interest in various applied fields. Valuable references are the books (Mallat, 1998, Percival and Walden, 2000, Vidakovic, 1999) and the survey paper (Antoniadis, 1997). For basics on wavelets, we refer the reader to Mallat (1998) or Misiti et al. (2003) for example.

The simplest considered model is of the following form:

$$X(t)=f(t)+\epsilon (t),t=1,\dots ,n\text{,}$$

where $(X(t){)}_{1⩽t⩽n}$ is observed, $(\epsilon (t){)}_{1⩽t⩽n}$ is a centered Gaussian white noise of unknown variance ${\sigma }^2$ and f is an unknown function to be recovered through the observations.

For a given orthogonal wavelet basis denoted by $\left({\left({\phi }_{J,k}\right)}_{k\in Z},{\left({\psi }_{j,k}\right)}_{1⩽j⩽J,k\in Z}\right)$ where $\psi$ is a wavelet, $\phi$ the associated scaling function, J a suitably chosen decomposition level and where $g_{j,k}(x)=2^{-j/2}g\left(2^{-j}x-k\right)$, wavelet denoising proceeds in three steps:

• Step 1: Compute the wavelet decomposition of the observed signal up to level J;

• Step 2: Threshold conveniently the wavelet detail coefficients;

• Step 3: Reconstruct a denoised version of the original signal, from the thresholded detail coefficients and the approximation coefficients, using the inverse wavelet transform.

Various strategies are available (see the survey paper Antoniadis et al., 2001) to perform this task and the asymptotic performance of the associated estimators is the minimax one up to a logarithmic factor, for large classes of functions simultaneously (let us mention that block thresholding, see Hall et al., 1999, could be used to remove the logarithmic factor). For simplicity and since only relative performance between the proposed multivariate procedures are of interest, we restrict our attention to the so-called universal threshold (introduced by Donoho and Johnstone, 1994) which is of the form $\widehat{\sigma }\sqrt{2\log (n)}$, where $\widehat{\sigma }$ is an estimator of $\sigma$ based on the detail coefficients at level 1 (the finest one). Such methods are effective because functions f belonging to various general classes are such that they admit a sparse wavelet representation (Kerkyacharian and Picard, 2000). So the energy of f is mainly concentrated in a few large wavelet coefficients which are adaptively selected by this procedure since the coefficients below the threshold are attributable to the additive noise.

This paper dedicated to a multivariate denoising procedure that takes into account the correlation structure of the noise, is organized in two main sections. Section 2 proposes a first denoising procedure which is a direct generalization of the one-dimensional strategy. The method is based on a change of basis followed by a classical one-dimensional soft-thresholding. This new procedure exhibits promising behavior on classical test signals and the associated estimator is found to be near minimax in the one-dimensional sense, for Besov balls. The change of basis is obtained from the diagonalization of a robust estimate of the noise covariance matrix given by the minimum covariance determinant estimator based on the matrix of finest details.

Section 3 first recalls the multiscale PCA proposed by Bakshi (1998) for statistical process control purposes. This scheme is discussed and then a second denoising procedure combining wavelets and PCA is proposed. The introduction of a PCA step try to take advantage of the deterministic relationships between the signals, leading to an additional denoising effect. It is then illustrated by some simulation examples and by an application to multichannel neural recordings.