- Email: [email protected]
Smoothly adjustable denoising using a priori knowledge
ARTICLE IN PRESS
Signal Processing 86 (2006) 2233–2242 www.elsevier.com/locate/sigpro
Smoothly adjustable denoising using a priori knowledge Patrick Celkaa,,1, Elly Gyselsb a
School of Engineering, Griffith University, Gold Coast, Queensland Australia Control and Signal Processing Section, Swiss Center for Electronics and Microtechnology, Rue Jaquet Droz 1, 2007 Neuchaˆtel, Switzerland
b
Received 7 April 2005; received in revised form 30 September 2005 Available online 16 November 2005
Abstract Assuming a signal is composed of information and noise, this paper presents a generic approach to denoising by mapping the noisy signal using a priori information about the signal to be retrieved. The method is based on a subspace decomposition of both the a priori information at disposal and the noisy signal, followed by shrinkage of both subspace coefficients and smooth mapping of first onto the second space. Our method propose 3 different ways of building the a priori knowledge: A model-base, averaging-based and recursive-based. The proposed method is particularly suited for low signal to noise ratio. The denoising methods are validated on synthetic electrocardiogram (ECG) signals and further assessed on real life ECG and visual brain evoked potentials using the wavelet transform. r 2005 Elsevier B.V. All rights reserved. Keywords: Subspace denoising; Wavelet transform; Biosignal processing; Electrocardiogram; Evoked potential
1. Introduction For several decades, scientists and engineers have developed strategies to reduce the noise from signals. Well known techniques are Wiener filtering [1,2], spectral analysis [3], model-based analysis [3], principal component analysis [4], wavelets [5], independent component analysis [6], nonlinear techniques [7], local subspace projection [8], probabilistic approaches [9], and combinations of these methods [10–13]. Only some of them, such as Corresponding author. Tel.: +61 7 5552 8459;
fax: +61 7 5552 8065. E-mail address: p.celka@griffith.edu.au (P. Celka). 1 The work has partly been done at the Swiss Center for Electronics and Microtechnology, and was supported by the Swiss IM2 NCCR Program.
principal component analysis (local or global) and wavelets, are suited for non-stationary signals. Most of these denoising techniques use a particular subspace decomposition operator H ¼ fH k gM k¼1 that maps the signal þ noise space S into ~ k , for k ¼ 1; . . . ; M, a particular set of subspaces S where noise and signal components are expected to be as separated as possible in a given normed space. ~ k , a weighting operator Pk Then, in each subspace S reduces the noise components. We will call the entire space weighting operator P ¼ fPk gM k¼1 . The weighting operation is often performed by thresh~ k or by olding of the components in the space S Wiener filtering [12]. The weighted components now ~ c . Finally, by inverse transform span the spaces S k M H1 ¼ fH 1 g , k k¼1 the signal is reconstructed in the signal space Sc 2 S.
0165-1684/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2005.10.003
ARTICLE IN PRESS 2234
P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
Fig. 1. Subspace denoising in S which shows the action of the subspace operator H and the weighting operator P.
A measurement signal xðnÞ (n is the discrete time index) in S can be written as xðnÞ ¼ sðnÞ þ ZðnÞ
with n ¼ 0; . . . ; N 1,
(1)
where sðnÞ is the information signal component we want to retrieve, and ZðnÞ is the noise signal component. The denoising scheme is displayed in Fig. 1. The following signals are generated through the denoising procedure: ~ xðnÞ ¼ HðxðnÞÞ, ~ s~ðnÞ ¼ PðxðnÞÞ,
ð2Þ ð3Þ
s^ðnÞ ¼ H1 ð~sðnÞÞ,
ð4Þ
where s^ðnÞ is the estimated information signal. The choice of the subspace operator H is most often an open question and depends on the application, the underlying generating system linearity, the stationarity, the available computing and memory resources and the signal statistics. Generally, orthogonal filter banks, such as principal components filter banks [10] or discrete orthogonal wavelet transforms [14], are preferred because of their numerous advantages such as perfect reconstruction of the signals, and fast implementation. The use of the wavelet decomposition allows for a natural subband analysis of the signals under inspection and a good joint time–frequency decomposition [5,15]. This implies that denoising using wavelet decomposition will lead to the best performance for non-stationary signals. The choice of the wavelet basis depends largely on the time and frequency behaviour of the information we want to extract. Denoising using wavelets usually requires soft or hard thresholding on the wavelet coefficients [5,16,17]. Determining the best threshold value is a matter of looking more closely at the signals at hand and the purpose of the denoising. We will use the discrete time wavelet subspace operator H ¼ W in
the following, but our approach can be implemented with any subspace decomposition. In some cases, we are not totally blind for performing the denoising. We may have some information given a priori such as the signal or noise statistics in the time or frequency domains. We know that the best denoising filter for stationary time series with a priori knowledge of the noise and signal statistics is the Wiener filter. Even if we have only an approximation of the signal and noise statistics at disposal, and the signals are nonstationary, we may do better than just blind denoising. Our approach is based on a smooth mapping of the subspace component of the given a priori signal sa ðnÞ and that of the signal xðnÞ. The smooth mapping allows the user to choose different levels of denoising. This is of great interest in medical applications where the signal-to-noise ratio (SNR) can be very low, such as evoked potential (EP) studies [18], electrocardiogram (ECG) and photoplethysmograph measurements recorded with wearable sensors [19]. The paper is organised as follows: Section 2 formulates the smooth mapping approach. Section 3 discusses the construction of the a priori knowledge signals. In Section 4, results based on synthetic and experimental signals are presented. Section 5 discusses the results and Section 6 concludes the paper. 2. Smooth mapping Let us assume that we have at disposal an a priori signal sa ðnÞ which contains most of the time and frequency information of the signal sðnÞ, and a measurement signal xðnÞ. We transform both signals xðnÞ and sa ðnÞ through the subspace operator H ¼ W ~ xðnÞ ¼ W ðxðnÞÞ, a s~ ðnÞ ¼ W ðsa ðnÞÞ.
ð5Þ ð6Þ
~ The signals xðnÞ and s~a ðnÞ are the discrete wavelet transforms of the signals xðnÞ and sa ðnÞ, and the subspace operator W is the N N wavelet transform matrix [14]. Eq. (5) can be rewritten in a matrix form as x~ ¼ W ðxÞ, s~ a ¼ W ðsa Þ,
ð7Þ ð8Þ
where the bold variables stand for the column vectors of samples from the corresponding signals’
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
spaces. The coefficients’ amplitude and positions in the vector s~ a are very close to the one we should find in s~ , the wavelet transform of the signal sðnÞ. Quiroga [18] has proposed to manually locate the most significant coefficients in s~ . The positions ni of the selected coefficients s~ðni Þ are stored in a set I M ¼ fni gK i¼1 . Using this set of positions I M , he proposed to locate the same coefficients in x~ and subsequently zeroed the other ones. This corresponds to our previously defined weighting matrix operator P which can be formulated as P ¼ DiagðpðkÞÞ
for k ¼ 1; . . . ; N
(9)
with pðkÞ ¼ 1 for k 2 I M , ¼ 0 for keI M ,
ð10Þ ð11Þ
where DiagðÞ is a diagonal matrix. Inverse transforming Px~ results in the denoised vector s^ ¼ W
1
~ Px.
(12)
Note that this weighting matrix is different from a simple thresholding technique because of the user manual selection. On the one hand, manually selecting the positions I M results in a supervised denoising and leads to optimal solutions in terms of expert knowledge. On the other hand, it is a drawback because of the subjectivity of the selection, the potential human mistakes from the ‘getting bored’ effect, and thus the lack of reproducibility. This method cannot be applied when a fully automatic denoising procedure is required. Also, this approach assumes that the signal sa ðnÞ is very similar to sðnÞ in the time and frequency domain, which may not always be the case. We propose a flexible approach which enables the user to adjust the denoising to his needs by adjusting very few parameters. First we locate the most significant coefficients in s~a by a wavelet soft thresholding operator PWs S [12], resulting in the set of positions I S for non-zero coefficients. Second, we locate the coefficients in x~ which correspond to the positions in I S and we apply a FIR Gaussian filter 2 2 GðkÞ ¼ ek =s , for k ¼ 0; . . . ; K 1, which is summarised in the PWx S operator. The Gaussian filtering in (16) allows the user to account for a possible time shift of the signal s~ compared to s~a , and is performed using a forward–backward approach for cancelling its phase effect. This results in the following two-step nonlinear weightings, coefficient selection and filtering operations:
2235
First step: The signal s~a is softly thresholded using ~a , s~aS ¼ PWs S s
(13)
where the weighting operator is given by Ws PWs S ¼ DiagðpS ðkÞÞ
(14)
and the soft thresholding nonlinear function expressed by t Ws a pS ðkÞ ¼ Yðj~s ðkÞj tÞ 1 a , (15) j~s ðkÞj where Y is the Heaviside function, t is the threshold. Second step: The resulting non-zero coefficients’ positions in s~aS are stored in the set I S . The set I S is used for the selection of the coefficients in x~ S as finally expressed by the following equation: ~ x~ S ¼ PWx S x,
(16)
where the selection and filtering operator is given by Wx PWx S ¼ GðkÞc DiagðpS ðkÞÞ
(17)
with c the convolution operator, and pWx S ðkÞ ¼ 1 ¼0
for k 2 I S , for keI S .
ð18Þ
The selection operation expressed in (18) allows for the extraction of small information bearing signal components embedded in x. Note that this selection procedure is different from hard or soft thresholding, and is nonlinear by nature. Finally, we apply a soft mixing of the coefficients s~aS and x~ S x~ Smooth ¼ ax~ S þ ð1 aÞ~saS ,
(19)
where a is the mixing parameter. The value of s in the Gaussian filter depends on how confident we are in the phase information contained in the a priori signal sa ðnÞ: the more confident we are, the smaller the width. The value of a reflects our confidence in the amplitude of the a priori signal sa ðnÞ. Using a ¼ 0 in the third step implies that x~ Smooth ¼ s~ aS , which means that we are fully confident in our a priori knowledge. Using a ¼ 1 results in x~ Smooth ¼ x~ S , which uses any a priori knowledge and signifies that we have no confidence at all in it. While the phase information can be of great importance in some applications such as brain EP (see [18] and references therein). We will focus on the amplitude in this work. Determining a is the main issue in this approach and is left to the expertise of the user, which renders this method flexible and easy to use. Yet, in the result section we present some insights
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
2236
about the performance of the method for different values of a which could help the choice of the user. The estimated information signal is recovered as follows s^ ¼ W 1 x~ Smooth ¼W
1
ðaPWs S Wx
ð20Þ þ ð1
a aÞPWs S W s Þ.
ð21Þ
3. Building the a priori 3.1. Averaging approach In most applications, users have some a priori knowledge about the nature of the system or signals under study. For instance, if one wants to perform the denoising of EP in the brain, the user may then record a few such signals sai ðnÞ, for i ¼ 1; . . . ; R, prior to the denoising. The a priori signal saA ðnÞ can then be obtained by averaging these recordings P a s ðnÞ a . (22) sA ðnÞ ¼ i i R In this particular case, the true EP saT ðnÞ is not known, but approximated by the averaging, which still is a noisy and distorted version of the true EP. The validity of the averaging approach is based on the assumption that the true signal to be recovered saT ðnÞ is deterministic, time-invariant, with added zero-mean noise. This method is demonstrated in Section 4.2.2.
i being the time index of the pattern. Such a signal can be, for example, an ECG and its PQRST complex the pattern [20,21] (see Fig. 3). Second, the user is making repetitive measurements sa ðn; iÞ on the system, with i being the repetition index. The first case involves time repetitiveness and each pattern has to be detected prior to building the a priori. Note that the repetitiveness does not need to be periodical or quasi-periodical. The second case concerns an experimental repetitiveness and each separate measurement can be used to build the a priori. In both cases, the recursive a priori signal approach is performed by an autoregressive moving average procedure saR ðn; iÞ ¼ lsaR ðn; i 1Þ þ ð1 lÞsa ðn; iÞ.
(23)
The parameter 0:5plo1 controls the recursive averaging memory. This approach is assessed in Section 4.1.2. 4. Results The method presented in the previous section is tested on synthetic and real life signals. The wavelet subspace operator W is used. The mother wavelet used in these experiments, quadratic biorthogonal B-splines [22], has been chosen to match the signal at hand: a set of biphasic smoothly varying waveforms concatenated in a so-called PQRST complex (see Fig. 3). The Gaussian window width was fixed to s ¼ 1 throughout this section. A wavelet decomposition up to 5 scales was used.
3.2. Model approach When a reliable model of the system or signals under inspection is at disposal, the a priori signal is given by saM ðnÞ. Again, the model is a distorted version of the true information signal saT ðnÞ. One advantage of this situation is that the signal saM ðnÞ is noiseless and can be produced by a time-varying system. However, the user is bound to its model and its consistency for the data at hand. This approach can also be used in combination with the averaging and the recursive methods. This approach is assessed in Section 4.1.1.
4.1. Synthetic signals We compare the recursive and model-based approaches with synthetic ECG data. The synthetic ECG is produced using the nonlinear dynamical model proposed by McSharry et al. [21] and provides us with saM ðnÞ. The parameters used in our simulation are those reported in [21], with a mean heart rate of 70 beats per minute, a lowfrequency to high-frequency band power ratio of LF=HF ¼ 0:5, and a standard deviation of the heart rate of 5 beats per minute. The sampling frequency is 512 Hz.
3.3. Recursive approach The last way to build the a priori signal is to use a recursive approach. There are two ways to implement the recursion. First the signal to be denoised is a repetition of a slowly time-varying pattern sa ðn; iÞ,
4.1.1. Model-based denoising The synthetic ECG signal is used then as a model of an ECG, to which noise has been added to simulate real ECG measurements. The added noise is built using real measurement noise from an ECG
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
signal recorded using electroactive polymer-based fabric electrodes similar to those used in [23]. Subject movement produces highly noisy data. Control of the SNR is thus made possible. Noisy ECG data can then be generated at will using the nonlinear dynamical model and Fourier-based amplitude corrected surrogate noise time series [24]. Our results are quantified by the Normalised Mean Square Error expressed in percent:
2237
Noisy ECG Original ECG Denoised ECG Wavelet
R P
Q
S T
1
(24)
where E½ stands for the mathematical expectation, and signals are assumed to have zero mean. Fig. 2 shows the averaged NMSE in function of the SNR and a when using the model-based denoising approach. One hundred Monte Carlo runs have been generated, for each value of SNR and a, to compute the averaged NMSE. It should be noted that the averaging a priori method of Section 3.1 would have had produced the same qualitative results. We can observe from Fig. 2 that for low values of the NMSE, when the SNR is decreasing, a should be increased to reach a given NMSE. As discussed previously, a is the confidence level of the user on the a priori signal. Thus, truthful processing of the data should use some boundaries of a, and thus will automatically constrain the minimum
α
0.75
E½ðsaM ðnÞ s^ðnÞÞ2 NMSE ¼ 100 , E½ðsaM ðnÞÞ2 ðnÞ þ E½^s2 ðnÞ
0.5
0.25
0
0.1
0.2
0.3
0.4 0.5 Time [s]
0.6
0.7
0.8
0.9
Fig. 3. A series of model-generated ECG signals showing the complex P-, Q-, R-, S- and T-waves, with an SNR ¼ 3 dB (solid thin), together with the noiseless ECG (dotted-dashed) and denoised signals (solid bold), for 5 values of a.
reachable NMSE. For a fixed value of a, we can observe that the NMSE curve decreases faster for high SNR than for low SNR, and that the slope of the curve is inversely proportional to a. Fig. 3 shows a series of synthetic noisy ECG ðSNR ¼ 3 dBÞ, for a ¼ 0; 0:25; 0:5; 0:75; 1, together with the model ECG signal saM ðnÞ used for denoising, and the denoised signal. We can observe how the different ECG waves (namely P-, Q-, R-, S-, and T-waves) are recovered during the denoising process as we gradually increase a from 0 to 1. Special attention should be drawn towards the extraction of the P-, Q- and S-waves which are small compared to the R- and T-waves.
Fig. 2. This figure shows 100 Monte Carlo NMSE averaged for different values of SNR and a. Nine different iso-NMSE are shown.
4.1.2. Recursive-based denoising For the recursive approach of Section 3.3, an ECG of 60 s has been generated with the same parameters as previously described with an SNR ¼ 3 dB. Each ith peak R-wave has been detected using a method similar to the one proposed by Afonso [25]. We have used a 5 s sliding window with half overlap for
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
2238
α=0.9
α=0.9
λ=1
λ=0.95
λ=0.95
λ=0.9
λ=0.9 NMSE
ECG
λ=1
λ=0.85
λ=0.85 λ=0.8
λ=0.8
23
23.5
24
24.5 Time [s]
25
0
25.5
Fig. 4. A series of denoised model-generated ECG signals ðSNR ¼ 3 dBÞ showing P-, Q-, R-, S- and T-waves for 5 values of l, together with the noisy ECG signal.
4.2.1. Electrocardiogram denoising ECG electroactive polymer-based fabric electrodes measurements similar to those presented in [23] have been performed for assessing the proposed denoising method. The sampling frequency is 250 Hz. We have used a 5 s sliding window with half overlap for processing the signals. An ith peak R-wave symmetrical window of 500 ms was applied yielding the signal sa ðn; iÞ for n ¼ 1; . . . ; 125 and the number of QRS complexes was 97. During the recording, the subject was asked to perform different slow movements such as squeezing shoulders, trunk rotation and deep breathing. The ECG measurements display various artefacts linked to these movements as shown in Fig. 5. Figs. 5 and 6 show different results of the recursive approach using different values of l and a ¼ 0:9. In Fig. 5, the NMSE values have been computed on each detected PQRST complex using a 20-averaged PQRST complex saA ðnÞ E½ðsaA ðnÞÞ2 ðnÞ þ E½^s2 ðnÞ
.
(25)
30 40 Time [s]
50
60
70
λ=1
α=0.9
λ=0.95 λ=0.9 λ=0.85
ECG
4.2. Measured signals
E½ðsaA ðnÞ s^ðnÞÞ2
20
Fig. 5. NMSE computed for each of the detected PQRST complex for 5 different values of l together with the measured ECG signal.
processing the signals. An ith peak R-wave symmetrical window of 500 ms was applied yielding the signal sa ðn; iÞ for n ¼ 1; . . . ; 256 and i ¼ 1; . . . ; 75. Fig. 4 shows 2.5 s of the denoised ECG signals for five values of l, together with the noisy ECG.
NMSE ¼ 100
10
32
34
36
38
40 42 Time [s]
44
46
48
50
Fig. 6. A series of denoised measured ECG signals for 5 values of l, together with the measured ECG signal.
The detection method is based on the work of Afonso [25]. The NMSE displays fluctuations due to changes in the PQRST pattern. These changes can be functionally related to the heart, or to artefacts due to the subject’s movements. Note that as l approach 1, the fluctuations approach zero (see (23)). These fluctuations in the NMSE could be used for detection of abnormal PQRST complexes when the signals are not perturbed by artefacts. Motionrelated artefact perturbations can be detected using other type of sensors, e.g. accelerometers [26,27].
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
Fig. 6 clearly shows how the l parameter influences the denoising results. A value of l ¼ 1 should be avoided for studies related to fine analysis of pattern changes because of the infinite memory effect of the recursive averaging method. Any value of lo0:5 should also be avoided due to the too short-range memory effect of the recursive averaging and poor denoising results.
25
Original 400 trials Averaged Smooth Wavelet cleaned
20
α=0.3
15 10 EEG channel #19
4.2.2. Brain evoked potentials denoising In this section, we present the results of denoising on visual brain evoked potential (VEP) which is a sub-class of EP. A visual stimulus was presented to a subject from which an electroencephalogram (EEG) was recorded. The EEG was recorded using the standard 10–20 system electrode placement and we have used 47 channels. Careful visual inspection of the signals was performed to discard any recording with a bad skin contact which produced heavily artifactual data. The signals were sampled at 500 Hz and lasted for 800 ms with 100 ms of pre-stimulus. R ¼ 400 trials were recorded and further averaged to provide the signal saA ðnÞ (see Section 3.1). Neuro-scientists and -clinicians often undertake a large collection of visual stimulus triggered VEP response signals, and then perform a statistical averaging of the stimulus locked VEP. In recent years, it has been pointed out that the statistical averaging can in fact deteriorate the single trial informational signal. Moreover, due to the natural non-stationarity of the signals and the experimental procedure, the statistical average cannot be applied. For this reason, methods for using single-trial EP have been developed. Due to the low SNR of the single trials, noise reduction methods have been developed with more or less success [18,28–30]. The commonality between these 4 techniques is the use of thresholding on wavelet coefficients to perform the actual denoising, and sometimes the manually selected set of good wavelet coefficients. In our approach we use an automatic coefficient selection ~ with or without in the a priori signal dual space S, thresholding its coefficients, and two adjustable parameters, a and s, which give the confidence in the a priori knowledge. The major difficulty is to differentiate the noise from the information signal in low SNR brain signals. We report here the results of our denoising procedure when using the averaging a priori method. The average EVP saA ðnÞ is still the ‘gold standard’ of EP analysis by neurologists and for this
2239
5 0 -5 -10 -15 -20 -0.1
0
0.1
0.2
0.3 0.4 Time [s]
0.5
0.6
Fig. 7. One VEP as measured together with a 400-trials average and the denoised VEP for a ¼ 0:3.
reason, our smooth mapping approach is perfectly suited to this type of application. Figs. 7–9 show one EVP as measured, together with the 400-trial averaged signal saA ðnÞ, and the denoised signal for a ¼ 0:3; 0:6; 0:9. To perform the denoising, we have used saA ðnÞ as a priori signal as proposed in Section 3.1. A bi-phasic wave pattern at about 150 ms and a large positive deflection at about 300 ms from the onset of the stimulus can be clearly distinguished on the denoised VEP. Figs. 7–9 illustrate the smoothing concept: while a is approaching 1, the denoised signal is closer to the a priori, while when it is decreased toward 0, the denoised signal is closer to the single trial. Note the significant difference between the single trial and the average signal indicating the large inter-trial variability. The user can choose a number of trials Ro400 to build the average a priori signal. However, due to the large variability of single trials SNR, we recommend to select the best VEP for building the a priori.
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
2240
25
25
Original 400 trials Averaged Smooth Wavelet cleaned
20
Original 400 trials Averaged Smooth Wavelet cleaned
20
α=0.9
15
15
10
10 EEG channel #19
EEG channel #19
α=0.6
5 0
5 0
-5
-5
-10
-10
-15
-15
-20
-20
-0.1
0
0.1
0.2
0.3 0.4 Time [s]
0.5
0.6
-0.1
0
0.1
0.2
0.3 0.4 Time [s]
0.5
0.6
Fig. 8. One VEP as measured together with a 400-trials average and the denoised VEP for a ¼ 0:6.
Fig. 9. One VEP as measured together with a 400-trials average and the denoised VEP for a ¼ 0:9.
5. Discussion
measured denoised signal in the subspace domain enable the fine tuning readily. While fine tuning is the core of the method, we have shown how the global measure of NMSE can guide the user towards well-chosen initial parameters. The main parameters are: (1) the mixing parameter a, (2) the Gaussian width of the subspace filter s and (3) if the recursive methodology is used, the forgetting autoregressive smoothing factor l. The s parameter has not been studied thoroughly in this work and will be addressed in a further study. The influence of s is particularly important when timing (phase information) is involved, e.g. evoked brain potentials and heart rate variability analysis. Applications to synthetic and measured ECG signals have shown the potential of the method, and its efficiency for the extraction of the P-, Q- and S-waves which are usually very small compared to the R- and T-waves. Experimental illustrations of the method using the average a priori on visual EP has proved to be efficient in extracting the information contained in the biosignals. Further neurophysiological and cardiovascular validation of the
The approach we have developed for denoising signals is flexible due to the restricted number of free parameters and thus allows the user to adapt the method to the purpose of her/his study. The proposed methodology stands between a pure automatic denoising procedure, where the user has no control of the algorithm, and a highly tuneable approach, in which the user can effectively fine-tune the algorithm. However, the latter can result in a subjective choice of the tuning parameters and also require a strong knowledge of signal processing. Our method is based on the availability or construction of a priori signals, and we have shown how to construct such a priori knowledge. The method can be used with any subspace approach with some slight modifications. The use of this a priori knowledge makes the method robust to low SNR and at the same time allows the extraction of small information patterns hidden in the noise, similarly to the method proposed by Quiroga [18]. The soft mixing between a priori signals and the
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
presented concepts will be persued with specialists in the field. This method is particularly suited when pure averaging destroys the information, due to the nonergodicity of the underlying dynamics, such as in evoked brain potentials analysis (see [18] and reference therein). In this situation, the use of a few averages as a priori knowledge, and the recursive approach can be very effective. In the situation where an effective model is available to the user, such as P300 studies, our approach can also be of interest. Another potential application is foreseen in synchronous P300-based brain computer interfaces where triggered averaged signals would allow for denoising and as such enhance of the extracted features from the brain electrical signals (see [31] for a review of brain computer interface technologies). The use of a quantitative denoising measure such as NMSE provides only a gross indication of the performance of the method. We do recommend caution utilising this measure when optimising the parameters of the denoising algorithm. It may indeed not meet the needs from experts, and is application dependent.
6. Conclusions We have proposed a new method for denoising signals based on a priori knowledge and subspace decomposition. We have implemented the method using wavelet subspace decomposition. The smoothly adjustable denoising methodology makes it attractive for professionals who are working with highly noisy signals, which vary with time and from measurement to measurement (non-stationary and non-ergodic). Three parameters, i.e. the smoothing parameter a, the recursive smoothing parameter l and the Gaussian filter window width s, control the denoising procedure which turns it very easy to use by non-specialist. Several examples using real life signals and synthetic signals have provided further insights into the method.
Acknowledgements The authors would like to thank Dr. K. Le for his careful reading of the manuscript, and the Geneva University Hospital, Functional Brain Mapping Laboratory, for providing the visual evoked potentials in the Swiss IM2 NCCR framework.
2241
References [1] M.B. Priestley, Non-linear and Non-stationary Time Series Analysis, Academic Press, New York, 1988. [2] L.L. Scharf, Statistical Signal Processing—Detection, Estimation, and Time Series Analysis, Addison-Wesley, Reading, MA, 1991. [3] S.V. Vaseghi, Advanced Signal Processing and Digital Noise Reduction, Wiley, New York, 1996. [4] R. Vautard, P. Yiou, M. Ghil, Singular-spectrum analysis: a toolkit for short, noisy chaotic signals, Physica D 58 (1992) 95–126. [5] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, New York, 1998. [6] T.-P. Jung, C. Humphries, T.W. Lee, M.J. McKeown, V. Iragui, Makeig, T.J. Sejnowski, Removing electroencephalographic artifacts from by blind source separation, Psychophysiology 37 (2000) 163–178. [7] H. Kantz, T. Schreiber, Nonlinear time series analysis, Nonlinear Science Series, vol. 7, Cambridge University Press, Cambridge, 1997. [8] R. Cawley, G. Hsu, Local-geometric-projection method for noise reduction in chaotic maps and flows, Phys. Rev. A 46 (1992) 3057–3082. [9] P.-F. Marteau, H.D.I. Abarbanel, Noise reduction in chaotic time series using scaled probabilistic methods, J. Nonlinear Sci. 1 (1991) 313–343. [10] S. Akkarakaran, P. Vaidyanathan, Results on principal component filter banks: colored noise suppression and existence issues, IEEE Trans. Inform. Theory 47 (2001) 1003–1020. [11] G. Matz, F. Hlawatsch, A. Raidl, Signal-adaptive robust time-varying Wiener filters: best subspace selection and statistical analysis, Proceedings of the IEEE ICASSP, vol. 1, 2001, pp. 3945–3948. [12] S.P. Ghael, A.M. Sayeed, R.G. Baraniuk, Improved wavelet denoising via empirical Wiener filtering, Proceedings of the SPIE, Mathematical Imaging, San Diego, vol. 3169, 1997, pp. 389–399. [13] R. Vetter, J.-M. Vesin, P. Celka, Ph. Renevey, J. Krauss, Automatic nonlinear noise reduction using local component analysis and MDL parameter selection, EMBS 2003, Cancun, Mexico, vol. 1, 2003, pp. 491–496. [14] M. Vetterli, J. Kovacevic, Wavelets and Subband Coding, Prentice-Hall, Englewood Cliffs, NJ, 1995. [15] L. Cohen, Time–Frequency Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1995. [16] D.L. Donoho, I. Johnstone, Adapting to unknown smoothness via wavelet shrinkage, J. Amer. Statist. Assoc. 90 (1995) 1200–1224. [17] D.L. Donoho, De-noising by soft-thresholding, IEEE Trans. Inform. Theory 41 (1995) 613–627. [18] R.Q. Quiroga, H. Garcia, Single-trial event-related potentials with wavelet denoising, Clin. Neurophysiol. 114 (2003) 376–390. [19] P. Celka, R. Vetter, Ph. Renevey, C. Verjus, V. Neuman, J. Luprano, J.-D. Decotignie, C. Piguet, Wearable biosensing: signal processing and communication architecture issues, J. Telecommun. Inform. Technol. (Special Issue on e-Health: Wearable devices, Personal Health Management Systems and Services based on biosensors), 2005, in preparation. [20] A.C. Guyton, J.E. Hall, Textbook of Medical Physiology, W.B. Saunders Company, London, 1984.
ARTICLE IN PRESS 2242
P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
[21] P.E. McSharry, G.D. Clifford, L. Tarassenko, L. Smith, A dynamical model for generating synthetic electrocardiogram signals, IEEE Trans. Biomed. Engrg. 50 (2003) 289–294. [22] A. Aldroubi, M. Unser, Wavelets in Medicine and Biology, CRC Press, Boca Raton, FL, 1996. [23] D. De Rossi, F. Carpi, F. Lorussi, A. Mazzoldi, R. Paradiso, E.P. Scilingo, A. Tognetti, Electroactive fabrics and wearable biomonitoring devices, AUTEX Res. J. 3 (2003) 180–185. [24] T. Schreiber, A. Schmitz, Surrogate time series, Physica D 142 (2000) 346–382. [25] V.X. Afonso, W.J. Tompkins, T.Q. Nguyen, S. Luo, Ecg beat detection using filter banks, IEEE Trans. Biomed. Engrg. 46 (1999) 192–202. [26] P. Renevey, R. Vetter, J. Krauss, P. Celka, Y. Depeursinge, Wrist-located pulse detection using IR signals, activity and nonlinear artifact cancellation, Proceedings of the IEEE EMBS 2001, 2001, pp. 1–4. [27] P. Renevey, R. Vetter, P. Celka, J. Krauss, Activity classification using HMM for improvement of wrist-located
[28]
[29]
[30]
[31]
pulse detection, Proceedings of the Biosignal 2002, 2002, pp. 192–196. A. Effern, K. Lehnertz, T. Schreiber, T. Grunwald, P. David, C.E. Elger, Nonlinear denoising of transient signals with application to event related potentials, Physica D 140 (2000) 257–266. A. Effern, K. Lehnertz, T. Grunwald, G. Ferna´ndez, P. David, C.E. Elger, Nonlinear denoising of transient signals with application to event related potentials, Psychophysiology 37 (2000) 859–865. A. Effern, K. Lehnertz, T. Grunwald, G. Ferna´ndez, P. David, C.E. Elger, Single trial analysis of event related potentials: nonlinear denoising with wavelets, Clin. Neurophysiol. 11 (2000) 2255–2263. T. Vaughan, W. Heetderks, L. Trejo, W. Rymer, M. Weinrich, M. Moore, A. Kbler, B. Dobkin, N. Birbaumer, E. Donchin, E. Wolpaw, J. Wolpaw, Brain-computer interface technology: a review of the second international meeting, IEEE Trans. Neural Systems Rehab. Engrg. 2 (2003) 94–109.
Signal Processing 86 (2006) 2233–2242 www.elsevier.com/locate/sigpro
Smoothly adjustable denoising using a priori knowledge Patrick Celkaa,,1, Elly Gyselsb a
School of Engineering, Griffith University, Gold Coast, Queensland Australia Control and Signal Processing Section, Swiss Center for Electronics and Microtechnology, Rue Jaquet Droz 1, 2007 Neuchaˆtel, Switzerland
b
Received 7 April 2005; received in revised form 30 September 2005 Available online 16 November 2005
Abstract Assuming a signal is composed of information and noise, this paper presents a generic approach to denoising by mapping the noisy signal using a priori information about the signal to be retrieved. The method is based on a subspace decomposition of both the a priori information at disposal and the noisy signal, followed by shrinkage of both subspace coefficients and smooth mapping of first onto the second space. Our method propose 3 different ways of building the a priori knowledge: A model-base, averaging-based and recursive-based. The proposed method is particularly suited for low signal to noise ratio. The denoising methods are validated on synthetic electrocardiogram (ECG) signals and further assessed on real life ECG and visual brain evoked potentials using the wavelet transform. r 2005 Elsevier B.V. All rights reserved. Keywords: Subspace denoising; Wavelet transform; Biosignal processing; Electrocardiogram; Evoked potential
1. Introduction For several decades, scientists and engineers have developed strategies to reduce the noise from signals. Well known techniques are Wiener filtering [1,2], spectral analysis [3], model-based analysis [3], principal component analysis [4], wavelets [5], independent component analysis [6], nonlinear techniques [7], local subspace projection [8], probabilistic approaches [9], and combinations of these methods [10–13]. Only some of them, such as Corresponding author. Tel.: +61 7 5552 8459;
fax: +61 7 5552 8065. E-mail address: p.celka@griffith.edu.au (P. Celka). 1 The work has partly been done at the Swiss Center for Electronics and Microtechnology, and was supported by the Swiss IM2 NCCR Program.
principal component analysis (local or global) and wavelets, are suited for non-stationary signals. Most of these denoising techniques use a particular subspace decomposition operator H ¼ fH k gM k¼1 that maps the signal þ noise space S into ~ k , for k ¼ 1; . . . ; M, a particular set of subspaces S where noise and signal components are expected to be as separated as possible in a given normed space. ~ k , a weighting operator Pk Then, in each subspace S reduces the noise components. We will call the entire space weighting operator P ¼ fPk gM k¼1 . The weighting operation is often performed by thresh~ k or by olding of the components in the space S Wiener filtering [12]. The weighted components now ~ c . Finally, by inverse transform span the spaces S k M H1 ¼ fH 1 g , k k¼1 the signal is reconstructed in the signal space Sc 2 S.
0165-1684/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2005.10.003
ARTICLE IN PRESS 2234
P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
Fig. 1. Subspace denoising in S which shows the action of the subspace operator H and the weighting operator P.
A measurement signal xðnÞ (n is the discrete time index) in S can be written as xðnÞ ¼ sðnÞ þ ZðnÞ
with n ¼ 0; . . . ; N 1,
(1)
where sðnÞ is the information signal component we want to retrieve, and ZðnÞ is the noise signal component. The denoising scheme is displayed in Fig. 1. The following signals are generated through the denoising procedure: ~ xðnÞ ¼ HðxðnÞÞ, ~ s~ðnÞ ¼ PðxðnÞÞ,
ð2Þ ð3Þ
s^ðnÞ ¼ H1 ð~sðnÞÞ,
ð4Þ
where s^ðnÞ is the estimated information signal. The choice of the subspace operator H is most often an open question and depends on the application, the underlying generating system linearity, the stationarity, the available computing and memory resources and the signal statistics. Generally, orthogonal filter banks, such as principal components filter banks [10] or discrete orthogonal wavelet transforms [14], are preferred because of their numerous advantages such as perfect reconstruction of the signals, and fast implementation. The use of the wavelet decomposition allows for a natural subband analysis of the signals under inspection and a good joint time–frequency decomposition [5,15]. This implies that denoising using wavelet decomposition will lead to the best performance for non-stationary signals. The choice of the wavelet basis depends largely on the time and frequency behaviour of the information we want to extract. Denoising using wavelets usually requires soft or hard thresholding on the wavelet coefficients [5,16,17]. Determining the best threshold value is a matter of looking more closely at the signals at hand and the purpose of the denoising. We will use the discrete time wavelet subspace operator H ¼ W in
the following, but our approach can be implemented with any subspace decomposition. In some cases, we are not totally blind for performing the denoising. We may have some information given a priori such as the signal or noise statistics in the time or frequency domains. We know that the best denoising filter for stationary time series with a priori knowledge of the noise and signal statistics is the Wiener filter. Even if we have only an approximation of the signal and noise statistics at disposal, and the signals are nonstationary, we may do better than just blind denoising. Our approach is based on a smooth mapping of the subspace component of the given a priori signal sa ðnÞ and that of the signal xðnÞ. The smooth mapping allows the user to choose different levels of denoising. This is of great interest in medical applications where the signal-to-noise ratio (SNR) can be very low, such as evoked potential (EP) studies [18], electrocardiogram (ECG) and photoplethysmograph measurements recorded with wearable sensors [19]. The paper is organised as follows: Section 2 formulates the smooth mapping approach. Section 3 discusses the construction of the a priori knowledge signals. In Section 4, results based on synthetic and experimental signals are presented. Section 5 discusses the results and Section 6 concludes the paper. 2. Smooth mapping Let us assume that we have at disposal an a priori signal sa ðnÞ which contains most of the time and frequency information of the signal sðnÞ, and a measurement signal xðnÞ. We transform both signals xðnÞ and sa ðnÞ through the subspace operator H ¼ W ~ xðnÞ ¼ W ðxðnÞÞ, a s~ ðnÞ ¼ W ðsa ðnÞÞ.
ð5Þ ð6Þ
~ The signals xðnÞ and s~a ðnÞ are the discrete wavelet transforms of the signals xðnÞ and sa ðnÞ, and the subspace operator W is the N N wavelet transform matrix [14]. Eq. (5) can be rewritten in a matrix form as x~ ¼ W ðxÞ, s~ a ¼ W ðsa Þ,
ð7Þ ð8Þ
where the bold variables stand for the column vectors of samples from the corresponding signals’
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
spaces. The coefficients’ amplitude and positions in the vector s~ a are very close to the one we should find in s~ , the wavelet transform of the signal sðnÞ. Quiroga [18] has proposed to manually locate the most significant coefficients in s~ . The positions ni of the selected coefficients s~ðni Þ are stored in a set I M ¼ fni gK i¼1 . Using this set of positions I M , he proposed to locate the same coefficients in x~ and subsequently zeroed the other ones. This corresponds to our previously defined weighting matrix operator P which can be formulated as P ¼ DiagðpðkÞÞ
for k ¼ 1; . . . ; N
(9)
with pðkÞ ¼ 1 for k 2 I M , ¼ 0 for keI M ,
ð10Þ ð11Þ
where DiagðÞ is a diagonal matrix. Inverse transforming Px~ results in the denoised vector s^ ¼ W
1
~ Px.
(12)
Note that this weighting matrix is different from a simple thresholding technique because of the user manual selection. On the one hand, manually selecting the positions I M results in a supervised denoising and leads to optimal solutions in terms of expert knowledge. On the other hand, it is a drawback because of the subjectivity of the selection, the potential human mistakes from the ‘getting bored’ effect, and thus the lack of reproducibility. This method cannot be applied when a fully automatic denoising procedure is required. Also, this approach assumes that the signal sa ðnÞ is very similar to sðnÞ in the time and frequency domain, which may not always be the case. We propose a flexible approach which enables the user to adjust the denoising to his needs by adjusting very few parameters. First we locate the most significant coefficients in s~a by a wavelet soft thresholding operator PWs S [12], resulting in the set of positions I S for non-zero coefficients. Second, we locate the coefficients in x~ which correspond to the positions in I S and we apply a FIR Gaussian filter 2 2 GðkÞ ¼ ek =s , for k ¼ 0; . . . ; K 1, which is summarised in the PWx S operator. The Gaussian filtering in (16) allows the user to account for a possible time shift of the signal s~ compared to s~a , and is performed using a forward–backward approach for cancelling its phase effect. This results in the following two-step nonlinear weightings, coefficient selection and filtering operations:
2235
First step: The signal s~a is softly thresholded using ~a , s~aS ¼ PWs S s
(13)
where the weighting operator is given by Ws PWs S ¼ DiagðpS ðkÞÞ
(14)
and the soft thresholding nonlinear function expressed by t Ws a pS ðkÞ ¼ Yðj~s ðkÞj tÞ 1 a , (15) j~s ðkÞj where Y is the Heaviside function, t is the threshold. Second step: The resulting non-zero coefficients’ positions in s~aS are stored in the set I S . The set I S is used for the selection of the coefficients in x~ S as finally expressed by the following equation: ~ x~ S ¼ PWx S x,
(16)
where the selection and filtering operator is given by Wx PWx S ¼ GðkÞc DiagðpS ðkÞÞ
(17)
with c the convolution operator, and pWx S ðkÞ ¼ 1 ¼0
for k 2 I S , for keI S .
ð18Þ
The selection operation expressed in (18) allows for the extraction of small information bearing signal components embedded in x. Note that this selection procedure is different from hard or soft thresholding, and is nonlinear by nature. Finally, we apply a soft mixing of the coefficients s~aS and x~ S x~ Smooth ¼ ax~ S þ ð1 aÞ~saS ,
(19)
where a is the mixing parameter. The value of s in the Gaussian filter depends on how confident we are in the phase information contained in the a priori signal sa ðnÞ: the more confident we are, the smaller the width. The value of a reflects our confidence in the amplitude of the a priori signal sa ðnÞ. Using a ¼ 0 in the third step implies that x~ Smooth ¼ s~ aS , which means that we are fully confident in our a priori knowledge. Using a ¼ 1 results in x~ Smooth ¼ x~ S , which uses any a priori knowledge and signifies that we have no confidence at all in it. While the phase information can be of great importance in some applications such as brain EP (see [18] and references therein). We will focus on the amplitude in this work. Determining a is the main issue in this approach and is left to the expertise of the user, which renders this method flexible and easy to use. Yet, in the result section we present some insights
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
2236
about the performance of the method for different values of a which could help the choice of the user. The estimated information signal is recovered as follows s^ ¼ W 1 x~ Smooth ¼W
1
ðaPWs S Wx
ð20Þ þ ð1
a aÞPWs S W s Þ.
ð21Þ
3. Building the a priori 3.1. Averaging approach In most applications, users have some a priori knowledge about the nature of the system or signals under study. For instance, if one wants to perform the denoising of EP in the brain, the user may then record a few such signals sai ðnÞ, for i ¼ 1; . . . ; R, prior to the denoising. The a priori signal saA ðnÞ can then be obtained by averaging these recordings P a s ðnÞ a . (22) sA ðnÞ ¼ i i R In this particular case, the true EP saT ðnÞ is not known, but approximated by the averaging, which still is a noisy and distorted version of the true EP. The validity of the averaging approach is based on the assumption that the true signal to be recovered saT ðnÞ is deterministic, time-invariant, with added zero-mean noise. This method is demonstrated in Section 4.2.2.
i being the time index of the pattern. Such a signal can be, for example, an ECG and its PQRST complex the pattern [20,21] (see Fig. 3). Second, the user is making repetitive measurements sa ðn; iÞ on the system, with i being the repetition index. The first case involves time repetitiveness and each pattern has to be detected prior to building the a priori. Note that the repetitiveness does not need to be periodical or quasi-periodical. The second case concerns an experimental repetitiveness and each separate measurement can be used to build the a priori. In both cases, the recursive a priori signal approach is performed by an autoregressive moving average procedure saR ðn; iÞ ¼ lsaR ðn; i 1Þ þ ð1 lÞsa ðn; iÞ.
(23)
The parameter 0:5plo1 controls the recursive averaging memory. This approach is assessed in Section 4.1.2. 4. Results The method presented in the previous section is tested on synthetic and real life signals. The wavelet subspace operator W is used. The mother wavelet used in these experiments, quadratic biorthogonal B-splines [22], has been chosen to match the signal at hand: a set of biphasic smoothly varying waveforms concatenated in a so-called PQRST complex (see Fig. 3). The Gaussian window width was fixed to s ¼ 1 throughout this section. A wavelet decomposition up to 5 scales was used.
3.2. Model approach When a reliable model of the system or signals under inspection is at disposal, the a priori signal is given by saM ðnÞ. Again, the model is a distorted version of the true information signal saT ðnÞ. One advantage of this situation is that the signal saM ðnÞ is noiseless and can be produced by a time-varying system. However, the user is bound to its model and its consistency for the data at hand. This approach can also be used in combination with the averaging and the recursive methods. This approach is assessed in Section 4.1.1.
4.1. Synthetic signals We compare the recursive and model-based approaches with synthetic ECG data. The synthetic ECG is produced using the nonlinear dynamical model proposed by McSharry et al. [21] and provides us with saM ðnÞ. The parameters used in our simulation are those reported in [21], with a mean heart rate of 70 beats per minute, a lowfrequency to high-frequency band power ratio of LF=HF ¼ 0:5, and a standard deviation of the heart rate of 5 beats per minute. The sampling frequency is 512 Hz.
3.3. Recursive approach The last way to build the a priori signal is to use a recursive approach. There are two ways to implement the recursion. First the signal to be denoised is a repetition of a slowly time-varying pattern sa ðn; iÞ,
4.1.1. Model-based denoising The synthetic ECG signal is used then as a model of an ECG, to which noise has been added to simulate real ECG measurements. The added noise is built using real measurement noise from an ECG
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
signal recorded using electroactive polymer-based fabric electrodes similar to those used in [23]. Subject movement produces highly noisy data. Control of the SNR is thus made possible. Noisy ECG data can then be generated at will using the nonlinear dynamical model and Fourier-based amplitude corrected surrogate noise time series [24]. Our results are quantified by the Normalised Mean Square Error expressed in percent:
2237
Noisy ECG Original ECG Denoised ECG Wavelet
R P
Q
S T
1
(24)
where E½ stands for the mathematical expectation, and signals are assumed to have zero mean. Fig. 2 shows the averaged NMSE in function of the SNR and a when using the model-based denoising approach. One hundred Monte Carlo runs have been generated, for each value of SNR and a, to compute the averaged NMSE. It should be noted that the averaging a priori method of Section 3.1 would have had produced the same qualitative results. We can observe from Fig. 2 that for low values of the NMSE, when the SNR is decreasing, a should be increased to reach a given NMSE. As discussed previously, a is the confidence level of the user on the a priori signal. Thus, truthful processing of the data should use some boundaries of a, and thus will automatically constrain the minimum
α
0.75
E½ðsaM ðnÞ s^ðnÞÞ2 NMSE ¼ 100 , E½ðsaM ðnÞÞ2 ðnÞ þ E½^s2 ðnÞ
0.5
0.25
0
0.1
0.2
0.3
0.4 0.5 Time [s]
0.6
0.7
0.8
0.9
Fig. 3. A series of model-generated ECG signals showing the complex P-, Q-, R-, S- and T-waves, with an SNR ¼ 3 dB (solid thin), together with the noiseless ECG (dotted-dashed) and denoised signals (solid bold), for 5 values of a.
reachable NMSE. For a fixed value of a, we can observe that the NMSE curve decreases faster for high SNR than for low SNR, and that the slope of the curve is inversely proportional to a. Fig. 3 shows a series of synthetic noisy ECG ðSNR ¼ 3 dBÞ, for a ¼ 0; 0:25; 0:5; 0:75; 1, together with the model ECG signal saM ðnÞ used for denoising, and the denoised signal. We can observe how the different ECG waves (namely P-, Q-, R-, S-, and T-waves) are recovered during the denoising process as we gradually increase a from 0 to 1. Special attention should be drawn towards the extraction of the P-, Q- and S-waves which are small compared to the R- and T-waves.
Fig. 2. This figure shows 100 Monte Carlo NMSE averaged for different values of SNR and a. Nine different iso-NMSE are shown.
4.1.2. Recursive-based denoising For the recursive approach of Section 3.3, an ECG of 60 s has been generated with the same parameters as previously described with an SNR ¼ 3 dB. Each ith peak R-wave has been detected using a method similar to the one proposed by Afonso [25]. We have used a 5 s sliding window with half overlap for
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
2238
α=0.9
α=0.9
λ=1
λ=0.95
λ=0.95
λ=0.9
λ=0.9 NMSE
ECG
λ=1
λ=0.85
λ=0.85 λ=0.8
λ=0.8
23
23.5
24
24.5 Time [s]
25
0
25.5
Fig. 4. A series of denoised model-generated ECG signals ðSNR ¼ 3 dBÞ showing P-, Q-, R-, S- and T-waves for 5 values of l, together with the noisy ECG signal.
4.2.1. Electrocardiogram denoising ECG electroactive polymer-based fabric electrodes measurements similar to those presented in [23] have been performed for assessing the proposed denoising method. The sampling frequency is 250 Hz. We have used a 5 s sliding window with half overlap for processing the signals. An ith peak R-wave symmetrical window of 500 ms was applied yielding the signal sa ðn; iÞ for n ¼ 1; . . . ; 125 and the number of QRS complexes was 97. During the recording, the subject was asked to perform different slow movements such as squeezing shoulders, trunk rotation and deep breathing. The ECG measurements display various artefacts linked to these movements as shown in Fig. 5. Figs. 5 and 6 show different results of the recursive approach using different values of l and a ¼ 0:9. In Fig. 5, the NMSE values have been computed on each detected PQRST complex using a 20-averaged PQRST complex saA ðnÞ E½ðsaA ðnÞÞ2 ðnÞ þ E½^s2 ðnÞ
.
(25)
30 40 Time [s]
50
60
70
λ=1
α=0.9
λ=0.95 λ=0.9 λ=0.85
ECG
4.2. Measured signals
E½ðsaA ðnÞ s^ðnÞÞ2
20
Fig. 5. NMSE computed for each of the detected PQRST complex for 5 different values of l together with the measured ECG signal.
processing the signals. An ith peak R-wave symmetrical window of 500 ms was applied yielding the signal sa ðn; iÞ for n ¼ 1; . . . ; 256 and i ¼ 1; . . . ; 75. Fig. 4 shows 2.5 s of the denoised ECG signals for five values of l, together with the noisy ECG.
NMSE ¼ 100
10
32
34
36
38
40 42 Time [s]
44
46
48
50
Fig. 6. A series of denoised measured ECG signals for 5 values of l, together with the measured ECG signal.
The detection method is based on the work of Afonso [25]. The NMSE displays fluctuations due to changes in the PQRST pattern. These changes can be functionally related to the heart, or to artefacts due to the subject’s movements. Note that as l approach 1, the fluctuations approach zero (see (23)). These fluctuations in the NMSE could be used for detection of abnormal PQRST complexes when the signals are not perturbed by artefacts. Motionrelated artefact perturbations can be detected using other type of sensors, e.g. accelerometers [26,27].
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
Fig. 6 clearly shows how the l parameter influences the denoising results. A value of l ¼ 1 should be avoided for studies related to fine analysis of pattern changes because of the infinite memory effect of the recursive averaging method. Any value of lo0:5 should also be avoided due to the too short-range memory effect of the recursive averaging and poor denoising results.
25
Original 400 trials Averaged Smooth Wavelet cleaned
20
α=0.3
15 10 EEG channel #19
4.2.2. Brain evoked potentials denoising In this section, we present the results of denoising on visual brain evoked potential (VEP) which is a sub-class of EP. A visual stimulus was presented to a subject from which an electroencephalogram (EEG) was recorded. The EEG was recorded using the standard 10–20 system electrode placement and we have used 47 channels. Careful visual inspection of the signals was performed to discard any recording with a bad skin contact which produced heavily artifactual data. The signals were sampled at 500 Hz and lasted for 800 ms with 100 ms of pre-stimulus. R ¼ 400 trials were recorded and further averaged to provide the signal saA ðnÞ (see Section 3.1). Neuro-scientists and -clinicians often undertake a large collection of visual stimulus triggered VEP response signals, and then perform a statistical averaging of the stimulus locked VEP. In recent years, it has been pointed out that the statistical averaging can in fact deteriorate the single trial informational signal. Moreover, due to the natural non-stationarity of the signals and the experimental procedure, the statistical average cannot be applied. For this reason, methods for using single-trial EP have been developed. Due to the low SNR of the single trials, noise reduction methods have been developed with more or less success [18,28–30]. The commonality between these 4 techniques is the use of thresholding on wavelet coefficients to perform the actual denoising, and sometimes the manually selected set of good wavelet coefficients. In our approach we use an automatic coefficient selection ~ with or without in the a priori signal dual space S, thresholding its coefficients, and two adjustable parameters, a and s, which give the confidence in the a priori knowledge. The major difficulty is to differentiate the noise from the information signal in low SNR brain signals. We report here the results of our denoising procedure when using the averaging a priori method. The average EVP saA ðnÞ is still the ‘gold standard’ of EP analysis by neurologists and for this
2239
5 0 -5 -10 -15 -20 -0.1
0
0.1
0.2
0.3 0.4 Time [s]
0.5
0.6
Fig. 7. One VEP as measured together with a 400-trials average and the denoised VEP for a ¼ 0:3.
reason, our smooth mapping approach is perfectly suited to this type of application. Figs. 7–9 show one EVP as measured, together with the 400-trial averaged signal saA ðnÞ, and the denoised signal for a ¼ 0:3; 0:6; 0:9. To perform the denoising, we have used saA ðnÞ as a priori signal as proposed in Section 3.1. A bi-phasic wave pattern at about 150 ms and a large positive deflection at about 300 ms from the onset of the stimulus can be clearly distinguished on the denoised VEP. Figs. 7–9 illustrate the smoothing concept: while a is approaching 1, the denoised signal is closer to the a priori, while when it is decreased toward 0, the denoised signal is closer to the single trial. Note the significant difference between the single trial and the average signal indicating the large inter-trial variability. The user can choose a number of trials Ro400 to build the average a priori signal. However, due to the large variability of single trials SNR, we recommend to select the best VEP for building the a priori.
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
2240
25
25
Original 400 trials Averaged Smooth Wavelet cleaned
20
Original 400 trials Averaged Smooth Wavelet cleaned
20
α=0.9
15
15
10
10 EEG channel #19
EEG channel #19
α=0.6
5 0
5 0
-5
-5
-10
-10
-15
-15
-20
-20
-0.1
0
0.1
0.2
0.3 0.4 Time [s]
0.5
0.6
-0.1
0
0.1
0.2
0.3 0.4 Time [s]
0.5
0.6
Fig. 8. One VEP as measured together with a 400-trials average and the denoised VEP for a ¼ 0:6.
Fig. 9. One VEP as measured together with a 400-trials average and the denoised VEP for a ¼ 0:9.
5. Discussion
measured denoised signal in the subspace domain enable the fine tuning readily. While fine tuning is the core of the method, we have shown how the global measure of NMSE can guide the user towards well-chosen initial parameters. The main parameters are: (1) the mixing parameter a, (2) the Gaussian width of the subspace filter s and (3) if the recursive methodology is used, the forgetting autoregressive smoothing factor l. The s parameter has not been studied thoroughly in this work and will be addressed in a further study. The influence of s is particularly important when timing (phase information) is involved, e.g. evoked brain potentials and heart rate variability analysis. Applications to synthetic and measured ECG signals have shown the potential of the method, and its efficiency for the extraction of the P-, Q- and S-waves which are usually very small compared to the R- and T-waves. Experimental illustrations of the method using the average a priori on visual EP has proved to be efficient in extracting the information contained in the biosignals. Further neurophysiological and cardiovascular validation of the
The approach we have developed for denoising signals is flexible due to the restricted number of free parameters and thus allows the user to adapt the method to the purpose of her/his study. The proposed methodology stands between a pure automatic denoising procedure, where the user has no control of the algorithm, and a highly tuneable approach, in which the user can effectively fine-tune the algorithm. However, the latter can result in a subjective choice of the tuning parameters and also require a strong knowledge of signal processing. Our method is based on the availability or construction of a priori signals, and we have shown how to construct such a priori knowledge. The method can be used with any subspace approach with some slight modifications. The use of this a priori knowledge makes the method robust to low SNR and at the same time allows the extraction of small information patterns hidden in the noise, similarly to the method proposed by Quiroga [18]. The soft mixing between a priori signals and the
ARTICLE IN PRESS P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
presented concepts will be persued with specialists in the field. This method is particularly suited when pure averaging destroys the information, due to the nonergodicity of the underlying dynamics, such as in evoked brain potentials analysis (see [18] and reference therein). In this situation, the use of a few averages as a priori knowledge, and the recursive approach can be very effective. In the situation where an effective model is available to the user, such as P300 studies, our approach can also be of interest. Another potential application is foreseen in synchronous P300-based brain computer interfaces where triggered averaged signals would allow for denoising and as such enhance of the extracted features from the brain electrical signals (see [31] for a review of brain computer interface technologies). The use of a quantitative denoising measure such as NMSE provides only a gross indication of the performance of the method. We do recommend caution utilising this measure when optimising the parameters of the denoising algorithm. It may indeed not meet the needs from experts, and is application dependent.
6. Conclusions We have proposed a new method for denoising signals based on a priori knowledge and subspace decomposition. We have implemented the method using wavelet subspace decomposition. The smoothly adjustable denoising methodology makes it attractive for professionals who are working with highly noisy signals, which vary with time and from measurement to measurement (non-stationary and non-ergodic). Three parameters, i.e. the smoothing parameter a, the recursive smoothing parameter l and the Gaussian filter window width s, control the denoising procedure which turns it very easy to use by non-specialist. Several examples using real life signals and synthetic signals have provided further insights into the method.
Acknowledgements The authors would like to thank Dr. K. Le for his careful reading of the manuscript, and the Geneva University Hospital, Functional Brain Mapping Laboratory, for providing the visual evoked potentials in the Swiss IM2 NCCR framework.
2241
References [1] M.B. Priestley, Non-linear and Non-stationary Time Series Analysis, Academic Press, New York, 1988. [2] L.L. Scharf, Statistical Signal Processing—Detection, Estimation, and Time Series Analysis, Addison-Wesley, Reading, MA, 1991. [3] S.V. Vaseghi, Advanced Signal Processing and Digital Noise Reduction, Wiley, New York, 1996. [4] R. Vautard, P. Yiou, M. Ghil, Singular-spectrum analysis: a toolkit for short, noisy chaotic signals, Physica D 58 (1992) 95–126. [5] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, New York, 1998. [6] T.-P. Jung, C. Humphries, T.W. Lee, M.J. McKeown, V. Iragui, Makeig, T.J. Sejnowski, Removing electroencephalographic artifacts from by blind source separation, Psychophysiology 37 (2000) 163–178. [7] H. Kantz, T. Schreiber, Nonlinear time series analysis, Nonlinear Science Series, vol. 7, Cambridge University Press, Cambridge, 1997. [8] R. Cawley, G. Hsu, Local-geometric-projection method for noise reduction in chaotic maps and flows, Phys. Rev. A 46 (1992) 3057–3082. [9] P.-F. Marteau, H.D.I. Abarbanel, Noise reduction in chaotic time series using scaled probabilistic methods, J. Nonlinear Sci. 1 (1991) 313–343. [10] S. Akkarakaran, P. Vaidyanathan, Results on principal component filter banks: colored noise suppression and existence issues, IEEE Trans. Inform. Theory 47 (2001) 1003–1020. [11] G. Matz, F. Hlawatsch, A. Raidl, Signal-adaptive robust time-varying Wiener filters: best subspace selection and statistical analysis, Proceedings of the IEEE ICASSP, vol. 1, 2001, pp. 3945–3948. [12] S.P. Ghael, A.M. Sayeed, R.G. Baraniuk, Improved wavelet denoising via empirical Wiener filtering, Proceedings of the SPIE, Mathematical Imaging, San Diego, vol. 3169, 1997, pp. 389–399. [13] R. Vetter, J.-M. Vesin, P. Celka, Ph. Renevey, J. Krauss, Automatic nonlinear noise reduction using local component analysis and MDL parameter selection, EMBS 2003, Cancun, Mexico, vol. 1, 2003, pp. 491–496. [14] M. Vetterli, J. Kovacevic, Wavelets and Subband Coding, Prentice-Hall, Englewood Cliffs, NJ, 1995. [15] L. Cohen, Time–Frequency Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1995. [16] D.L. Donoho, I. Johnstone, Adapting to unknown smoothness via wavelet shrinkage, J. Amer. Statist. Assoc. 90 (1995) 1200–1224. [17] D.L. Donoho, De-noising by soft-thresholding, IEEE Trans. Inform. Theory 41 (1995) 613–627. [18] R.Q. Quiroga, H. Garcia, Single-trial event-related potentials with wavelet denoising, Clin. Neurophysiol. 114 (2003) 376–390. [19] P. Celka, R. Vetter, Ph. Renevey, C. Verjus, V. Neuman, J. Luprano, J.-D. Decotignie, C. Piguet, Wearable biosensing: signal processing and communication architecture issues, J. Telecommun. Inform. Technol. (Special Issue on e-Health: Wearable devices, Personal Health Management Systems and Services based on biosensors), 2005, in preparation. [20] A.C. Guyton, J.E. Hall, Textbook of Medical Physiology, W.B. Saunders Company, London, 1984.
ARTICLE IN PRESS 2242
P. Celka, E. Gysels / Signal Processing 86 (2006) 2233–2242
[21] P.E. McSharry, G.D. Clifford, L. Tarassenko, L. Smith, A dynamical model for generating synthetic electrocardiogram signals, IEEE Trans. Biomed. Engrg. 50 (2003) 289–294. [22] A. Aldroubi, M. Unser, Wavelets in Medicine and Biology, CRC Press, Boca Raton, FL, 1996. [23] D. De Rossi, F. Carpi, F. Lorussi, A. Mazzoldi, R. Paradiso, E.P. Scilingo, A. Tognetti, Electroactive fabrics and wearable biomonitoring devices, AUTEX Res. J. 3 (2003) 180–185. [24] T. Schreiber, A. Schmitz, Surrogate time series, Physica D 142 (2000) 346–382. [25] V.X. Afonso, W.J. Tompkins, T.Q. Nguyen, S. Luo, Ecg beat detection using filter banks, IEEE Trans. Biomed. Engrg. 46 (1999) 192–202. [26] P. Renevey, R. Vetter, J. Krauss, P. Celka, Y. Depeursinge, Wrist-located pulse detection using IR signals, activity and nonlinear artifact cancellation, Proceedings of the IEEE EMBS 2001, 2001, pp. 1–4. [27] P. Renevey, R. Vetter, P. Celka, J. Krauss, Activity classification using HMM for improvement of wrist-located
[28]
[29]
[30]
[31]
pulse detection, Proceedings of the Biosignal 2002, 2002, pp. 192–196. A. Effern, K. Lehnertz, T. Schreiber, T. Grunwald, P. David, C.E. Elger, Nonlinear denoising of transient signals with application to event related potentials, Physica D 140 (2000) 257–266. A. Effern, K. Lehnertz, T. Grunwald, G. Ferna´ndez, P. David, C.E. Elger, Nonlinear denoising of transient signals with application to event related potentials, Psychophysiology 37 (2000) 859–865. A. Effern, K. Lehnertz, T. Grunwald, G. Ferna´ndez, P. David, C.E. Elger, Single trial analysis of event related potentials: nonlinear denoising with wavelets, Clin. Neurophysiol. 11 (2000) 2255–2263. T. Vaughan, W. Heetderks, L. Trejo, W. Rymer, M. Weinrich, M. Moore, A. Kbler, B. Dobkin, N. Birbaumer, E. Donchin, E. Wolpaw, J. Wolpaw, Brain-computer interface technology: a review of the second international meeting, IEEE Trans. Neural Systems Rehab. Engrg. 2 (2003) 94–109.