Single-trial evoked potential estimation using wavelets

Computers in Biology and Medicine 37 (2007) 463 – 473 www.intl.elsevierhealth.com/journals/cobm

Single-trial evoked potential estimation using wavelets夡 Zhisong Wanga , Alexander Maierb , David A. Leopoldb , Nikos K. Logothetisc , Hualou Lianga,∗ a School of Health Information Sciences, University of Texas Health Science Center at Houston, 7000 Fannin, Suite 600, Houston, TX 77030, USA b Unit on Cognitive Neurophysiology and Imaging, National Institute of Health, Building 49, Room B2J-45, MSC-4400, 49 Convent Dr., Bethesda, MD

20892, USA c Max Planck Institut für biologische Kybernetik, Spemannstrasse 38, 72076 Tübingen, Germany

Abstract In this paper we present conventional and translation-invariant (TI) wavelet-based approaches for single-trial evoked potential estimation based on intracortical recordings. We demonstrate that the wavelet-based approaches outperform several existing methods including the Wiener filter, least mean square (LMS), and recursive least squares (RLS), and that the TI wavelet-based estimates have higher SNR and lower RMSE than the conventional wavelet-based estimates. We also show that multichannel averaging significantly improves the evoked potential estimation, especially for the wavelet-based approaches. The excellent performances of the wavelet-based approaches for extracting evoked potentials are demonstrated via examples using simulated and experimental data. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Wavelet; Single-trial; Evoked potential; Local field potential (LFP); Multiunit activity (MUA); Translation-invariant; Structure-from-motion (SFM); LMS; RLS; Wiener

1. Introduction With less signal attenuation and reduced interferences and noise, intracortical recordings generally offer higher signal-tonoise ratios (SNR) than scalp recordings such as the electroencephalogram (EEG). Two types of signals, spikes and local field potential (LFP), can be derived from intracortical extracellular microelectrode recordings. Spiking activity has been extensively studied in brain research. However, spikes only provide information about the outputs of a small number of neurons in a brain area. LFP, on the other hand, arises largely from dendritic activity of populations of neurons and thereby is dominated by the excitatory synaptic inputs to a larger cortical area as well as intra-areal local processing. One major component of LFP is the evoked potential, which reflects coordinated neural ensemble activity associated with an external event. Evoked potentials offer an important source of information to study the neural basis of perception and behavioral choice. However, in 夡 This work was supported in part by the NIH R03 MH067776, R01 NS054314, R01 MH072034, and the Whitehall Foundation. ∗ Corresponding author. Tel.: +1 713 500 3914; fax: +1 713 500 3929. E-mail address: [email protected] (H. Liang).

0010-4825/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2006.08.011

addition to the evoked potential, LFP also contains potentials caused by background activity, which are unrelated to any specific event. Extracting evoked potentials from LFP is an essential task in both cognitive research and clinical applications. Perhaps the most widely used approach to remove the background activity is ensemble averaging, which averages across trials the potentials measured over repeated presentations of a stimulus. To achieve this goal, a large number of trials is generally required. However, methods requiring fewer trials are desirable to save acquisition time and avoid nuisance for the subjects. In addition, averaging also assumes that there are no variations at all in amplitude, latency or waveform of the evoked potential across many trials. This strong assumption may not be valid in practice [1]. Whenever this assumption is violated, important information regarding the dynamics of the cognitive process is lost. In [2], an optimal time-varying Wiener filter was proposed to obtain an estimate of the evoked potential for each trial. Here a very large filtering matrix with the same dimension as the number of the data samples in each trial was devised to minimize the mean squared error of the evoked potential estimate. If the evoked potential varies randomly from trial to trial, the Wiener filter becomes a minimum variance estimator and its evoked

464

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473

potential estimate approximates the expected evoked potential [3]. However, such an implementation of the Wiener filter requires that the number of trials be much bigger than the number of data samples in each trial to avoid illconditioning of the autocorrelation matrix. Generally speaking, this necessitates even more trials than those required by ensemble averaging due to the great amount of data samples. In [3], the authors considered time-varying adaptive filters such as least mean square (LMS) and recursive least squares (RLS), which overcomes the difficulty of the Wiener filter in [2] by choosing a filter order smaller than the number of data samples in each trial. In fact, we can also implement the Wiener filter on each trial in a similar way as that of LMS and RLS [4]. LMS, RLS and the Wiener filter all require a reference signal for each trial, which is generally calculated as the ensemble average of all trials except the trial being considered. Therefore, they suffer from the trial-by-trial variability effect since they rely on averaging multiple trials to obtain a good reference signal for each single trial. Fourier transform based spectral analysis of evoked potential began more than three decades ago (see e.g., [5]). However, because the frequency components defined by the Fourier transform have infinite time support, it cannot track well the transient variations of evoked potentials. Short-time Fourier transform (STFT) provides a means of joint time–frequency analysis by applying moving windows to the signal and Fourier transform the signal within each window [6]. The problem with STFT, however, is that signals with different spectral components are treated with the same frequency resolution since a constant window length (hence a constant time resolution) is chosen throughout the analysis. A suitable window length needs to be chosen to balance the trade-off between the frequency and time resolution, and at the same time, satisfy the underlying assumption that the signal in each of the moving windows is stationary. In comparison with Fourier transform-based techniques, wavelet transform-based methods offer new capacity and advantages. Having superior resolution in both time and frequency domains, the wavelet transform is particularly useful for the analysis of non-stationary signals. In fact, the wavelet transform has become indispensable and widely used in many applications including communications, speech and audio processing, image processing, geophysics, finance, medicine, and neuroscience. In particular, neurophysiological signals are generally non-stationary and therefore wavelet transform-based approaches are more suitable than the Fourier transform-based methods in extracting information from neurophysiological signals. In the 1990s, wavelet transform was introduced for evoked potential analysis of EEG [7–9]. Recently, the wavelet transform was applied for EEG evoked potential extraction by choosing a few wavelet coefficients [10]. One disadvantage of the method is that it is implemented in a supervised manner requiring a priori knowledge of the time and frequency ranges of the evoked potential. Although such knowledge is well-known in EEG, it is not available in intracortical recordings such as LFP. Our goal is to extract evoked potentials from intracortical recordings. Since the SNR of intracortical recordings is high, evoked potentials will be wavelet decomposed with large wavelet coefficients, whereas the ongoing background activ-

ity will be decomposed with small coefficients. Therefore, the evoked potential can be estimated by thresholding the wavelet coefficients using the thresholds proposed by Donoho et al. [11,12]. In addition to this conventional wavelet-based evoked potential estimator, we also consider the translation-invariant (TI) wavelet-based evoked potential estimator. The former may suffer from pseudo-Gibbs phenomena near the discontinuities, whereas the latter can overcome this problem [13]. In this paper we present conventional and TI wavelet-based approaches for single-trial evoked potential estimation. We demonstrate that the wavelet-based approaches outperform several existing methods including the Wiener filter, LMS, and RLS. We compare the TI wavelet-based and conventional wavelet-based estimates and demonstrate that the former has higher SNR and lower root mean square error (RMSE) than the latter. We also show that multichannel averaging significantly improves the evoked potential estimation, especially for the wavelet-based approaches. The rest of the paper is organized as follows. In Section 2, we present the materials and wavelet-based approaches. Simulated and experimental examples comparing the performances of the Wiener filter, LMS, RLS, conventional and TI wavelet-based approaches are given in Section 3. Finally, Section 4 contains the conclusions. 2. Materials and methods 2.1. Subjects and neurophysiological recordings Electrophysiological recordings were performed in a healthy adult male rhesus monkey. After behavioral training was complete, occipital recording chambers were implanted and a craniotomy was made. Intracortical recordings were conducted with a multielectrode array while the monkey was viewing structurefrom-motion (SFM) stimuli, which consisted of an orthographic projection of a transparent sphere that was covered with randomly distributed dots on its entire surface. Stimuli rotated for the entire period of presentation, giving the appearance of threedimensional structure. The monkey was required to indicate the choice of rotation direction (clockwise or counterclockwise) by pushing one of two levers. Correct responses for disparitydefined stimuli were acknowledged with application of a fluid reward. In the case of fully ambiguous stimuli, where no correct response can be externally defined, the monkey was rewarded by chance. The recording site was the middle temporal (MT) area of the monkey’s visual cortex, which is associated with visual motion processing. LFP was obtained by filtering the collected data between 1 and 500 Hz. Multiunit activity (MUA) was obtained by filtering the data between 500 Hz and 1 kHz. MUA was further rectified to compute power as a function of time. 2.2. Wavelet transform It is well known that the Fourier transform decomposes a signal in bases of sines and cosines, which have infinite durations. The wavelet transform, on the other hand, decomposes a

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473

signal in bases of the scaled and shifted versions of a mother wavelet, which is a special waveform with finite duration. At the low scales the mother wavelet is contracted to match the highfrequency components. At the high scales the mother wavelet is stretched to match the low-frequency components. Therefore, the wavelet transform can obtain simultaneous time and frequency localizations and is a powerful complement to the Fourier transform. While the wavelet transform is still limited to the general time–frequency resolution trade-off inherent to spectral analysis, it can offer higher temporal resolution at lower frequencies and vice versa, which suits well the 1/f spectral profile of evoked potentials. The continuous wavelet transform (CWT) of a onedimensional signal x(t) with respect to a mother wavelet (t) is defined as [14]    +∞ 1 t −b Wx (a, b) = √ x(t) dt, (1) a |a| −∞ where a, b (a, b ∈ R and a  = 0) are the scale and translation parameters, respectively. As its name implies, CWT has continuous scales and translations, which result in information redundancy and high computational load. In practice, the scale and translation parameters are discretized for fast implementation. In fact, the scales on the powers of two {2l }l∈Z are sufficient and this will yield the dyadic wavelet transform:    +∞ 1 t −n Wx (l, n) = √ dt. (2) x(t) 2l 2l −∞ For a wavelet centered at time zero and frequency f0 , the wavelet coefficients Wx (l, n) characterizes the signal x(t) around the time n and frequency 2−l f0 . Using dyadic wavelet transform and multiresolution analysis, a signal can be decomposed into a coarse approximation and several details at different scales (levels). 2.3. Wavelet-based evoked potential estimation Intracortical recordings such as LFP contain a mixture of evoked potentials and ongoing background activity. The former is the signal of interest to be estimated and the latter is the “noise” to be suppressed. Our goal is to remove the noise and recover the signal of interest. Since the SNR of the intracortical recordings is high, evoked potential will be decomposed with large wavelet coefficients, whereas the ongoing background activity will be decomposed with small coefficients. Therefore, evoked potential can be estimated by thresholding the wavelet coefficients. In addition to this conventional wavelet-based evoked potential estimator, we also consider the TI wavelet-based evoked potential estimator. The former may suffer from pseudo-Gibbs phenomena near the discontinuities, whereas the latter overcomes this problem [13]. To obtain the TI wavelet-based evoked potential estimate, we shift the data, threshold the shifted data, unshift the thresholded data, and then average the results for all shifting. The method can be efficiently implemented in O(N log2 (N )) floating point operations (flops) using a TI table, which is similar in structure to a wavelet packet table [13].

465

Since the TI wavelet-based method is applied to each trial independently, it will not be affected by other trials and is useful for single-trial evoked potential estimation. As shown in the examples later on, the TI wavelet-based evoked potential estimates have higher SNR and lower RMSE than the conventional wavelet-based estimates. Wavelet-based evoked potential estimation can be summarized as a three-step procedure including wavelet decomposition, nonlinear thresholding, and inverse wavelet reconstruction. Thresholding is a critical step for separating the signal from noise and it is dependent on the noise model. Most wavelet thresholding methods in the literature were proposed for the white Gaussian noise model. In [15], temporally correlated noise was studied and level-dependent thresholding was proposed. It was shown that there tends to be little or no correlation between the wavelet coefficients at different levels and the wavelet coefficients at each level are approximately white noise. Hence level-dependent thresholding follows naturally. For the TI wavelet-based evoked potential estimator, the threshold at the jth level has the following form [15]:  (3) j = j 2 loge (N log2 (N )) with j = MADj /0.6745,

(4)

where MADj is median absolute value estimated on the wavelet coefficients at the jth level. For the conventional wavelet-based evoked potential estimator, the threshold at the jth level can be written as  j = j 2 loge (N ). (5) Two kinds of thresholding methods are popular: hard thresholding and soft thresholding. Hard thresholding sets all the wavelet coefficients with the magnitudes less than the threshold to zero and maintains the other wavelet coefficients. Soft thresholding also sets the wavelet coefficients with the magnitudes less than the threshold to zero. However, it reduces the remaining coefficients in magnitude by the threshold. Since the estimates by hard thresholding may contain noisy spikes, we stick to soft thresholding in the paper. Let wj,k be the kth wavelet coefficient at the jth level. Soft thresholding is defined as  0 |wj,k | j , wˆ j,k = wj,k − j wj,k > j , (6) wj,k + j wj,k < − j . If multichannel measurements are available, spatial filtering can be performed to further improve the evoked potential estimation. If the electrodes are close to each other, the evoked potential of all channels will be well-aligned in time due to the peculiar relationship between cortical architecture and mass action signals such as LFP, i.e., the functional architecture of visual cortex that comprises local clustering of similar neuronal responses taken together with the spatial summation of microscopic currents changes to macroscopic voltage changes. This is especially the case in many intracortical microelectrode array recordings. Therefore, the multichannel single-trial evoked

466

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473

potential estimate can be obtained by averaging the singlechannel single-trial evoked potential estimates. By doing this, the contribution from the phase-locked single-trial evoked potential is maintained, while incoherent noise contributions are reduced. It is worthy to mention that multichannel averaging differs from single-channel ensemble averaging in that trial-bytrial variability will not affect the evoked potential estimate of the former. As shown in the examples later on, multichannel evoked potential estimates have higher SNR and lower RMSE than the single-channel estimates. 3. Results In this section, we provide simulated and experimental examples to compare the single-trial evoked potential estimation performances of the Wiener filter, LMS, RLS, conventional and TI wavelet-based approaches. For the wavelet-based methods, if not explicitly specified, mother wavelets are chosen from the Daubechies wavelet family due to its properties such as compact support, high vanishing moments for a given support width, and minimum-phase scaling functions [16]. We use level-dependent thresholding to deal with temporally correlated noise. Soft thresholding is employed to yield smooth estimates. The filter order for the Wiener filter, LMS and RLS is chosen to be 30. Also a reference signal need to be chosen for the Wiener filter, LMS and RLS. For the single-channel case, the reference signal is calculated as the ensemble average of all trials except the trial being considered. For the multichannel case, the reference signal is calculated as the average of all multichannel trials except the trial from the channel being considered. Note that the wavelet-based methods do not require a reference signal.

between the true evoked potential s(n), and the estimated evoked potential sˆ (n). The RMSE is computed as  RMSE =

N 1  (s(n) − sˆ (n))2 N

1/2 ,

(8)

n=1

where N is the number of data samples in each trial. The second measure is the SNR, which is computed as   N 2 n=1 s (n) SNR = 10 log10 N (dB). (9) 2 n=1 (s(n) − sˆ (n)) In the first example, we have a single channel with 20 trials, all of which contain the same evoked potential. Table 1 compares the LMS, RLS, Wiener filter, conventional waveletbased and TI wavelet-based approaches in terms of the mean and standard deviation (SD) of the RMSEs and SNRs of the evoked potential estimates of all trials. It is obvious from looking at the table that the wavelet-based methods obtain the most accurate estimates among all approaches. In addition, the TI wavelet-based method outperforms the conventional waveletbased method. Fig. 2(a) corresponds to a single trial of the noisy raw signal. Fig. 2(b)–(f) show the evoked potential estimates for the single trial obtained via the LMS, RLS, Wiener filter, conventional wavelet-based and TI wavelet-based methods, respectively. Note that the wavelet-based methods have much more smooth estimates than the other methods. Smoothness is an expected property of the evoked potential based on 8 6

3.1. Simulation results In the simulation, we use a specific waveform for the evoked potential as shown in Fig. 1. The background neural activity is simulated as a mixture of colored noise e1 (n) and white Gaussian noise e2 (n), which change across trials. The following autoregressive model is used for the colored noise e1 (n) [3]: e1 (n) = 1.5080e1 (n − 1) − 0.1587e1 (n − 2) − 0.3109e1 (n − 3) − 0.0510e1 (n − 4) + u(n),

(7)

where u(n) is a zero-mean white Gaussian noise. The total SNR is 6 dB. We generate 20 trials of simulated data. Each trial contains 1500 data samples. Several merit measures are employed to assess the performances of different approaches. The first measure is the RMSE

Amplitude

4 2 0 -2 -4 -6 200

400

600

800

1000

1200

1400

Time samples

Fig. 1. Simulated evoked potential.

Table 1 Comparison of the different approaches for a single channel in the absence of evoked potential variation Approach

LMS

RLS

Wiener filter

Conventional wavelet

TI wavelet

Mean of RMSEs SD of RMSEs Mean of SNRs (dB) SD of SNRs (dB)

0.6072 0.0475 11.7242 3.7832

0.4914 0.0404 13.5634 5.4005

0.5084 0.0427 13.2719 5.2287

0.4396 0.0541 14.6339 8.3765

0.4100 0.0475 15.2254 8.8836

8

8

6

6

4

4 Amplitude

Amplitude

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473

2 0

2 0

-2

-2

-4

-4

-6

-6 200 400 600 800 1000 1200 1400 Time samples

200 400 600 800 1000 1200 1400

(b) 8

6

6

4

4

2 0

2 0

-2

-2

-4

-4

-6

-6 200 400 600 800 1000 1200 1400 Time samples

200 400 600 800 1000 1200 1400 Time samples

(d)

8

8

6

6

4

4 Amplitude

Amplitude

Time samples

8

Amplitude

Amplitude

(a)

(c)

467

2 0

2 0

-2

-2

-4

-4

-6

-6 200 400 600 800 1000 1200 1400

(e)

Time samples

200 400 600 800 1000 1200 1400

(f)

Time samples

Fig. 2. Simulated results for a single trial: (a) noisy signal; (b) LMS denoised signal; (c) RLS denoised signal; (d) Wiener filter denoised signal; (e) conventional wavelet denoised signal; and (f) TI wavelet denoised signal.

Table 2 Comparison of the conventional wavelet-based and TI wavelet-based approaches when mother wavelets are chosen from different wavelet families for a single channel in the absence of evoked potential variation Approach

Conventional wavelet

TI wavelet

Wavelet family

Daubechies

Symlets

Coiflets

Daubechies

Symlets

Coiflets

Mean of RMSEs SD of RMSEs Mean of SNRs (dB) SD of SNRs (dB)

0.4396 0.0541 14.6339 8.3765

0.4368 0.0545 14.6868 8.3172

0.4410 0.0546 14.6073 8.3651

0.4100 0.0475 15.2254 8.8836

0.4100 0.0474 15.2249 8.8805

0.4138 0.0482 15.1460 8.8075

the reliable evoked potential estimation results obtained by ensemble averaging over a large number of trials. Table 2 compares the wavelet-based and TI wavelet-based approaches in terms of the mean and SD of the RMSEs and SNRs of the evoked potential estimates of all trials, when mother wavelets

are chosen from the Daubechies, Symlets and Coiflets wavelet family, respectively. The similar results in RMSEs and SNRs demonstrate that the wavelet-based and TI wavelet-based approaches are not sensitive to the particular choice of the wavelet family.

468

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473

Table 3 Comparison of the different approaches for a single channel in the presence of evoked potential variation Approach

LMS

RLS

Wiener filter

Conventional wavelet

TI wavelet

Mean of RMSEs SD of RMSEs Mean of SNRs (dB) SD of SNRs (dB)

0.7353 0.1764 10.4973 5.9064

0.6430 0.2034 12.0032 8.3075

0.7059 0.2201 11.1979 7.5818

0.4396 0.0541 14.6339 8.3764

0.4100 0.0475 15.2260 8.8836

Table 4 Comparison of the different approaches for 20 channels LMS

RLS

Wiener filter

Conventional wavelet

TI wavelet

RMSE SNR (dB)

0.5343 12.7590

0.2788 18.4074

0.3107 17.4660

0.0936 27.8875

0.0859 28.6259

8

8

6

6

4

4 Amplitude

Amplitude

Approach

2 0

2 0

-2

-2

-4

-4

-6

-6 200 400 600 800 1000 1200 1400 Time samples

200 400 600 800 1000 1200 1400

(b)

Time samples

8

8

6

6

4

4 Amplitude

Amplitude

(a)

2 0

2 0

-2

-2

-4

-4

-6

-6 200 400 600 800 1000 1200 1400 Time samples

200 400 600 800 1000 1200 1400

(d) 8

6

6

4

4

2 0

2 0

-2

-2

-4

-4

-6

-6 200 400 600 800 1000 1200 1400

(e)

Time samples

8

Amplitude

Amplitude

(c)

Time samples

200 400 600 800 1000 1200 1400

(f)

Time samples

Fig. 3. Simulated results for the single trials from multichannels: (a) noisy signal; (b) LMS denoised signal; (c) RLS denoised signal; (d) Wiener filter denoised signal; (e) conventional wavelet denoised signal; and (f) TI wavelet denoised signal.

150

150

100

100

50

50 Amplitude

Amplitude

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473

0

0

-50

-50

-100

-100

-150

-150 100 200 300 400 500 600 700 800 900 Time samples

100 200 300 400 500 600 700 800 900

(b)

Time samples

150

150

100

100

50

50 Amplitude

Amplitude

(a)

0 -50 -100

0 -50 -100

-150

-150 100 200 300 400 500 600 700 800 900 Time samples

100 200 300 400 500 600 700 800 900

(d) 150

100

100

50

50

0 -50 -100

0 -50 -100

-150

-150 100 200 300 400 500 600 700 800 900

(e)

Time samples

150

Amplitude

Amplitude

(c)

469

Time samples

100 200 300 400 500 600 700 800 900

(f)

Time samples

Fig. 4. Experimental results for a single trial of LFP: (a) noisy signal; (b) LMS denoised signal; (c) RLS denoised signal; (d) Wiener filter denoised signal; (e) conventional wavelet denoised signal; and (f) TI wavelet denoised signal.

Next, we consider another single channel example where there is a random latency jitter in the evoked potential across the 20 trials. The jitter has a Gaussian distribution with zero mean and SD equal to the duration of 30 samples. The noise for each trial is the same as the corresponding trial in the first example. Latency jitter does exist in practice. For example, internally induced latency jitter occurs when the monkey judges the onset times of the same external stimulus differently at different trials. Table 3 compares the LMS, RLS, Wiener filter, conventional wavelet-based and TI wavelet-based approaches in terms of the mean and SD of the RMSEs and SNRs of the evoked potential

estimates of all trials. The wavelet-based methods maintain the same RMSE as in the first example. On the other hand, significant performance degradations arise in the Wiener filter, LMS and RLS due to the poor reference signal. In the third example, we have 20 channels, each representing a single trial. Table 4 compares the Wiener filter, LMS, RLS, conventional wavelet-based and TI wavelet-based approaches in terms of the mean and SD of the RMSE and SNR of the evoked potential estimates. Comparing Tables 1 and 4, we see that multichannel averaging significantly improves the evoked potential estimates, especially for the wavelet-based approaches

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473 200

200

150

150

100

100

50

50 Amplitude

Amplitude

470

0 -50

0 -50

-100

-100

-150

-150

-200

-200 -250

-250 100 200 300 400 500 600 700 800 900

(a)

Time samples

100 200 300 400 500 600 700 800 900

(b)

Time samples

Fig. 5. Experimental results for the single trials from multichannels of LFP: (a) noisy signal; and (b) TI wavelet denoised signal.

with an improved SNR of more than 10 dB. Again, the TI wavelet-based approach performs the best among all the methods. Fig. 3(a) shows the multichannel average of a single trial of the noisy signal. Fig. 3(b)–(f) show the multichannel average of the evoked potential estimates for the single trial obtained via the LMS, RLS, Wiener filter, conventional wavelet-based and TI wavelet-based methods, respectively. 3.2. Experimental results In the experiment, LFP and MUA data from seven channels are collected simultaneously. Fig. 4(a) corresponds to a single trial of the “noisy” LFP signal from a single channel. Fig. 4(b)–(f) show the evoked potential estimates for the single trial obtained via the LMS, RLS, Wiener filter, conventional wavelet-based and TI wavelet-based methods, respectively. We see that wavelet-based methods have much more smooth estimates than the other methods. Fig. 5(a) shows the single trials of the noisy LFP signals for all seven channels. The raw signal of each channel is denoted by a different color. Fig. 5(b) shows the evoked potential estimates of the single trials for all seven channels obtained via the TI wavelet-based method. Note that the evoked potential estimates are well-aligned in time and hence multichannel averaging is possible. Fig. 6(a) shows the multichannel average of a single trial of the noisy LFP signal. Fig. 6(b)–(f) show the multichannel average of the evoked potential estimates for the single trial obtained via the LMS, RLS, Wiener filter, conventional wavelet-based and TI wavelet-based methods, respectively. Clearly, there is an improved visualization of the evoked potential estimates by using multichannel average. Fig. 7(a) and (b) show the time–frequency distributions for a single trial of a noisy raw LFP signal and its TI waveletbased evoked potential estimate, respectively. The SFM stimuli appeared at 0 s. The same axis range for the amplitude is used here. Event-related activity is clearly noticeable in the time–frequency distribution of the TI wavelet-based evoked potential estimate, whereas such activity can hardly be seen

from that of the raw signal. Therefore, we conclude that the TI wavelet-based method can recover the evoked potential. Fig. 8(a) shows the ensemble average of 20 trials of a rectified single trial MUA signal for the preferred rotation direction and non-preferred rotation direction, respectively. The terms “preferred” and “non-preferred” refer to the neuronal activity when the monkey viewed the bistable SFM stimuli. The significantly larger neuronal response corresponds to the preferred rotation direction of the neuron, while the significantly smaller neuronal response corresponds to the non-preferred rotation direction of the neuron. The monkey indicated the rotation direction by pushing one of two levers. Fig. 8(b) shows the ensemble average of 20 trials of the TI wavelet-based evoked potential estimates for the preferred rotation direction and non-preferred rotation direction, respectively. The SFM stimuli appeared at 0 s. It is hard to discriminate the neural response of the two conditions in the noisy MUA signal, while there is a clear difference between the two conditions in the TI wavelet-based evoked potential estimates. Therefore, the TI wavelet-based method can help improve the prediction of the monkey’s behavioral choice.

4. Conclusions In this paper we have shown how to apply conventional wavelet-based and TI wavelet-based approaches for single-trial evoked potential estimation. We have shown that wavelet-based approaches outperform several existing methods including the Wiener filter, LMS, and RLS. In addition, unlike the other methods, wavelet-based approaches are immune to the effect of trial-by-trial evoked potential variations. We have compared the TI wavelet-based and conventional wavelet-based estimates and demonstrated that the former has higher SNR and lower RMSE than the latter. We have also shown that multichannel averaging significantly improves the evoked potential estimation, especially for the wavelet-based approaches. Simulated and experimental examples have demonstrated the effectiveness of the wavelet-based approaches for extracting evoked potentials.

150

150

100

100

50

50 Amplitude

Amplitude

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473

0

0

-50

-50

-100

-100

-150

-150 100 200 300 400 500 600 700 800 900 Time samples

100 200 300 400 500 600 700 800 900

(b)

Time samples

150

150

100

100

50

50 Amplitude

Amplitude

(a)

0

0

-50

-50

-100

-100 -150

-150 100 200 300 400 500 600 700 800 900 Time samples

100 200 300 400 500 600 700 800 900

(d) 150

100

100

50

50

0

0

-50

-50

-100

-100 -150

-150

100 200 300 400 500 600 700 800 900

100 200 300 400 500 600 700 800 900

(e)

Time samples

150

Amplitude

Amplitude

(c)

471

Time samples

(f)

Time samples

Fig. 6. Experimental results for the multichannel average of a single trial of LFP: (a) noisy signal; (b) LMS denoised signal; (c) RLS denoised signal; (d) Wiener filter denoised signal; (e) conventional wavelet denoised signal; and (f) TI wavelet denoised signal.

5. Summary With less signal attenuation and reduced interferences and noise, intracortical local field potential (LFP) offers a higher signal-to-noise ratio (SNR) than scalp electroencephalograms (EEG). One major component of LFP is the evoked potential, which is a neural signal that reflects coordinated neural ensemble activity associated with an external event. The evoked potential offers an important source of information to study

the relation between neural activity and behavior. However, in addition to the evoked potential, LFP also contains potentials caused by ongoing background activity, which are not related to the specific event. Extracting evoked potentials from LFP is an essential task in both cognitive research and clinical applications. In this paper we have shown how to apply conventional wavelet-based and translation-invariant (TI) wavelet-based approaches for single-trial evoked potential estimation. We have

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473

70 60 50 40 30 20 10 0 -10 -20 -30

450 Frequency (Hz)

400 350 300 250 200 150 100 50 0 -0.2

0

(a)

0.2

0.4

0.6

0.8

400 350 300 250 200 150 100 50 0 -0.2

0

(b)

Time (s)

70 60 50 40 30 20 10 0 -10 -20 -30

450 Frequency (Hz)

472

0.2

0.4

0.6

0.8

Time (s)

Fig. 7. Experimental results for the time–frequency plot of a single trial of LFP: (a) noisy signal; and (b) TI wavelet denoised signal.

Non-preferred

12

Preferred

11 10 9 8 7 6 5 4 -0.4 -0.2

(a)

13

0

0.2

0.4

0.6

0.8

1

1.2

Ensemble average of neural response

Ensemble average of neural response

13

Time (s)

(b)

Preferred

11 10 9 8 7 6 5 4 -0.4 -0.2

1.4

Non-preferred

12

0

0.2 0.4 0.6 0.8

1

1.2 1.4

Time (s)

Fig. 8. Experimental results for the ensemble average of 20 trials of MUA: (a) noisy signal; and (b) TI wavelet denoised signal.

shown that wavelet-based approaches outperform several existing methods including the Wiener filter, least mean square (LMS), and recursive least squares (RLS). In addition, unlike the other methods, wavelet-based approaches are immune to the effect of trial-by-trial evoked potential variations. We have compared the TI wavelet-based and conventional waveletbased estimates and demonstrated that the former has higher SNR and lower root mean square error (RMSE) than the latter. We have also shown that multichannel averaging significantly improves the evoked potential estimation, especially for the wavelet-based approaches. Simulated and experimental examples have demonstrated the effectiveness of the wavelet-based approaches for extracting evoked potentials. References [1] W. Truccolo, K.H. Knuth, A.S. Shah, S.L. Bressler, C.E. Schroeder, M. Ding, Estimation of single-trial multi-component ERPs: differentially variable component analysis (dVCA), Biol. Cybern. 89 (2003) 426–438. [2] K.-B. Yu, C.D. McGillem, Optimum filters for estimating evoked potential waveforms, IEEE Trans. Biomed. Eng. 30 (1983) 730–737. [3] X.H. Yu, Z.Y. He, Y.S. Zhang, Time-varying adaptive filters for evoked potential estimation, IEEE Trans. Biomed. Eng. 41 (11) (1994) 1062– 1071. [4] S. Haykin, Adaptive Filter Theory, Prentice-Hall, Englewood Cliffs, NJ, 2001. [5] D.O. Walter, A posteriori Wiener filtering of average evoked response, Electroencephalogr. Clin. Neurophysiol. 27 (1969) 61–70.

[6] L. Cohen, Time Frequency Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1995. [7] E.A. Bartnik, K.J. Blinowska, P.J. Durka, Single evoked potential reconstruction by means of wavelet transform, Biol. Cybern. 67 (2) (1992) 175–181. [8] O. Bertrand, J. Bohorquez, J. Pernier, Time-frequency digital filtering based on an invertible wavelet transform: an application to evoked potentials, IEEE Trans. Biomed. Eng. 41 (1) (1994) 77–88. [9] R.Q. Quiroga, M. Schürmann, Functions and sources of evoked EEG alpha oscillations studied with the wavelet transform, Clin. Neurophysiol. 110 (1999) 643–654. [10] R.Q. Quiroga, H. Garcia, Single-trial event-related potentials with wavelet denoising, Clin. Neurophysiol. 114 (2003) 376–390. [11] D.L. Donoho, I.M. Johnstone, Ideal spatial adaptation via wavelet shrinkage, Biometrika 81 (1994) 425–455. [12] D.L. Donoho, De-noising by soft-thresholding, IEEE Trans. Inform. Theory 41 (1995) 613–627. [13] R.R. Coifman, D.L. Donoho, Translation-Invariant De-Noising, in: A. Antoniadis, G. Oppenheim (Eds.), Wavelets and Statistics, Lecture Notes in Statistics, Springer, New York, 1995. [14] C. Chui, An Introduction to Wavelets, Academic Press, San Diego, CA, 1992. [15] I.M. Johnstone, B.W. Silverman, Wavelet threshold estimators for data with correlated noise, J. R. Statist. Soc. Ser. B (Statist. Methodol.) 59 (1997) 319–351. [16] I. Daubechies, Recent results in wavelet applications, Proc. SPIE Int. Soc. Opt. Eng. 3391 (1998) 2–9. Zhisong Wang received the B.E. degree in electrical engineering from Tongji University, Shanghai, China, in 1998. He received the M.Sc. and Ph.D.

Z. Wang et al. / Computers in Biology and Medicine 37 (2007) 463 – 473 degrees in electrical engineering from the University of Florida, in 2003 and 2005, respectively. From 2001 and 2005, he was a Research Assistant in the Department of Electrical and Computer Engineering at the University of Florida, Gainesville. Since 2005, he has been a Postdoctoral Fellow in the School of Health Information Sciences at the University of Texas Health Science Center at Houston. His current research interests include biomedical signal processing, cognitive and computational neuroscience, wavelet analysis, time series analysis, pattern recognition, machine learning, array signal processing, and spectral estimation. He is a member of Sigma Xi and Phi Kappa Phi. Alexander Maier studied neurobiology at the University of Munich (Germany) and received his Ph.D. from the University of Tuebingen (Germany) for work at the nearby Max Planck Institute of Biological Cybernetics. He currently works at the National Institute of Mental Health in Bethesda, MD as a Visting Fellow. David A Leopold attained a B.S. in biomedical engineering from Duke University in 1991. He subsequently received his Ph.D. from Baylor College of Medicine in 1997, where he studied neurophysiological mechanisms of multistable perception. After postdoctoral work at the Max Planck Institute for Biological Cybernetics in Tübingen, Germany, he joined the National Institutes of Health in 2004. At the NIH, he is presently the chief of the Unit on Cognitive Neurophysiology and Imaging and the director of the Neurophysiology Imaging Core Facility. His researches combines electrophysiological and fMRI techniques to study mechanisms of visual perception.

473

Nikos K Logothetis is the Director of the Neurophysiology Division of the Max Planck Institute for Biological Cybernetics in Tübingen, Germany. He obtained a B.S. in mathematics at the University of Athens, a B.S. in biology at the University of Thessaloniki, and his Ph.D. in human neurobiology from the Ludwig Maximillians University in Munich. In 1985, he moved to the Brain and Cognitive Sciences Department of the Massachusetts Institute of Technology; in 1990, he joined the faculty of the Division of Neuroscience at Baylor College of Medicine. Seven years later, he moved to the Max Planck Institute to continue work on the physiological mechanisms that underly visual perception and object recognition. Since 1992, he has also been Adjunct Professor of Neurobiology at the Salk Institute in San Diego, and, since 1995, Adjunct Professor of Ophthalmology at the Baylor College of Medicine, Houston. His recent work includes the application of functional imaging techniques to monkeys, and measurement of how the functional magnetic resonance imaging signal relates to neural activity. Hualou Liang studied signal processing at Dalian University of Technology, Dalian, China, and received Ph.D. in Physics from Chinese Academy of Sciences, Beijing, China. He was postdoctoral researcher in Tel-Aviv University, Tel-Aviv, Israel, Max-Planck-Institute for Biological Cybernetics, Tuebingen, Germany, and the Center for Complex Systems and Brain Sciences, Florida Atlantic University, Boca Raton, USA. He is presently an Associate Professor in University of Texas Health Science Center, Houston, USA. His research experience throughout the years ranges from areas in biomedical signal processing to cognitive and computational neuroscience. He teaches graduate interdisciplinary courses including “Biomedical Signal Processing” and “Computational Biomedicine”.