{'fraction_non_alphanumeric': 0.029994962678023537, 'fraction_numerical': 0.012318541924256994, 'mean_word_length': 4.780307040762308, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Obtaining single stimulus evoked potentials with Wavelet Denoising

We present a method for the analysis of electroencephalograms (EEG). In particular, small signals due to stimulation, so called evoked potentials, have to be detected in the background EEG. This is achieved by using a denoising implementation based on the wavelet decomposition. One recording of visual evoked potentials, and recordings of auditory evoked potentials from 4 subjects corresponding to different age groups are analyzed. We find higher variability in older individuals. Moreover, since the evoked potentials are identified at the single stimulus level (without need of ensemble averaging), this will allow the calculation of better resolved averages. Since the method is parameter free (i.e. it does not need to be adapted to the particular characteristics of each recording), implementations in clinical settings are imaginable.

However, in many cases the success of such procedure is limited, especially when not many trials can be obtained or when characteristics of the evoked potentials change from trial to trial.
Several methods have been proposed in order to filter averaged EPs (see [12]). The success of such methods would imply the need of less number of trials and would eventually allow the extraction of single trial evoked potentials from the background EEG. Although averaging has been used since the middle 50's, up to now none of these attempts has been successful in obtaining single trial EPs, at least in a level that they could be applied to different type of EPs and that they could be implemented in clinical settings. Most of these approaches involves Wiener filtering [19,10] (or a "minimum mean square error filter" based on auto and cross-correlations in [14]) and have the common drawback of considering the signal as a stationary process. Since EPs are transient responses related with specific time and frequency locations, such time-invariant approaches are not likely to give optimal results. For this reason, de Weerd and coworkers [7,8] introduced a time-varying generalization for filtering averaged evoked potentials. The time-variant Wiener filter they proposed is clearly more suitable for the analysis of evoked potentials but with the caveats that such filter bank implementation does not give a perfect reconstruction, and that it is based on the Fourier Transform (therefore the signal being decomposed in bases of sines and cosines with the drawbacks that this imposes, as we will describe later).
These limitations, as well as the ones related with time-invariant methods can be solved by using the wavelet formalism. The Wavelet Transform is a time-frequency representation proposed first in [11], that has an optimal resolution both in the time and frequency domains and has been successfully applied to the study of EEG-EP signals [16,18,6,3]. The objective of the present study is to follow an idea originally proposed in [1], and to present a very straightforward method based on the Wavelet Transform to obtain the evoked responses at the single trial level. The key point in the denoising of evoked potentials is how to select in the wavelet domain the activity representing the signal (the EPs) and then eliminate the one related with noise (the background EEG). In fact, the main difference between our implementation and previous related approaches is in the way that the wavelet coefficients are selected. Briefly, such choice should consider latency variations between the single trial responses and it should not introduce spurious effects in the time range where the EPs are expected to occur. In this respect, the denoising implementation we propose will allow the study of variability between single trials, an information that could have high physiological relevance.
The paper is organized as follows: In section 2 the data sets to be analyzed are described. Section 3 gives the mathematical background of the Wavelet Transform and the multiresolution decomposition. The implementation of the denoising of EPs is presented in section 4, together with a discussion of the advantages of the present implementation in comparison with previous approaches. The application of the method to visual and auditory EPs is shown in sections 5 and 6, respectively. Finally, in section 7 the conclusions are drawn. 4 We will study evoked potential recordings obtained with an oddball paradigm upon two different stimulus modalities: visual evoked potentials (VEP) and auditory evoked potentials (AEP). The robustness of the method is stressed by the fact that although an (almost) identical implementation was used for both data sets, these were taken in different laboratories and under different recording settings.

Visual evoked potentials
Visual evoked potentials from one normal subject were obtained by using a checkerboard light pattern. Two different visual stimuli were presented in a pseudo-random order (oddball paradigm): 75% of the stimuli were the so called "non-target" (a color reversal of the checks) and the other 25% were the deviant stimuli or "target" (also a color reversal but with a half check displacement of the pattern). The subject was instructed to ignore the nontarget stimuli and to count the number of appearances of the target ones (see [16] for more details on the experimental setup). Scalp recordings were obtained from the left occipital (O1) electrode (near to the location of the visual primary sensory area) with linked earlobes reference. Sampling rate was 250 Hz and after bandpass filtering in the range 0.1 − 70Hz, ∼ 2sec of data (256 data pre-and 256 data post-stimulation) were saved on a hard disk. Inter-stimulus intervals varied randomly between 2.5 − 3.5 sec. The recording session consisted of 200 stimuli presentations. From a total of 50 target stimuli, after rejection of trials with artifacts (contamination of the recording with spurious activity; e.g. blinking) 30 were selected for further analysis.
Two evoked responses can be observed with this paradigm (see average response on the uppermost left plot in figure 2): First, a sensory related positive peak at about 100 ms after stimulation (P100) followed by a negative rebound (N200), that appear both upon non-target and target stimuli. Second, a positive peak at about 400-500 ms after stimulation (P300) appearing only upon target stimuli and related with the cognitive process of recognizing these stimuli as deviant. As with VEP, in the case of AEPs two main responses appear: the N100-P200 peaks (polarity is reversed in comparison to VEP) and a P300 response (but appearing earlier than in the VEP case, at about 300ms after stimulation).

From Fourier to Wavelets
The Fourier Transform of a given signal is defined as its inner product with complex exponential (sines and cosines) functions of different frequencies.
It allows a better visualization of the periodicities of the signal, especially when several frequencies are superposed. However, Fourier Transform gives no information about the time location of these periodicities and it further requires stationarity of the signal.
By tapering ("windowing") the complex exponential mother functions of the Fourier Transform, the Short Time Fourier Transform (STFT) gives a time evolution of the frequencies that can be obtained just by sliding the windows throughout the signal. The STFT gives an optimal time-frequency representation, but a critical limitation appears when windowing the data due to the uncertainty principle (see [4]). If the window is too narrow, the frequency resolution will be poor, and if the window is too wide, the time localization will be not so precise. Data involving slow processes will require wide windows while a narrow window will be more suitable for data with fast transients (high frequency components). Then, owing to its fixed window size, the STFT is not suitable for analyzing signals involving different range of frequencies, as in the case of EPs.
Grossmann and Morlet [11] introduced the Wavelet Transform in order to overcome this problem. The main advantage of wavelets is that they have a varying window size, being wide for low frequencies and narrow for the high ones, thus leading to an optimal time-frequency resolution in all the frequency ranges [4]. Furthermore, owing to the fact that windows are adapted to the transients of each scale, wavelets do not need stationarity.
The Continuous Wavelet Transform (CWT) of a signal x(t) ∈ L 2 (R) is defined as the inner product between the signal and the wavelet functions where C a,b are the wavelet coefficients and ψ a,b (t) are dilated (contracted) and shifted versions of a unique wavelet function ψ(t) (a, b are the scale and translation parameters, respectively). The CWT gives a decomposition of x(t) in different scales, tending to be maximum at those scales and time locations where the wavelet best resembles x(t). Contracted versions of ψ a,b (t) will match high frequency components of x(t) and on the other hand, dilated versions will match low frequency oscillations.
The CWT maps a signal of one independent variable t onto a function of two independent variables a, b. This procedure is redundant and not efficient for algorithm implementations. In consequence, it is more practical to define the Wavelet Transform only at discrete scales a and discrete times b by choosing the set of parameters {a j = 2 j , b j,k = 2 j k}, with j, k ∈ Z. We obtain then the discrete wavelet family that under appropriate conditions [4] forms a basis of L 2 (R), each wavelet function having a good localization in the time and frequency domains. In analogy with eq.(1) we define the Dyadic Wavelet Transform as For well defined wavelets it can be inverted, thus giving the reconstruction of x(t) whereψ is the dual of ψ (in case of orthogonal wavelets,ψ and ψ are identical). For this study, we choose quadratic biorthogonal B-Splines [5] as mother functions (see Fig.1) due to their similarity with the evoked responses. B-Splines are piecewise polynomials that form a base in L 2 [4,21]. We remark the following properties that make them optimal in signal analysis (see [21,4,5] for details): they are (anti-) symmetric, smooth, they have a nearly optimal time-frequency resolution and they have compact support.

Multiresolution decomposition
The information given by the dyadic Wavelet Transform can be organized according to a hierarchical scheme called multiresolution analysis [13]. If we denote by W j the subspaces of L 2 generated by the wavelets ψ j,k for each level j, the space L 2 can be decomposed as a direct sum of the subspaces W j , Let us define the subspaces The subspaces V j are a multiresolution approximation of L 2 and they are generated by scalings and translations of a single function φ j,k = φ( 2 −j n − k ) called the scaling function (see proof in [13]). Then, for the subspaces V j we have the complementary subspaces W j , namely: Let us suppose we have a discretely sampled signal x(n) ≡ a 0 (n) with finite energy. We can successively decompose it with the following recursive scheme where the terms a j (n) = k a j−1 (k) φ j,k (n) ∈ V j give the coarser representation of the signal and d j (n) = k a j−1 ψ j,k (n) ∈ W j give the details for each scale j = 0, 1, · · · , N. For any resolution level N > 0, we can write the decomposition of the signal as This method gives a decomposition of the signal that can be implemented with very efficient algorithms 1 due to the recursiveness of the decomposition (eq. (9)). Moreover, Mallat [13] showed that each detail (d j ) and approximation signal (a j ) can be obtained from the previous approximation a j−1 via a convolution with (FIR or truncated IIR) high-pass and low-pass filters, respectively.
De Weerd and coworkers [7,8]  ). 2 However we should remark that the multiresolution decomposition gives a perfect reconstruction (disregarding border effects) because the mother functions ψ j,k (t) of eq. (3) form a basis of L 2 (R). This is not to be confounded with the requirement of orthogonality as claimed in [3]. 4 Denoising of evoked potentials 4.1 Implementation of the method As proposed by Donoho [9], the conventional definition of denoising implies a thresholding criterion in the wavelet domain. The signal is recovered from noisy data just by setting to zero those wavelet coefficients below a certain threshold (hard denoising) or with the use of a smoother transformation (soft denoising). However, this procedure is not optimal for recovering the evoked potentials because these ones are of the order or even smaller than the background EEG. Therefore, instead of using a thresholding criterion, we im- we would not enough eliminate the background EEG activity in order to recognize the EPs (in the extreme case of keeping all coefficients we will just reconstruct the original signal). In this respect we propose to use test signals, such as an spontaneous EEG (or a pre-stimulus EEG segment as we will show in fig.5) in order to check for eventual spurious interpretations due to an unfortunate selection of the coefficients.
We heuristically found the selection of those coefficients remarked in black traces in Fig.2 an optimal compromise between EP resolution and sensitivity to variations between trials. As we will show in the next sections, they allow the visualization of the single trial EPs and also they cover a reasonable time and scale ranges where the EPs are physiologically expected to occur.
We should also mention that the wavelet coefficients to be kept could be smoothed by using, e.g. soft thresholding [9] in order to decrease border effects. Although the hard thresholding we used can introduce spurious border fluctuations (see for example in Fig.2  In summary the method consist in the following steps: 1. The averaged EP is decomposed by using the wavelet multiresolution decomposition. 2. The wavelet coefficients not correlated with the average evoked potential (but also considering a time range in which single trial EPs are expected to occur) are identified and set to zero.
3. The inverse transform is applied, thus obtaining a denoised signal.

Comparison with previous works
The criteria for choosing which wavelet coefficients are correlated with the signal and which ones with noise is the key feature of the different denoising implementations. This is in fact our main difference with previous approaches. As already mentioned, Donoho [9] proposed a thresholding criteria, something not suitable for separating evoked potentials from the background EEG activity. The aim of those methods using time-variant [19,10], timeinvariant [7,8] and wavelet based [3] Wiener filtering was to clean averaged EPs rather than to obtain single trial EPs 4 . Therefore, they did not consider latency variations between trials, one of the most important features obtained from single trial analysis in comparison with ensemble averaging.
A similar remark is applicable to the work of Demiralp et al [6]. These authors correlated one single wavelet coefficient with the P300 response, then using its sign for discriminating between trials with and without P300. By selecting trials in this way, they succeeded in obtaining better averages of the P300.
Mc Gillen and Aunon [14] and Bartnik et al [1], instead, aimed at obtaining single trial evoked potentials. The first authors proposed a filter based in auto-and cross-correlations that however, due to its time-invariance, is not optimal for the analysis of EPs. This is the main explanation for the spurious differences between the original and the filtered signals observed in their examples (figs. 3 and 4 in ref. [14]). Bartnik and coworkers introduced a denoising implementation similar to the one presented by us. However, based on correlation and discriminant analysis (using as input the wavelet coefficients from the average and the single trials) they proposed an automatic criterion for finding which wavelet coefficients best distinguishes the EP from the EEG. In comparison with our approach, this criterion has the the following drawbacks: First, it is not appropriate if more than one peak is present, as e.g. in the case of the P100-N200 and P300 responses to be shown in the next section. Second, it is in principle not sensitive to variations between trials (i.e. variations in latency will be cancelled in the average and consequently, they will be disregarded by the correlation analysis). Third, a denoising based in such correlation criterion is likely to alter the morphology of the EPs. In particular, in the example shown by Bartnik et al the denoised signals look much smoother than the original ones, in many cases worsening the recognition and time localization of the EPs (see fig.7 in ref. [1]). Note that with denoising (black curves) we can distinguish the evoked responses P100-N200 and P300 in most of the trials. Note also that these responses are not easily identified in the original signal due to their similarity with the ongoing EEG. We already can observe some variability between the trials, e.g. in the trials #2,#13 the EPs are practically not present. For easier visualization, in Fig.4 we plot the single trial evoked potentials with and without denoising (left and right side plot, respectively) by using contour plots. In the denoised plot we observe between 100−200ms a yellow/red pattern followed by a blue one corresponding to the P100-N200 peaks. These responses remain more or less stable during the whole recording session. The more unstable and wider yellow/red pattern at about 400−600ms correspond to the cognitive related P300 responses. Noteworthy, all these patterns are more difficult to recognize in the original signal (left plot).
In order to check that the better visualization obtained after denoising is not just an artifact of the method, in Fig.5 we show the results of the same analysis applied to an spontaneous EEG signal (without EPs). As the EEG signal, we took for each trial of the VEP recording set the data segments corresponding to the second previous to stimulation (we assume we can disregard pre-stimulus effects, such as expectation). In contrast to the previous figure we only observe some randomly appearing spots that do not form a coherent pattern as when the EPs are present.
Auditory target and non-target evoked potentials of a typical subject (427) are shown in Figs.6 and 7, respectively. In both figures we observe a blue pattern at about 100ms after stimulation corresponding to the N100, followed by a yellow/red pattern corresponding to the P200. Upon target stimuli we can also identify a yellow pattern at about 300ms corresponding to the P300 responses. The P300 is better defined in the first 5 trials, later it appears more spread and unstable.
We further study the variability of the N100 peaks upon non-target stimuli (variability of P300 upon target stimuli can be studied in a similar way).
We first obtained from the denoised responses the amplitude and latency of the N100 peaks , which where localized as the maximum in a time window between 50 − 180ms after stimulation. We remark that this is hard to implement in the original signal due to the presence of the ongoing EEG. We then calculated the mean and standard deviation of the peak amplitude and latency. The first three columns of Table 1 show the results for the four subjects studied and the control EEG signal (without EPs) already described in subject, which will be discussed later). However, we observe a higher variability both in amplitude and latency for the older subjects. These variations are even larger for the control EEG signal. We remark that our present goal is to illustrate how we can obtain information not seen in the EP averages rather than making any physiological conclusion. In this context, a detailed quantitative analysis of the variability between single trials will be reported elsewhere with a larger database of subjects and statistical validation.
We also calculated the cross-correlations between the single trials and the average for each subject, thus having a measure of "how good" the single trials resemble the average EP. Then we computed selective averages using only those trials whose cross-correlation with the average was larger than 0.4 5 . Although we will show selective averages for the original and denoised signals, cross-correlations were only computed with the denoised ones in order to avoid effects of random correlations of the background EEG. As implemented in [20] for pre-filtered trials, for the selected trials we calculated lag corrected averages, i.e. by shifting each trial ∆ i = t av − t i data points, where t av is the latency of the average and t i the one of each single sweep. The last two columns of table 1 reports the mean cross-correlation values and the number of trials selected for each subject. As expected, the EEG signal has a lower mean cross-correlation than the EPs.
In Fig.8  holds for the shifted averages, especially in the third case due to its larger latency variation (see Table 1). In the forth subject as in the case of the EEG, neither the selective and the shifted averages show a better identification of the EPs in comparison with the background activity. Note also that this is expected from their low cross-correlation values. Again, a more reliable dis-tinction between subjects with a very large EP variability and control EEGs would require the statistical analysis of a larger sample, something that it is beyond the scope of the present study. 21 We presented a method for extracting evoked potentials from the background EEG. It is not limited to the study of EPs/EEGs and similar implementations can be used for recognizing transients even in signals with low signal to noise ratio.