Time-Frequency Data Reduction for Event Related Potentials: Combining Principal Component Analysis and Matching Pursuit

Joint time-frequency representations offer a rich representation of event related potentials (ERPs) that cannot be obtained through individual time or frequency domain analysis. This representation, however, comes at the expense of increased data volume and the difficulty of interpreting the resulting representations. Therefore, methods that can reduce the large amount of time-frequency data to experimentally relevant components are essential. In this paper, we present a method that reduces the large volume of ERP time-frequency data into a few significant time-frequency parameters. The proposed method is based on applying the widely-used matching pursuit (MP) approach, with a Gabor dictionary, to principal components extracted from the time-frequency domain. The proposed PCA-Gabor decomposition is compared with other time-frequency data reduction methods such as the time-frequency PCA approach alone and standard matching pursuit methods using a Gabor dictionary for both simulated and biological data. The results show that the proposed PCA-Gabor approach performs better than either the PCA alone or the standard MP data reduction methods, by using the smallest amount of ERP data variance to produce the strongest statistical separation between experimental conditions.


I. Introduction
Event-related potential (ERP) signals measured at the scalp are produced by partial synchronization of neuronal field potentials across the cortex [1].This synchronization mediates the 'top-down' and 'bottom-up' communication both within and between brain areas and has particular importance during the anticipation of and attention to stimuli or events.Event related potentials (ERPs) are obtained by averaging EEG signals recorded over multiple trials or epochs time-locked to the particular stimulus.ERP signal analysis has proven to be effective in assessing the brain's current functional state and reflect many pathological processes (e.g., [2], [3], [4], [5], [6]).
Typically, ERP analysis is performed in the time domain, where the amplitudes and latencies of prominent peaks in the averaged potentials are usually measured and correlated with information processing mechanisms.However, this conventional approach has two major shortcomings.First, it is well-known that ERPs are transient and non-stationary signals.Second, ERPs generally contain multiple overlapping processes operating across different time and frequency ranges.A primary approach to this problem has been to utilize time-frequency signal representations to detect transient activity and to disentangle overlapping processes.Several methods exist to fulfill this goal including wavelet and wavelet packet decomposition [7], [8], [9], [10], [11], sparse signal representations using overcomplete dictionaries (such as matching pursuit [12], [13] and basis pursuit [14]), Cohen's class of time-frequency distributions [15], [16] and the recently introduced high resolution time-frequency distributions [17], [18], [19], [20].
Wavelet transforms have been successfully applied to the analysis of evoked potentials in a variety of studies [7], [4], [21].They have been shown to be advantageous over the Fourier transform, since the time varying frequency information can be observed.However, wavelets have well-known limitations in terms of time-frequency resolution tradeoff, i.e., at high frequencies, the temporal resolution is high whereas the frequency resolution is low and vice versa for low frequencies.Sparse representations such as matching pursuit and basis pursuit aim to find a 'best' fit to the given signal in terms of the elements of a redundant family of functions, called the dictionary [13], [12].The 'best' fit to the given signal is quantified through both the mean square error between the representation and the actual signal and the sparseness of the representation, i.e., the number of elements of the dictionary used in the representation should be minimal.This approach has the advantage of offering a fully quantitative description of the ERPs by parameterizing the time-frequency plane at the expense of being computationally expensive.Cohen's class of distributions provides advantages over the other time-frequency representations in that it accurately characterizes the physical time-frequency properties of a signal, e.g., energy and marginals, and yields uniformly high resolution over the entire time-frequency plane [15], [22].Recently, time-frequency distributions with improved resolution and concentration around the instantaneous frequency have been introduced such as the reassigned time-frequency representations, higher order polynomial distributions, and complex-lag distributions [20], [17], [23].Although these methods improve the resolution of the representations, they come at the expense of increased computational complexity and in some cases losing some of the desirable properties such as the marginals.Moreover, these distributions have been shown to be the most effective for polynomial phase signals, whereas ERPs have been shown to be well represented by damped sinusoids [24], thus the improvement provided by these more complex distributions would be minimal.For these reasons, in this paper we will focus on the Cohen's class of distributions, in particular the Reduced Interference Distributions.
The high resolution provided by Cohen's class of time-frequency distributions come at the expense of increased data.The application of these distributions to large sets of ERP data has tended to rely on a time-frequency region of interest (TF-ROI; region of interest on the TF surface) to define activity for evaluation.Therefore, there is a growing need for data reduction and feature extraction methods for reducing the three dimensional time-frequency surfaces of ERPs to a few parameters.The problem of feature extraction and data reduction has been traditionally addressed using parametric and non-parametric methods.Parametric approaches include sparse representations using overcomplete dictionaries [13], [12], [14], [25], [26], [27], [28], extraction of features from the time-frequency distributions such as the energy in different frequency bands, computation of higher order joint moments [29], [30] and entropy [31].Non-parametric data reduction methods, on the other hand, include datadriven multivariate component analysis such as the application of matrix factorization methods to time-frequency distributions.These methods include the non-negative matrix factorization (NMF) [32], [33], [34], singular value decomposition (SVD) [35], independent component analysis (ICA) [1], [36] and principal component analysis (PCA) [37], [38], [39], [40] to extract time-frequency features for classification purposes or for reducing the timefrequency surfaces to a few meaningful components.The application of the matrix factorization approaches have been mostly limited to decomposing a single time-frequency matrix into significant time and frequency components to reduce the dimensionality and extract features for consequent classification [34], [39].However, in ERP analysis there is a need for multivariate processing, i.e., it is important to extract components that describe a collection of signals, such as those collected over multiple channels or multiple subjects.The principal component analysis of time-frequency vectors representing multiple subjects described in [37], [41] addresses this issue by extracting time-frequency principal components over a collection of ERP waveforms.At this point, it is also important to motivate the use of PCA over other data factorization methods.PCA is a multivariate technique that seeks to uncover latent variables responsible for patterns of covariation in the data set and has been used widely for time domain ERP data description and reduction [42], [43].It is commonly applied to the covariance of the data matrix and is thus similar to SVD in the extracted components.PCA does not make any strong assumptions about the data, unlike NMF which imposes non-negativeness, with the only assumption being that the observations are linear functions of the extracted components which is a common assumption in ERP analysis.ICA has been proposed as a promising alternative to PCA for ERP data reduction [1], [36].However, recent comparisons of PCA with ICA for ERP data analysis indicates that ICA suffers from the component "splitting" problem, i.e., components that should not be separated are split into multiple components, and that it is more suitable for spatial decompositions rather than temporal ones [44], [45].Further, ICA has been most commonly applied to time-domain ERP signal representations, and its use with timefrequency ERP representations has not been well-validated.For these reasons, in the current paper we use PCA as the first step in our data reduction algorithm.
In this paper, we address the data extraction and reduction problem in the time-frequency plane by combining parametric and non-parametric methods in a non-stationary setting.The ultimate goal is to find time-frequency components that are common to a large set of ERP data and that can summarize the relevant activity in terms of a few parameters.We introduce a new data reduction method based on applying matching pursuit decomposition to the timefrequency domain principal components to further reduce the information from the principal components and to fully quantify the time-frequency parameters of ERPs.Since the principal components extracted from ERP time-frequency surfaces are well-localized in time and frequency, we propose quantifying them in terms of well-known compact signals, Gabor logons 1 , on the time-frequency plane.Even though there are various choices for the basis functions that can be used to decompose a given signal, in this paper Gabor logons are chosen for representing time-frequency structure of ERP signals for two major reasons.First, it is known that these functions achieve the lower bound of the uncertainty principle (time-bandwidth product) and have been described as the 'elementary signals' on the timefrequency plane [22], [46].Second, the parameters of the Gabor logons are well-suited for identifying between transient versus oscillatory brain activity as well as separating between overlapping time-frequency events with varying duration or frequency oscillation.They have been widely used in time-frequency representation of ERP signals [47], [48], [49], in particular EEG phenomena including sleep spindles [50], [13] and epileptic seizures [51].An algorithm similar to matching pursuit is developed in the time-frequency plane to determine the best set of logons that describe each ERP time-frequency principal component [12].Fitting Gabor logons to the extracted principal components offers three potential benefits.First, decomposing the principal components (PCs) into a few logons would capture the major activity described by that principal component while at the same time 1 In this paper, 'Gabor logons' and 'logons' will be used interchangeasbly.
serve as a tool of denoising, i.e. removing the unwanted noise or activity that may exist in the principal component.Second, insofar as a single logon can characterize the primary activity in the experimental manipulations for each principal component, this would offer evidence that the principal components approach is efficient at extracting compact timefrequency representations.Finally, the extracted logons offer an important unit of analysis in their own right, in that they are maximally compact by definition.The proposed methods are compared to both parametric and nonparametric data reduction methods in the timefrequency plane, namely, the standard matching pursuit algorithm [12], [52] and PCA in terms of efficiency, computational complexity and the effectiveness in describing the experimental effects in the data.To evaluate these methods, we employ both biological [41] and simulated data [37], that have been previously evaluated using the PCA on the timefrequency plane approach.
The rest of this paper is organized as follows.Section II gives a brief review of timefrequency distributions and various matching pursuit approaches.Section III introduces the data reduction method proposed in this paper, combining principal component analysis with matching pursuit on the time-frequency plane.Section IV details the data analyzed in this paper and presents the results of applying the proposed method to both simulated data and ERP signals.A comparison with different time-frequency data reduction methods is also given in this section.Finally, Section V concludes the paper and discusses the major contributions.

II. Background A. Time-Frequency Distributions
A bilinear time-frequency distribution (TFD), C(t, ω), from Cohen's class can be expressed as 2 [22]: where ϕ(θ, τ) is the kernel function in the ambiguity domain (θ, τ), x(t) is the signal, and t and ω are the time and the frequency variables, respectively.Some of the most desired properties of TFDs are the energy preservation, the marginals, and the reduced interference.For bilinear time-frequency distributions, cross-terms occur when the signal is multicomponent, i.e. if then , where C x i , x i and C x i , x j refer to the auto-terms and cross-terms, respectively.The cross-terms will introduce time-frequency structures that do not correspond to the time-frequency spectrum of the actual signal.For this reason, in this paper we will use reduced interference distributions (RIDs) that concentrate the energy across the auto-terms, satisfy the energy preservation and the marginals [53].

B. Matching Pursuit
The matching pursuit algorithm, originally proposed by Mallat, aims at obtaining the 'best' linear representation of a signal in terms of functions, {g i } i=1, 2,…,N (sometimes referred to as atoms), from an overcomplete dictionary, D, using an iterative search algorithm [12].
1. Define the 0th order residual as R 0 x = x, D ̃ = D and set k = 0.
2 All integrals are from −∞ to ∞ unless otherwise stated.

2.
For the kth order residual, R k x, select the best atom such that the inner product between the residual and the atom is maximized (2) 3. Compute the residue R k+1 x as (3) and go back to step 2 until a predetermined stopping criterion is achieved.The stopping criterion can either be a pre-selected number of atoms to describe the signal or a percentage of energy of the original signal described by the selected atoms.After M iterations, the following linear representation is obtained: This procedure converges to x in the limit, i.e. , and preserves signal energy.

C. Simultaneous Matching Pursuit
The principle of MP can easily be generalized to the simultaneous decomposition of multiple signals, X = (x 1 , x 2 , …, x r ), into atoms from the same overcomplete dictionary, D. This approach is sometimes referred to as the multichannel matching pursuit or the multivariate matching pursuit (MMP) algorithm in literature since it is usually applied to multiple signals collected over multiple channels or sensors [52], [54], [55], [56].In this paper, we will refer to this method as the simultaneous matching pursuit (SMP) to avoid any confusions since the method will be applied to multiple ERPs from different subjects and not from multiple channels.This algorithm can be described as follows: 1. Define for each signal x l the 0-th order residual as R 0 x l = x l and set D ̃ = D, k = 0.

2.
For the kth order residual, R k x l , select the best atom such that the sum of the squared inner products between the atom and the residual from each signal is maximized 3. Compute the residue R k+1 x l for each signal:

III. PCA-Gabor Method
Ideally, a time-frequency domain ERP data reduction method will faithfully reproduce established time and frequency based findings (i.e.peaks in the time domain such as P300 or summaries of frequency activity such as alpha), but also allow a more complex view of these phenomena using the joint time-frequency information available in the TFDs.The decomposition method used in this paper is based on two stages of consecutive data reduction.The first stage is a direct extension of PCA into the joint time-frequency domain and the second stage is the parametrization of the time-frequency principal components using a matching pursuit type algorithm.

A. PCA on the Time-Frequency Plane
The first stage of the algorithm extends principal component analysis to the time-frequency plane as follows: 1. Compute the time-frequency distribution of each ERP waveform from multiple subjects, x i , 1 ≤ i ≤ L: (8) where ψ is the discrete-time kernel in the time and time-lag domain and x i (n) is the ith ERP waveform.In this paper, the binomial kernel, given by for , is used as the time-frequency kernel.

2.
Given L ERP waveforms, rearrange the time-frequency surfaces into vectors and form the matrix where λ j is the eigenvalue of each principal component PC j .The principal components determine the span of the time-frequency space.

5.
Rotate the principal components using varimax rotation [57].Varimax rotation is an orthogonal transform that rotates the principal components such that the variance of the factors is maximized.This rotation improves the interpretability of the principal components.

6.
Rearrange each principal component into a time-frequency surface to obtain the ERP components in the time-frequency domain.
After the principal components on the time-frequency plane are extracted, they are ordered based on their eigenvalues and the most significant ones are used in the following parametrization stage.The number of principal components to keep is determined based on a normalized energy threshold.

B. Matching Pursuit on the Time-Frequency Plane
In this section, we introduce a matching pursuit type algorithm in the time-frequency domain to further parameterize the ERP time-frequency surfaces.The goal is to be able to describe the principal components using a compact set of time-frequency parameters using Gabor logons as the dictionary elements.The proposed algorithm is similar to the original matching pursuit [12] and the discrete Gabor decomposition [58], except that it is directly implemented in the time-frequency domain rather than in the time domain.This implementation is preferred over the standard MP for two reasons.First, the principal components are already in the time-frequency domain, and inverting them back to the time domain would increase the computational complexity.Second, this offers a way of directly modeling the time-frequency energy distribution.
An overcomplete dictionary of Gabor logons on the time-frequency plane is constructed by computing the time-frequency distribution of discrete time atoms where σ is the scale parameter, n 0 and k 0 are the discrete time and frequency shift parameters, respectively, and K is the total number of frequency samples.The elements of the dictionary, D, are the binomial TFDs of these atoms, G i (n, k).The number of elements in the dictionary are determined by the range of n 0 , k 0 and σ.In this paper, n 0 = 1, …, N, where N is the total number of time samples, k 0 = 1, …, K, where K is the total number of frequency samples, and σ = {1, 2, 4, …, 2 ⌊log 2 N⌋−1 }.
The proposed greedy search algorithm is an extension of the orthogonal matching pursuit (OMP) described in [59] to the time-frequency domain.The orthogonal matching pursuit adds a least-squares minimization to each step of MP to obtain the best approximation over the atoms that have already been chosen.This revision significantly improves the convergence speed of the algorithm.For a given time-frequency matrix, TFD, the search for logons that best describe the surface can be summarized as: 1. Initialize the residue as R 0 = TFD and set l = 0 and D ̃ = D.

2.
At the lth iteration, find the Gabor logon over the whole overcomplete dictionary, i.e. over all (n 0 , k 0 , σ), that has the largest inner product with the residue timefrequency surface, R l (11)

3.
Compute the approximation at the lth step, A l , as: where A ∈ span{G i , i = 1, 2, …, l}.This problem is solved using a least squares optimization approach.
4. Subtract the approximation, A l , from the residue to compute the new residue timefrequency distribution at the l + 1th iteration 6. Go back to step 2 until a pre-determined number of atoms is selected or the normalized mean squared error (NMSE) between the TFD and the approximation at the lth iteration is below a pre-determined threshold, i.e. ( NMSE is a measure of how close the approximation from the dictionary is to the original time-frequency distribution.Since the mean square error is normalized by the energy of the original TFD, it is always between 0 and 1.

A. Description of Biological Data
The biological data used in this paper has been previously presented utilizing PCA on the time-frequency plane approach, and thus we will only detail the relevant parameters here.The reader is directed to the previously published paper for greater detail [41].The sample consisted of twins in the Minnesota Twin Family Study (MTFS), a longitudinal and epidemiological investigation of the origins and development of substance use disorders and related psychopathology.All male and female twin participants for whom ERP data were available from the study's psychophysiological assessment served as subjects for this investigation.This sample combined subjects from the two age cohorts of the MTFS: subjects in one cohort were 17 years old at intake whereas subjects in the other were approximately 11 years old at intake.Data for this younger cohort came from a follow-up assessment conducted when subjects were approximately 17 years old.The sample thus comprised 2,068 17-year-old adolescents in all (mean age = 17.7;SD = 0.5; range = 16.7 to 20.0).
A visual oddball task was used.Each of the 240 stimuli comprising this task was presented on a computer screen for 98 ms, with the inter-trial interval (ITI) varying randomly between 1 and 2 s.A small dot, upon which subjects were instructed to fixate, appeared in the center of the screen during the ITI.On two-thirds of the trials, participants saw a plain oval to which they were instructed not to respond.On the remaining third of the trials, participants saw a superior view of a stylized head, depicting the nose and one ear.These stylized heads served as 'target' stimuli.Participants were instructed to press one of two response buttons attached to each arm of their chair to indicate whether the ear was on the left side of the head or the right.Half of these target trials consisted of heads with the nose pointed up, such that the left ear would be on the left side of the head as it appeared to the subject (easy discrimination).Half consisted of heads rotated 180 degrees so that the nose pointed down, such that the left ear would appear on the right side of the screen and the right ear would appear on the left side of the ear (hard discrimination).
For each trial, 2 s of EEG, including a 500 ms pre-stimulus baseline, were collected at a sampling rate of 256 Hz.EEG data were recorded from three parietal scalp locations, one on the midline (Pz) and one over each hemisphere (P3 and P4).Consistent with the previous report, only data from the Pz electrode is reported here.Similarly, although ERPs to standard (frequent) stimuli were collected, they were not analyzed for the current paper; target condition responses serve as the basis for all decompositions and analyses presented.Therefore, the analysis in this paper focuses on data reduction for ERPs collected across multiple subjects from a single channel.However, the methods developed can easily be extended to single subject and/or multiple channel data.
Principal component decompositions were employed to evaluate the proposed approach.For the purposes of this study, decompositions for condition-averaged data were conducted on narrow time and frequency ranges, to focus on lower frequency delta and theta activity.Condition-averaged ERPs were constructed separately for easy and hard discrimination conditions.These included frequencies ranging from 0 to 5.75 Hz and time ranging from stimulus onset to 1000 ms post-stimulus.The range was narrowed to focus on the timefrequency range containing the majority of variance: theta, delta and low frequency activity.

B. Description of Simulated Data
Two simulated datasets were employed in the current paper.As with the biological data, these datasets were employed previously with PCA approach alone [37].Briefly, the two sets included are: 3-logons and 3-logons with noise.All simulated sets were 100 Hz sampled signals of 1000 ms, with the first and last 100 ms discarded after the TFD is computed to remove edge effects.The first simulated dataset contains 3 logons with clearly separated time and frequency centers: 30Hz/100 ms, 20 Hz/400 ms, and 10 Hz/700 ms.For 3-logons with noise, noise was added at the 4dB signal to noise level.In all simulations, each signal entered was assigned to a different simulated topographical region, to simulate the activity from different brain areas.To accomplish this, the signals were divided into 63 simulated channels creating a 7 × 9 grid within which differential weightings could be applied.Each signal entered, i.e. each logon, was weighted by a 4 × 4 grid differentially located within the overall 7 × 9 grid.The differential loadings were implemented to simulate a signal with more focal activity that decays in topographic space.The simulated datasets each contained 7560 total waveforms, comprised of 120 trials by 63 electrodes.

C. Results
In this section, we will present the results of applying the PCA-Gabor method on both simulated data and ERP signals described above.In the PCA-Gabor analysis, we will focus on extracting the 'best' logon fit to the principal component surface.Extracting the 'best' logon offers a way of parameterizing the PCs using the time, frequency and scale parameters of the logon as well as serving as a denoising tool since the 'best' logon will focus on representing the actual signal energy as opposed to background noise.Through this analysis, we will show the effectiveness of the PCA-Gabor method both as a modeling/data reduction tool and a denoising tool.The proposed method also offers an alternative to previous ERP studies that use matching pursuit to decompose each signal individually [13].This analysis has the disadvantage of being computationally expensive and extracting a large number of logons to represent a collection of signals.Since a comparison between matching pursuit at the individual signal level and the PCA-Gabor method would not be helpful due the number of logons extracted being much larger for the MP, we compare the PCA-Gabor method to the simultaneous matching pursuit with Gabor dictionary (SMP-Gabor) and to previous results obtained by the PCA method on the time-frequency plane [37].
1) Analysis of Simulated Data-The different methods were first evaluated for simulated data made up of Gabor logons.For this analysis, decompositions from two simulated datasets containing 3-logons with and without additive white noise were used.The PCA approach involved selecting the three time-frequency PCs with the highest eigenvalues.The PCA-Gabor approach extracted the 'best' Gabor logon for each of the three PCs yielding three logons.Finally, the SMP-Gabor method extracted the best three logons that explained the whole data set.For the 3-logons without noise all of the methods explained more than 95% variance of the data set with the PCA performing the best (Table I).The logons extracted from PCA-Gabor and SMP-Gabor were identical (see Fig. 1) explaining exactly the same amount of data variance indicating that under ideal conditions PCA-Gabor performs as well as the standard SMP-Gabor.
For 3-logons in 4dB noise, similarly, three components were extracted by each of the algorithms, i.e., three PCs with PCA, three logons fitted to the PCs with PCA-Gabor, and three logons fitted to the whole dataset with SMP-Gabor.The extracted components were evaluated in terms of the amount of signal variance they captured by projecting the components onto the original 3-logon dataset.From Table I, it can be seen that PCA-Gabor captures the most amount of signal variance with PCA coming in second.The logons extracted from the SMP-Gabor method can only explain 38.20% of the total signal variance since the algorithm focuses on extracting components that capture the most amount of common variance in the data (whereas PCA is covariance-based), which in this case corresponds to noise.Fig. 1 illustrates how the logons extracted by SMP-Gabor become wider in time and less correlated with the actual logons for the noisy data.This figure also shows that PCA-Gabor acts as a denoising mechanism reducing the noise in the PCs and thus representing more of the signal.
2) Analysis of Biological Data-For the biological data, first we will compare the variance characterized using the different approaches.We extract 11 PCs using PCA, the 11 logons extracted from these PCs using PCA-Gabor, and 11 logons that best explain the energy of the whole dataset using the SMP-Gabor method.Once the different components are extracted, they are projected onto each of the 8328 condition-averaged ERP waveforms.
For the three methods compared in this paper, PCA, PCA-Gabor and SMP-Gabor, PCA explained most of the data variance with 91%.For SMP-Gabor, the variance explained was 81%, whereas for PCA-Gabor it was 70%.
While the PCA explains the most overall variance in the data, and PCA-Gabor the least, additional analysis is needed to evaluate the methods in terms of experimentally relevant variance.To accomplish this, the three methods were compared for statistical separation using three common variables: sex, reaction time, and task difficulty.Because the activity extracted from the three methods covers much of the same time-frequency range, it is expected that the three methods should provide similar statistical effects.Statistical evaluation was conducted using a repeated-measures general linear model (GLM) including sex, reaction time, and task difficulty.A separate GLM was conducted for the sets of 11 components from each method.The design was Sex (male/female) by RT (reaction time; continuous) by task Difficulty (easy/hard; the within-subjects repeated measure).These main effects were highly significant for all three methods, confirming that similar experimentally relevant activity was extracted.Partial eta-squared values for the three methods are summarized in Table II.Here, for sex, RT, and difficulty, the nominal order of the amount of experimentally relevant variance in the statistical effects was the same, largest in the PCA-Gabor, next in the PCA, and the least in the SMP-Gabor.
By comparing the data reduction methods in terms of both overall data variance as well as experimentally relevant variance, stronger inferences can be made about how well the methods perform.In particular, the PCA-Gabor method captured the largest amount of experimentally relevant variance, while using the least amount of overall data variance.This is analogous to the results of the simulated data with noise, where the PCA-Gabor method extracted the most signal power in terms of experimentally relevant variance, while excluding the largest amount of noise power in terms of experimentally irrelevant variance.Thus, in these terms, the PCA-Gabor method was the most optimal among the three methods.
Finally, it is important to compare the different approaches in terms of computational complexity.All of the algorithms are run on a PC with Pentium 4 processor at 2 GHz using MATLAB 7.0, and evaluated after generating the time-frequency surfaces.The PCA-Gabor method took 11.6 seconds including the time to find the principal components (3.2 seconds) and to search for the best logon fit for the resulting PCs (8.4 seconds).The simultaneous matching pursuit on the other hand took 1276 seconds.Thus, in terms of computational complexity, PCA approach was the fastest, followed closely by the PCA-Gabor method.Trailing by a large margin is the SMP-Gabor method, which was computationally expensive due to the core search algorithm required.
Table III summarizes the key properties of the three methods compared in this section in terms of their data dependence, time and frequency parametrization, computational efficiency for the ERP data set and whether the resulting decomposition is based on explaining the most variance or covariance in the data.

D. Discussions
Several overall trends in the results are important to detail.First, the PCA-Gabor characterized more experimental variance than the PCA, with less of the overall raw data variance.This suggests that the Gabor decomposition of the PCA represents the relevant information obtained in the PCA, supporting the view that the activity extracted by PCA largely contains activity that conforms to Gabor constraints.Second, because the PCA-Gabor explains nominally more experimentally relevant variance, and outperforms the SMP-Gabor, while using less of the raw data variance than either, it supports the contention that this approach produces a more optimal Gabor decomposition of the collection of signals than the standard matching pursuit.Finally, it is interesting to specifically consider the fact that the PCA-Gabor method explains the least amount of data variance compared to the other two methods.The components extracted by PCA explain most of the data variance since PCA is designed to maximize the variance explained and extracts components that are orthogonal to each other.The PCA-Gabor method, on the other hand, approximates the energy of each principal component with a single logon and thus, the total variance explained is lower than the original principal components.However, this method has the advantage of retaining the signal variance and getting rid of the noise variance, thus acting as an effective denoising method, while also parameterizing the time-frequency surfaces.The third method, SMP-Gabor, explains more of the data variance compared to PCA-Gabor but has less experimental condition sensitivity (e.g., statistical significance).This increased variance and reduced sensitivity can be explained by looking at the Gabor logons extracted from the biological data by PCA-Gabor and SMP-Gabor methods shown in Fig. 2.Although the two methods extract some common logons, SMP-Gabor method emphasizes the low frequency activity.The first three logons extracted by this method are low frequency logons with a large time spread.The major reason for this is that the SMP-Gabor method operates entirely on the variance, and thus focuses the most on the high-amplitude (i.e., variance) low-frequency area of the surface.The PCA approaches, on the other hand, operate on covariance, which focuses more on activity that is functionally related (i.e.covaries).This point is also supported by the Gabor decomposition of the grand average of the 8328 waveforms given in Fig. 3. Table IV compares the parameters of the 11 logons extracted from the 11 PCs using the PCA-Gabor method, from the 8238 TFD surfaces using the SMP-Gabor method and from the grand average of the 8238 waveforms using standard MP-Gabor.This table indicates that there are some commonalities between the SMP-Gabor and MP-Gabor on the grand average surface since they extract similar logons describing the low frequency activity, e.g., logons 1-3 are almost identical.

V. Conclusions
In this paper, a time-frequency data reduction method combining a nonparametric datadriven approach, principal component analysis, with a parametric approach, matching pursuit with a Gabor dictionary, was presented.Using the proposed method, it was possible to characterize large amounts of ERP data with a small number of time-frequency parameters.This joint application of PCA with Gabor decomposition offered several advantages over individual PCA and Gabor decomposition.First, compared to PCA the proposed method improves the SNR of the extracted components, i.e., performs denoising, while simultaneously parameterizing the time-frequency surfaces and offering a succinct representation of the data set.Second, the application of Gabor decomposition onto the principal components instead of the actual data helps to extract parameters that represent the covariation among observations rather than characterize the average energy across observations.This property of PCA-Gabor becomes especially important when there is considerable noise in the data since standard matching pursuit algorithms will focus on fitting parameters to capture the most amount of energy, which in this case may be noise components.This phenomenon exhibited itself in the analysis of the biological data as PCA-Gabor most effectively differentiated between the experimental conditions with the least amount of data variance, or in other words capturing the least amount of noise, compared to the other two methods.This was mainly because the extracted logons explained the main effects described by the principal components with higher signal-to-noise ratio (most experimentally relevant variance).Future work will focus on the extensions of the proposed methods to different data factorization approaches such as the ICA.For ERP data collected over multiple channels, spatial ICA may be used as an alternative to PCA and the proposed data reduction method can be applied onto the independent components.Future work will also evaluate the Gabor parameters in relation to well-known cognitive ERP events such as P300, as well as ERP events with known specific neurological origins, such as anterior cingulate cortex activation as measured in the error-related negativity (ERN) paradigm.A comparison of the principal components analysis (PCA), Gabor logons extracted from the principal components (PCA-Gabor) and Gabor logons extracted by Simultaneous Matching Pursuit using a Gabor dictionary (SMP-Gabor) from two simulated datasets (3-logons with no noise and 3-logons with noise).A comparison of the principal components analysis (PCA), Gabor logons extracted from the principal components (PCA-Gabor) and Gabor logons extracted by Simultaneous Matching Pursuit using a Gabor dictionary (SMP-Gabor) for the ERP dataset.

4 .
Set k = k + 1, D ̃ = D ̃ \ g k , and go back to step 2 until a predetermined stopping criterion is achieved.The stopping criterion can either be a pre-selected number of atoms to describe the collection of signals or an average percentage of energy of the original signals described by the selected atoms.After M iterations, the following linear representation is obtained for each signal:(7)

3 . 4 .
Compute the covariance matrix, Σ = XX T .Decompose the covariance matrix using principal component analysis NIH-PA Author Manuscript NIH-PA Author Manuscript NIH-PA Author Manuscript(10)

Fig. 3 .
Fig. 3.11 Gabor logons extracted from the grand average of ERP signals.
scale parameter σ) of the 11 logons extracted by PCA-Gabor from the 11 PCs, SMP-Gabor from the 8238 TFD surfaces and MP-Gabor from the grand average of the 8238 time-