Sylvie Dufour

Laurent M Sachs

The chicken and egg situation of thyroid hormone and glucocorticoid signaling during postembryonic development

Keywords: Thyroid hormones, glucocorticoids, post-embryonic development, hormonal crosstalk

Post-embryonic development refers to the period of development after embryogenesis. This normally includes extensive growth, patterning and morphogenesis. In mammals including human, the perinatal period is marked by post-embryonic development regulated by a hormonal tsunami. Indeed, multiple hormones regulate fetal growth, brain development and organ maturation in preparation for extra-uterine life.

Two types of these hormones, thyroid hormones (THs) and glucocorticoids (GCs) play prominent and widespread roles during this period. They share common key-target regions (hypothalamus and pituitary) for the control of other hormone synthesis. They regulate common and specific processes, and interact at multiple levels to affect their targeted physiological processes. Their actions are mediated by ligand-activated transcription factors, which are members of the nuclear receptor superfamily and therefore have structural functional similarity. Genetic and environmental perturbations of TH and/or GC signaling may cause morphological and physiological impairments. These disruptions may also predispose individuals to diseases of adulthood, including neurologic, metabolic and cardiovascular diseases. All these significant health consequences motivate efforts to look for treatment options to avoid or mitigate potential pathologies. Unfortunately, endocrine and molecular mechanisms pertinent to the perinatal period are not well elucidated. In order to understand mechanisms in developmental endocrinology, animal models are required. Since such physiological processes are difficult to study in uterus-enclosed mammalian embryos, especially due to confounding maternal effects, alternative models can be attractive. The Anuran amphibian model is a pertinent choice. In this issue of Endocrinology, the manuscript from Shewade et al. [START_REF] Shewade | Corticosterone Is Essential for Survival Through Frog Metamorphosis[END_REF] presents elegant work that manages to address in amphibian the question of the interaction between TH and GCs signaling during post-embryonic development.

As for birth in mammals, larval metamorphosis which occurs in anuran amphibians and some fish, represents a remarkable post embryonic developmental crisis. Frogs and humans have THs and GCs which peak at metamorphosis and birth, respectively, have comparable roles for these hormones in many target organs, have conserved hormone receptors and mechanisms of gene regulation. In addition to highly conserved regulations, roles and action mechanisms of THs and GCs, frog is also a valuable vertebrate model because metamorphosis in their post-embryonic developmental period is independent of maternal endocrine influence. Amphibian metamorphosis is an ecological life transition accompanied by spectacular morphological, physiological and behavioral changes (such as diet, locomotor organs, respiratory system, etc.). Investigations by various authors revealed that the concerted responses to THs in different tissues towards sometimes opposite fates: apoptosis (e.g., resorption of the tail or the gills), proliferation (e.g., limb morphogenesis), and remodeling (e.g., the intestines or the central nervous system) require GCs [START_REF] Sachs | Insufficiency of Thyroid Hormone in Frog Metamorphosis and the Role of Glucocorticoids[END_REF]. In amphibians as in all vertebrates, GCs are also key regulators of intermediary metabolism, and the wellknown mediators of the stress response playing the role of interface with environmental signals. In Anuran amphibians, GCs increase THs signaling to accelerate metamorphic progression and allow escape from an unfriendly environment. Last but not least, metamorphosis corresponds to the transition from an aquatic environment to air breathing, that can be seen as a reminiscent of emergence onto land of tetrapod ancestors, and be also compared to the transition from foetal life in amniotic liquid to first air breathing of the newborn.

Since the past few years, TALEN or CRISPR-Cas gene knock out methodologies have been applied to decipher the specific roles played out by THs receptors, TR alpha and beta, in the regulation of the various metamorphic processes in Xenopus tropicalis [START_REF] Buchholz | Dual function model revised by thyroid hormone receptor knockout frogs[END_REF][START_REF] Sakane | Functional analysis of thyroid hormone receptor beta in Xenopus tropicalis founders using CRISPR-Cas[END_REF]. To address the requirement for GCs in metamorphosis, TALEN knock out of GCs receptor (GR) was recently performed, and revealed that GR is necessary for survival of Xenopus tropicalis throughout metamorphosis [START_REF] Sterner | Glucocorticoid receptor is required for survival through metamorphosis in the frog Xenopus tropicalis[END_REF]. In the work presented by the same group, Shewade et al. [START_REF] Shewade | Corticosterone Is Essential for Survival Through Frog Metamorphosis[END_REF] further investigated the role of GCs, by using TALEN gene disruption to knock out pomc (proopiomelanocortin) gene in Xenopus tropicalis. The pomc gene encodes a polypeptide hormone precursor, itself not functional, but that undergoes extensive tissue-specific posttranslational modification to generate a range of smaller biologically active peptides. In corticotroph cells of the anterior pituitary, the cleavage of the pro-peptide leads to one of the major end products, the adrenocorticotrophin (ACTH). ACTH is involved in the stimulation of the adrenal gland and in the cortical secretion of GCs. Shewade et al. made use of this property to prevent GCs production by knocking out the pomc gene. As expected, mutant tadpoles have a reduced plasma corticosterone level at metamorphosis. These tadpoles have also reduced expression of GCs-response genes and TH-response genes. Importantly, their rates of growth and development is reduced with death inevitably occurring towards the end of metamorphosis during tail resorption. The death of pomc KO mutants, like the GR KO mutants, further supports the critical role of GCs for survival through frog metamorphosis. Death can be prevented by exogenous corticosterone treatment but also by THs treatment. This ability of exogenous THs to overcome death in pomc mutants by itself indicates a strong interaction between the two hormones to ensure the completion of all transformations through metamorphosis. However, the physiological reasons of the death in GR or pomc KO mutants still remain to be addressed, as many hypotheses could be advanced, such as insufficient development of crucial organs (i.e. maturation of the lungs, a critical step at birth in human), dyssynchrony of organs/tissues transformation during metamorphosis, or shortage in metabolites and energy supplies highly essential through metamorphosis.

In human, as in frogs, POMC-deficiency results in hypo-cortisolism. At human birth, if untreated, it leads to death mainly as a consequence of lung immaturity. In premature infants, to remedy ineffective pulmonary function, prenatal GC administration is routinely used to accelerate lung maturation. Obesity, adrenal insufficiency and red hair are the phenotypes found when the pomc gene is affected in humans [START_REF] Krude | Severe early-onset obesity, adrenal insufficiency and red hair pigmentation caused by POMC mutations in humans[END_REF]. In mouse, homozygote deficiency of pomc gene is associated with obesity, increased food intake, adiposity, defects in adrenal development, and altered pigmentation [START_REF] Yaswen | Obesity in the mouse model of pro-opiomelanocortin deficiency responds to peripheral melanocortin[END_REF]. Serum corticosterone are below detectable levels. However, unlike frogs, lack of all derived POMC peptides is compatible with some survival throughout prenatal development, with still partial embryonic lethality associated with POMC deficiency. More interesting, the pomc -/-mice show lowered serum levels of thyroxine [START_REF] Challis | Mice lacking pro-opiomelanocortin are sensitive to high-fat feeding but respond normally to the acute anorectic effects of peptide-YY(3-36)[END_REF], implicating central melanocortin signaling in the regulation of the hypothalamicpituitary-thyroidal axis. Unfortunately, Shewade et al. do not contribute information on the issue of GCs actions on the hypothalamic-pituitary-thyroidal axis and THs production. Their surprising result showing that supplementation with THs allows the pomc knock-out tadpoles to survive suggests that GCs act on the state of thyroid axis or signaling. It seems that THs and GCs can regulate each other depending on molecular, cellular and physiological context. This cross-regulation can amplify hormone signals, regulate hormone activities, and coordinate hormone actions. In line with their crosstalk and synergism in the regulation of crucial steps of post-embryonic development, a concomitant rise in THs and GCs is observed during perinatal/post-embryonic period in mammals, hatching in birds and reptiles, as well as fish and amphibian metamorphosis. The cross-talk acts on different steps of the THs and GCs signaling pathways. At the level of the hypothalamus and pituitary, it involves corticotropin releasing hormone (CRH) itself, which in fish, amphibians, reptiles and birds, stimulates not only the release of ACTH as in all vertebrates, but also that of thyroidstimulating hormone (TSH), a regulation that might have been lost in mammals. At central and peripheral levels, when GCs synergize with THs to accelerate tadpole metamorphosis, GCs increase TR binding capacity, increase conversion of active THs from pro-hormone, and decrease the degradation of THs. GCs treatment also upregulate TR expression. Conversely, THs are involved in GCs receptor regulation, increasing or decreasing expression depending on the studied tissues. The situation is therefore similar to the chicken or the egg causality dilemma with the difficulty of sequencing actions where each seems to depend on others being done first. It seems that TH acts before GC and vice versa, that TH controls GC synthesis and vice versa, as well as TH action controls GC action and vice versa. All these interactions / modes of regulation deserve to be researched in mammals.

While focus has traditionally been put on the role of THs in amphibian metamorphosis, these new advances highlight the role of GCs and the cross-talk with THs. Involvement of GCs and THs cross-talk appears as a common feature to fish and amphibian metamorphosis, bird hatching and mammal birth. This suggests that THs and CGs interaction in the regulation of postembryonic developmental crisis, has an ancient origin and has been largely conserved throughout vertebrate evolution. From an evolutionary and adaptative perspective, the corticotropic stress axis with GCs, may be the best candidate to have first linked environmental cues to triggering of metamorphosis and change of ecophase, in vertebrates. On the other side, the major role of the thyrotropic axis in the regulation of development in vertebrates may lead to the opposite hypothesis. In extant species, the current knowledge does not allow either to clearly define the temporal sequence of THs and GCs actions and interactions in the regulation of post-embryonic development, a dilemma similar to the chicken and egg metaphor.

This being said, these observations come with a lot of new challenges and questions to fully understand the respective role of THs and GCs and their cross-talk during post-embryonic period, that can be addressed with relevant non-mammalian models, opening new research perspectives to potentially avoid or mitigate health consequences of disturbed THs or GCs signaling.

Acknowledgments

We thank Dr A. Virmani for proofreading the manuscript. Address all correspondence and requests for reprints to: Laurent M. Sachs, PhD, Unité Mixte de Recherche 7221 Physiologie Moléculaire et Adaptation, Département Adaptation du Vivant, Centre National de la Recherche Scientifique, Muséum National d'Histoire Naturelle, CP32, 57 Rue Cuvier, Paris, Cedex 05, 75231 France. E-mail: sachs@mnhn.fr.


Marco Congedo
email: marco.congedo@gmail.com

Cédric Gouy-Pailler

Christian Jutten

On the Blind Source Separation of Human Electroencephalogram by Approximate Joint Diagonalization of Second Order Statistics

Keywords: Blind Source Separation (BSS), Independent Component Analysis (ICA), Approximate Joint Diagonalization (AJD), Electroencephalography (EEG), Volume Conduction, Fourier Cospectra

On the blind source separation of human electroencephalogram by approximate joint diagonalization of second order statistics

Introduction

Recent studies on human electroencephalogram (EEG) are based on the theory of brain volume conduction. It is well established that the generators of brain electric fields recordable from the scalp are macroscopic post-synaptic potentials created by assemblies of pyramidal cells of the neocortex (Speckmann and Elger, 2005). Pyramidal cells are aligned and oriented perpendicularly to the cortical surface. Their synchrony is possible thanks to a dense net of local horizontal connections (mostly <1mm). At recording distances larger than about three/four times the diameter of the synchronized assemblies the resulting potential behaves as if it were produced by electric dipoles; all higher terms of the multipole expansion vanish and we obtain the often invoked dipole approximation (Lopes Da Silva and Van Rotterdam, 2005;Nunez and Srinivasan, 2006, Ch. 3). Three physical phenomena are important for the arguments we advocate in this study. First, unless dipoles are moving there is no appreciable delay in the scalp sensor measurement (Lopes da Silva and Van Rotterdam, 2005). Second, in brain electric fields there is no appreciable electro-magnetic coupling (magnetic induction) in the frequencies up to about 1MHz, thus the quasi-static approximation of Maxwell equations holds throughout the spectrum of interest (Nunez and Srinivasan, 2006, p. 535-540). Finally, for source oscillations below 40Hz it has been verified experimentally that capacitive effects are also negligible, implying that potential difference is in phase with the corresponding generator (Nunez and Srinivasan, 2006, p. 61). These phenomena strongly support the superposition principle, according to which the relation between neocortical dipolar fields and scalp potentials may be approximated by a system of linear equations [START_REF] Sarvas | Basic Mathematical and Electromagnetic Concepts of the Biomagnetic Inverse Problem[END_REF]. Whether this is a great simplification, we need to keep in mind that it does not hold true for all cerebral phenomena. Rather, it does at the macroscopic spatial scale we are interested in here.

A common approach to the study of human EEG is to describe patterns in space and time and link empirical findings with anatomical and physiological knowledge. The problem is characterized by high temporal resolution (about 1ms) and low spatial resolution (several cm 3 ). For example, it has been estimated that without time averaging about 60 million contiguous neurons must be synchronously active as to produce observable scalp potentials (Nunez and Srinivasan, 2006, p. 21). Such a cluster would realistically extend over several cm 2 of cortical gyral surface, whereas disentangling fields emitted by cortical functional units may require much higher precision. Because of volume conduction, scalp EEG potentials describe a mixture of the fields emitted by several dipoles extending over large cortical areas. Practically, in order to improve the spatial resolution it is often necessary to trade in the temporal one operating some form of temporal averaging. In summary, the path followed by much of current EEG research is to "isolate" in space and time the generators of the observed EEG as much as possible, counteracting the mixing caused by volume conduction and maximizing the signal-to-noise ratio (SNR).

Over the years we have assisted to the development of several classes of methods to improve the spatial specificity. Those include, among others, surface and cortical Laplacian [START_REF] Nunez | Electric Field of the Brain[END_REF], equivalent dipole fitting [START_REF] Mosher | Multiple dipole modeling and localization from spatio-temporal MEGdata[END_REF] and distributed minimum norm (model-driven) or minimum variance (data-driven) inverse solutions [START_REF] Greenblatt | Local Linear Estimators for the Bioelectromagnetic Inverse Problem[END_REF][START_REF] Da | Functional Localization of Brain Sources using EEG and/or MEG data: Volume Conductor and Source Models[END_REF]. Targeted attempts include sparsification approaches [START_REF] Gorodnitsky | Neuromagnetic source imaging with FOCUSS: a recursive weighted minimum norm algorithm[END_REF][START_REF] Cotter | Sparse solutions to linear inverse problems with multiple measurement vectors[END_REF] and spatial filters known as beamformers [START_REF] Rodríguez-Rivera | MEG and EEG Source Localization in Beamspace[END_REF][START_REF] Congedo | Subspace Projection Filters for Real-Time Brain Electromagnetic Imaging[END_REF]. Surface Laplacian methods apply a spatial high-pass filtering to the scalp potential by estimating their second spatial derivative. They tend to overemphasize high spatial frequency and radial (to the scalp surface) dipolar fields. Inverse solutions seek source localization in a chosen solution space and rely on geometrical models of the head tissue. Unfortunately, the accurate description of EEG volume conduction is complicated by inhomogeneity (resistivity varies with type of tissue) and anisotropy (resistivity varies in different directions); therefore source localization methods are inevitably undermined by geometrical modeling error.

Another approach that persists in EEG literature is blind source separation (BSS). First studied in our laboratory during the first half of the 80's (Ans et al., 1985;[START_REF] Hérault | Space or time adaptive signal processing by neural network models[END_REF] BSS has enjoyed considerable interest worldwide only a decade later, inspired by the seminal papers of Jutten and Hérault (1991), Comon (1994) and [START_REF] Bell | An Information-Maximization Approach to Blind Separation and Blind Deconvolution[END_REF]. BSS has today greatly expanded encompassing a wide range of engineering applications such as speech enhancement, image processing, geophysical data analysis, wireless communication and biological signal analysis [START_REF] Hyvärinen | Independent Component Analysis[END_REF][START_REF] Cichocki | Adaptive Blind Signal and Image Processing[END_REF][START_REF] Choi | Blind source separation and independent component analysis: A review[END_REF]. Such ubiquity springs from the "blind" nature of the BSS problem formulation: no knowledge of volume conduction or of source waveform is assumed. The problem may be attacked from several perspectives; several hundred BSS algorithms have been proposed over the last 20 years with more added on every year. Typically, such methods are based on the cancellation of second order statistics (SOS) and/or of higher (than two) order statistics (HOS). Their commonality resides in the assumption of a certain degree of source spatial independence, which is precisely modeled by the cancellation of those statistics. Both HOS and SOS have been employed with success in EEG. They are today established for denoising/artifact rejection [START_REF] Vigário | Extraction of ocular artifacts from EEG using independent component analysis[END_REF][START_REF] Jung | Removal of eye activity artifacts from visual event-related potentials in normal and clinical subjects[END_REF][START_REF] Vorobyov | Blind Noise Reduction for Multisensory Signals using ICA and Subspace Filtering, with Applications to EEG Analysis[END_REF][START_REF] Iriarte | Independent component analysis as a tool to eliminate artifacts in EEG: a quantitative study[END_REF][START_REF] Joyce | Automatic removal of eye movement and blink artifacts from EEG data using blind component separation[END_REF][START_REF] Kierkels | A model-based objective evaluation of eye movement correction in EEG recordings[END_REF][START_REF] Fitzgibbon | Removal of EEG noise and artifact using blind source separation[END_REF][START_REF] Frank | Automated protocol for evaluation of electromagnetic component separation (APECS): Application of a framework for evaluating statistical methods of blink extraction from multichannel EEG[END_REF][START_REF] Halder | Online artifact removal for braincomputer interfaces using support vector machines and blind source separation[END_REF][START_REF] Phlypo | Removing ocular movement artefacts by a joint smoothened subspace estimator[END_REF][START_REF] Romero | A comparative study of automatic techniques for ocular artifact reduction in spontaneous EEG signals based on clinical target variables: A simulation case[END_REF][START_REF] Crespo-Garcia | Muscle artifact removal from human sleep EEG by using independent component analysis[END_REF], improving brain computer interfaces (Qin et al., 2005;[START_REF] Serby | An improved P300-Brain-Computer Interface[END_REF][START_REF] Wang | Extracting Rhythmic Brain Activity for Brain-Computer Interfacing through Constrained Independent Component Analysis[END_REF][START_REF] Dat | Feature Selection Based on Fisher Ratio and Mutual Information Analyses for Robust Brain Computer Interface[END_REF][START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF] and for increasing the SNR of single-trial time-locked responses [START_REF] Cao | Independent component analysis for anaveraged single-trial MEG data decomposition and single-dipole localization[END_REF][START_REF] Sander | Single Evoked Somatosensory MEG Responses Extracted by Time Delayed Decorrelation[END_REF][START_REF] Lemm | Enhancing the Signal-to-Noise Ratio of ICA-Based Extracted ERPs[END_REF][START_REF] Tang | Contrasting single-trial ERPs between experimental manipulations: Improving differentiability by blind source separation[END_REF][START_REF] Guimaraes | Single-trial classification of MEG recordings[END_REF][START_REF] Zeman | Independent component analysis and clustering improve signal-to-noise ratio for statistical analysis of event-related potentials[END_REF]). Yet, it appears that only four of the many existing algorithms have repeatedly occurred in EEG literature. They are known as FastICA [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF], JADE [START_REF] Cardoso | Blind beamforming for non-Gaussian signals[END_REF], InfoMax (Bell and Bejnowsky, 1995) and SOBI [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF]. FastICA, InfoMax and JADE are ICA (HOS) methods, while SOBI is a SOS method. JADE and SOBI are solved by approximate joint diagonalization [START_REF] Cardoso | Blind beamforming for non-Gaussian signals[END_REF]Pham, 2001 b;[START_REF] Yeredor | Non-orthogonal joint diagonalization in the least-squares sensewith application in blind source separation[END_REF][START_REF] Ziehe | A Fast Algorithm for Joint Diagonalization with Non Orthogonal Transformations and its Application to Blind Source Separation[END_REF][START_REF] Vollgraf | Quadratic Optimization for Simultaneous Matrix Diagonalization[END_REF][START_REF] Li | Nonorthogonal Joint Diagonalization Free of Degenerate Solution[END_REF][START_REF] Fadaili | Nonorthogonal Joint Diagonalization/Zero Diagonalization for Source Separation Based on Time-Frequency Distributions[END_REF]Dégerine and Kane, 2007), a powerful algebraic tool which allows promising extensions that we will consider in this study.

The BSS problem for the brain

For N scalp sensors and M≤N EEG dipolar fields with fixed location and orientation in the analyzed time interval, the linear BSS model simply states the superposition principle discussed . On the other hand, biological noise (extra-cerebral artifacts such as eye movements and facial muscle contractions) and environmental noise (external electromagnetic interference) may obey a mixing process as well, thus they are generally modeled as components of ( ) t s , along with cerebral ones.

Notice that while biological and environmental noise can be identified as separated components of ( ) t s , hence removed, source estimation will be affected by the underlying cerebral background noise propagating with the same coefficients as the signal [START_REF] Belouchrani | Blind Source Separation Based on Time-Frequency Signal Representations[END_REF]. Second, the assumption of time-invariance of the mixing process in (1.0) must apply only locally. The demixing matrix is assumed fixed for a given temporal interval, but may be allowed to change (slowly) across successive intervals (Pham 2001 a;[START_REF] Li | Blind estimation of channel parameters and source components for EEG signals: a sparse factorization approach[END_REF]. Such a model allows changes in the location and orientation of dipole layers over time. The assumptions underlying model (1.0) are crucial for the success of the source separation, thus will be reconsidered in more details in the discussion.

A suitable class of solutions to the brain BSS problem

To tackle problem (1.1) assuming knowledge of sensor measurement only we need to reduce the number of admissible solutions. In this paper we are interested in weak restrictions converging

toward condition ˆ( ) ( ) t t = s Gs , (1.2)
where s(t) holds the time-course of the true (unknown) source processes and the system matrix

= ≈ G BA ΛP (1.3)
approximates a signed scaling (a diagonal matrix Λ) and raw permutation (P). Equation (1.2) is obtained substituting (1.0) in (1.1) ignoring the noise term in the former. Whether condition (1.2) may be satisfied is a problem of identifiability, which establish the theoretical ground of BSS theory [START_REF] Tong | Blind Estimation of correlated source signals[END_REF]Tong et al., 1991 a, b;[START_REF] Tong | Waveform-Preserving Blind Estimation of Multiple Independent Sources[END_REF][START_REF] Cardoso | Blind Signal Separation: Statistical Principles[END_REF][START_REF] Pham | Blind Separation of Instantaneous Mixtures of Non Stationary Sources[END_REF][START_REF] Pham | Exploiting source non stationary and coloration in blind source separation[END_REF][START_REF] Theis | Uniqueness of Complex and Multidimensional Independent Component Analysis[END_REF]. In turn, matching condition (1.2) implies that we can recover faithfully the source waveform out of a scale (including sign) and permutation (order) indeterminacy. The idea suits EEG well, since the waveform bears meaningful physiological and clinical information. Notice the correspondence between the m th source, its separating vector (m th row of B ) and its scalp spatial pattern (mixing vector), given by the m th column of ˆˆ+ = A B . Hereafter superscript + indicates the Moore-Penrose pseudo-inverse. The mono-dimensionality of those vectors and their sign/energy indeterminacy implies the explicit modeling of the orientation and localization parameters of the m th source, but not its moment. This is also the case of inverse solutions with good source localization performance [START_REF] Greenblatt | Local Linear Estimators for the Bioelectromagnetic Inverse Problem[END_REF]. Nonetheless, we can evaluate the relative energy of each source sorting them by decreasing value of explained variance, as illustrated in appendix (D). Moreover, linearity allows switching back from the source space into the sensor space. Substituting (1.1) into

(1.0) and dropping the noise term in the latter yields BSS filtering

ˆˆ' ( ) ( ) ( ) t t t = = v ARs ARBv , (1.4)
where R is a diagonal matrix with m th diagonal element equal to 1 if the m th component is to be retained and equal to 0 if it is to be removed.

Different approaches for solving the source separation problem

It has been known for a long time that in general the BSS problem cannot be solved for sources that are Gaussian, independent and identically distributed (iid) [START_REF] Darmois | Analyse générale des liaisons stochastiques[END_REF]. The iid condition implies that each sample of the source components is statistically independent from the others and that they all follow the same probability distribution. Therefore, in order to solve the BSS problem the sources must be either (1) possibly iid, but non-Gaussian or (2) not iid. In case (1), one assumes that at most one source is Gaussian and that they are all mutually statistically independent.

The mutual independence assumption (spatial independence of all pair-wise sources) should not be confused with the iid condition (temporal independence of successive samples within each source process). Actually, the iid condition implies that no temporal information is used, thus the method is efficient regardless the temporal dependence of sources. Those methods are known as independent component analysis (ICA) (Jutten and Hérault, 1991;Comon, 1994Comon, , 1999;;[START_REF] Hyvärinen | Independent Component Analysis[END_REF].

ICA requires higher order statistics (HOS), explaining why it may succeed only if at most one source has Gaussian distribution: in fact Gaussian distributions are fully defined by their statistics up to the second order (SOS). The idea of ( 2) is to break the non-Gaussianity assumption. This can be done by assuming that source components are all pair-wise uncorrelated and that either (a) within each source component the successive samples are temporally correlated1 , [START_REF] Tong | Blind Estimation of correlated source signals[END_REF][START_REF] Molgedey | HG Separation of a Mixture of Independent Signals using Time Delayed Correlations[END_REF][START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF][START_REF] Ziehe | TDSEP-an efficient algorithm for blind separation using time structure[END_REF] or (b) samples in successive time intervals do not have the same statistical distribution, i.e., they are non stationary [START_REF] Matsuoka | A neural net for blind separation of nonstationary signals[END_REF][START_REF] Souloumiac | Blind Source Detection and separation using second order nonstationarity[END_REF][START_REF] Choi | Blind Separation of nonstationary sources in noisy mixtures[END_REF][START_REF] Pham | Blind Separation of Instantaneous Mixtures of Non Stationary Sources[END_REF][START_REF] Choi | Second Order Nonstationary Source Separation[END_REF].

Provided that source components have non-proportional spectra or the time courses of their variance (energy) vary differently, one can show that second order statistics are sufficient for solving the source separation problem. Since second order statistics are sufficient, the method is able to separate also Gaussian sources, contrary to ICA.

In summary, ICA methods require a higher independence assumption (HOS independence), while SOS methods rely on a weaker uncorrelation assumption (SOS independence) coupled to the assumption that source components display unique spectral density signature (characteristic source coloration) and/or unique source energy variation signature (characteristic non-stationarity). We will see that different source energy in different experimental conditions is also a sufficient additional condition allowing separation using SOS. If these assumptions are fulfilled the separating matrix can be identified uniquely, thus source can be recovered regardless the true mixing process (uniform performance property: see for example [START_REF] Cardoso | Blind Signal Separation: Statistical Principles[END_REF] and regardless the iid condition for ICA or the Gaussian condition for SOS BSS.

A recent trend in the literature on EEG source separation is the design of algorithms blending spatial independence with other specific assumptions. Such an approach is called semi-blind source separation. Different priors have been introduced in the cost function pursuing spatial independence.

In a Bayesian framework [START_REF] Roberts | Independent component analysis: source assessment and separation, a Bayesian approach[END_REF] introduces priors on signal distributions, but those seems little useful in EEG since typically the distributions are not known a-priori. Temporal constraints have been introduced using ad hoc reference signals based on experimental stimulation or visual inspection of the observed mixtures [START_REF] James | Temporally constrained ICA: an application to artifact rejection in electromagnetic brain signal analysis[END_REF]Lu and Rajapakse, 2005). [START_REF] Zhang | Morphologically constrained ICA for extracting weak temporally correlated signals[END_REF] describes a method to extract reference signals from the observed waveforms designed to work even when the interesting source signal is not visible in the observed measurement. Spatial constraints have been employed to recover sources of interest with known spatial topographies [START_REF] Hesse | On semi-blind source separation using spatial constraints with applications in EEG analysis[END_REF][START_REF] Ille | Artifact correction of the ongoing EEG using spatial filters based on artifact and brain signal topographies[END_REF]. Spectral constraints have been also introduced to recover sources with known spectral content [START_REF] James | On the use of Spectrally Constrained ICA applied to single-channel Ictal EEG recordings within a Dynamical Embedding Framework[END_REF][START_REF] Wang | Extracting Rhythmic Brain Activity for Brain-Computer Interfacing through Constrained Independent Component Analysis[END_REF][START_REF] Barbati | Functional source separation applied to induced visual gamma activity[END_REF]. Finally, a flexible semi-blind ICA approach where one can incorporate priors of experimental or physiological origin as well has been proposed by [START_REF] Barbati | Functional source separation from magnetoencephalographic signals[END_REF]. All BSS methods cited so far rely on the spatial independence assumption. Other methods exist, relaxing definitely this assumption and imposing instead positivity of sources and mixtures (Lee andSeung, 1999, 2001) or source sparsity [START_REF] Gribonval | A survey of Sparse Component Analysis for blind source separation: principles, perspectives, and new challenges[END_REF][START_REF] Li | Blind estimation of channel parameters and source components for EEG signals: a sparse factorization approach[END_REF]. These methods will not be considered in this study despite their theoretical interest since so far their use with EEG has been marginal.

SOS vs. HOS: statistical considerations

Joyce et al. (2004) reports that successful separation of EEG data can be achieved using as few as 100 data points using SOBI (SOS) and 1000 or more using ICA (HOS) algorithms. This is a known advantage of SOS-based BSS methods; an higher statistical efficiency allows performing BSS on shorter time intervals, which is a safe strategy to prevent serious departures from the linear instantaneous model assumptions (see discussion). It has also been suggested that SOS estimations are more robust with respect to noise as compared to HOS estimations [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF][START_REF] Joyce | Automatic removal of eye movement and blink artifacts from EEG data using blind component separation[END_REF]. We notice that this is not necessarily true. For instance, the estimation of kurtosis is unaffected by white Gaussian noise, whereas this is not true for SOS estimations. However, SOS estimations are expected to be more robust with respect to outliers, in that their influence in ensemble average estimations is magnified by elevation to the second power, whereas for HOS it is magnified by elevation to the third and/or fourth power.

SOS vs. HOS: neurophysiological considerations

The hypothesis of spatial independence

Human neocortex is a prodigious net of local and global interconnections. There are about as many neurons (10 10 ) as cortico-cortical fibers connecting them in the 1-15 cm range (Nunez and Srinivasan, 2006, p. 7). Dense and sometimes distributed connections exist between the neocortex and sub-cortical structures as well. Therefore, one may ask if assuming independent time course of cortical cell assemblies is reasonable. It has been speculated that forcing independence of the BSS output may result in spurious source components [START_REF] Li | Blind estimation of channel parameters and source components for EEG signals: a sparse factorization approach[END_REF]. A study on induced visual gamma activity has also questioned the source non-Gaussianity assumption as required by ICA [START_REF] Barbati | Functional source separation applied to induced visual gamma activity[END_REF]. In practice, the BSS output is never exactly independent, but just as independent as possible and this may explain why BSS is useful with EEG. HOS and SOS BSS may be both conceived as spatial filters minimizing the dependence of the observed mixtures. For EEG data this is an effective way to counteract the effect of volume conduction. In fact, we have seen that the brain tissue behaves approximately as a linear conductor, thus observed potentials (mixtures) are more dependent than the generating dipolar fields. The fundamental difference between HOS and SOS BSS springs from the kind of statistical information they try to extract from the data. Thus, the critical question is how well the respective statistics inform about the actual source process, question to which we now turn.

Extra-cerebral Artifacts

An important problem with EEG recording is the contamination of extra-cerebral artifacts.

The The use of coloration has been suggested for separating EEG from eye movements [START_REF] Joyce | Automatic removal of eye movement and blink artifacts from EEG data using blind component separation[END_REF]. In most comparative studies [START_REF] Kierkels | A model-based objective evaluation of eye movement correction in EEG recordings[END_REF][START_REF] Fitzgibbon | Removal of EEG noise and artifact using blind source separation[END_REF][START_REF] Halder | Online artifact removal for braincomputer interfaces using support vector machines and blind source separation[END_REF][START_REF] Romero | A comparative study of automatic techniques for ocular artifact reduction in spontaneous EEG signals based on clinical target variables: A simulation case[END_REF] coloration have been found superior than non-Gaussianity for eye movements removal, but [START_REF] Phlypo | Removing ocular movement artefacts by a joint smoothened subspace estimator[END_REF] found inconclusive results (slightly favoring HOS) and Frank and Friskhoff (2007) concluded clearly in favor of HOS. We notice that the linear instantaneous model (1.0) does not admit rotating dipoles. Furthermore, the movements of the two eyes are extremely correlated, forming two spatially dependent dipoles. Insomuch, in principle linear BSS methods cannot resolve them. In effect, a method assuming instantaneous mixing will try to model an "average" big dipole covering both eyes, or in the middle of them and will try to explain the rotating orientation by two or three orthogonal components. This raises the problem of "misallocation of variance" and will result in sub-optimal separation [START_REF] Frank | Automated protocol for evaluation of electromagnetic component separation (APECS): Application of a framework for evaluating statistical methods of blink extraction from multichannel EEG[END_REF]. Based on these considerations we believe that a general BSS approach based on the spatial independence assumption is not optimal for separating EEG from eye movements, regardless the use of SOS or HOS. To accomplish this task semi-blind approaches appear more promising (e.g., [START_REF] Ille | Artifact correction of the ongoing EEG using spatial filters based on artifact and brain signal topographies[END_REF].

On the other hand, the separation of muscle contamination appears easier to treat. Generally both SOS and HOS methods perform correctly with them [START_REF] Crespo-Garcia | Muscle artifact removal from human sleep EEG by using independent component analysis[END_REF]), yet research on more targeted approaches is in progress [START_REF] Gasser | Correction of muscle artefacts in the EEG power spectrum[END_REF][START_REF] De Clercq | Canonical correlation analysis applied to remove muscle artifacts from the electroencephalogram[END_REF].

Spontaneous and induced EEG

Observed potentials are the summation of post-synaptic potentials over large cortical areas caused by trains of action potentials carried by afferent fibers. The action potentials come in trains/rest periods, resulting in sinusoidal oscillations of the scalp potentials, with negative shifts during the train discharges and positive shifts during rest. The periodicity of trains/rest periods are deemed responsible for high-amplitude EEG rhythms (oscillations) up to about 12Hz, whereas higher frequency (>12Hz) low-amplitude rhythms may result from sustained (tonic) afferent discharges (Speckmann and Elegr, 2005). There is no doubt that an important portion of spontaneous EEG activity is rhythmic, whence strongly colored (Niedermeyer, 2005 a;[START_REF] Steriade | Cellular Substrates of Brain Rhythms[END_REF]Buzsáki, 2006, Ch. 6, 7). Some rhythmic waves come in more or less short bursts. Typical examples are sleep spindles (7-14Hz) (Niedermeyer, 2005 b;[START_REF] Steriade | Cellular Substrates of Brain Rhythms[END_REF], frontal Theta (4-7Hz) and Beta (13-35Hz)

waves (Niedermeyer, 2005 a). Others are more sustained, as it is the case for slow Delta (1-2Hz)

waves during deep sleep stages III and IV (Niedermeyer, 2005 b), the Rolandic Mu rhythms (around 10Hz and 20Hz) and posterior Alpha rhythms (8-12Hz) (Niedermeyer, 2005 a). In all cases brain electric oscillations are not ever-lasting and one can always define time intervals when rhythmic activity is present and others when it is absent or substantially reduced. Such intervals may be precisely defined based on known reactivity properties of the rhythms. For example, in event-related synchronization/desynchronization (ERD/ERS: Pfurtscheller and Lopes da Silva, 2004), which are time locked, but not phase locked increases/decreases of the oscillating energy [START_REF] Steriade | Cellular Substrates of Brain Rhythms[END_REF], intervals may be defined before and after event onset. On the other hand event-related potentials (ERP: Lopes Da Silva, 2005 b), which are both time-locked and phase-locked can be further partitioned in several successive intervals comprising the different peaks. Such source energy variation signatures can be modeled precisely by SOS, as we will specify.

Transients

Another class of brain electric phenomena comprises transient waves such as spikes, sharp waves and spike-wave complexes in epileptic disorder (Niedermeyer, 2005 c), vertex waves during sleep (Niedermeyer, 2005 b), etc. Transients are characterized by abrupt and sometimes large potential shifts. In general, those are not naturally characterized by coloration, unless they results from the superposition of several continuous colored waves. This is the case, for example, of K-complexes observed during sleep (Niedermeyer, 2005 b), which are a superposition of a slow wave (<1Hz) and a

Delta wave (1-4Hz) [START_REF] Steriade | Cellular Substrates of Brain Rhythms[END_REF]. Nonetheless, transient activities are by definition spaced by intervals of inactivity, hence the difference between the energy in their active and inactive intervals (non stationarity) may be captured adequately by SOS statistics. However, due to their possible highly non-Gaussian nature, this kind of phenomena is naturally modeled by HOS statistics.

In summary, it appears that a wide variety of spontaneous and induced EEG phenomena are captured appropriately by SOS statistics, however for transient activity HOS may be better candidates.

So far SOS methods applied to EEG have concentrated mainly on coloration (e.g., SOBI). The validity of the coloration assumption for recovering actual EEG dipolar fields has received experimental support [START_REF] Tang | Validation of SOBI components from High-density EEG[END_REF][START_REF] Sutherland | Reliable detection of bilateral activation in human primary somatosensory cortex by unilateral median nerve stimulation[END_REF][START_REF] Van Der Loo | Correlation between Independent Components of scalp EEG and intra-cranial EEG (iEEG) time series[END_REF]. We have contended that source energy variation over time is a ubiquitous property of EEG and it should be exploited besides coloration. This is the focus of SOS time-frequency approaches, which are well established in other technical fields [START_REF] Belouchrani | Blind Source Separation Based on Time-Frequency Signal Representations[END_REF][START_REF] Pham | Exploiting source non stationary and coloration in blind source separation[END_REF][START_REF] Choi | Second Order Nonstationary Source Separation[END_REF][START_REF] Bousbia-Salah | A one step time-frequency blind identification[END_REF]. Here we pursue further this path in the context of experimental and clinical EEG data.

Approximate joint diagonalization

The class of SOS BSS methods we are considering is consistently solved by approximate joint diagonalization algorithms [START_REF] Cardoso | Blind beamforming for non-Gaussian signals[END_REF]Pham, 2001 b;[START_REF] Yeredor | Non-orthogonal joint diagonalization in the least-squares sensewith application in blind source separation[END_REF]Ziehe et al., 2004;[START_REF] Vollgraf | Quadratic Optimization for Simultaneous Matrix Diagonalization[END_REF][START_REF] Li | Nonorthogonal Joint Diagonalization Free of Degenerate Solution[END_REF][START_REF] Fadaili | Nonorthogonal Joint Diagonalization/Zero Diagonalization for Source Separation Based on Time-Frequency Distributions[END_REF]Dégerine and Kane, 2007). Given a set of matrices {Q 1, Q 2 ,… }, the AJD seeks a matrix B such that the products in {Q 1, Q 2 ,… } are chosen so as to hold in the off-diagonal entries statistics describing some form of dependence among the sensor measurement channels; then the AJD will vanish those terms resulting in linear combination vectors (the rows of B ) extracting "independent" components from the observed mixture via (1.1). More particularly, the joint diagonalization is applied on matrices that change according to the assumptions about the source. They are those changes, when available, that provide enough information to solve the BSS problem. If the source process is assumed to be colored, one may consider lagged covariance matrices of signals. If the source process is assumed to be non stationary between blocks of data, one may consider covariance matrices estimated on different time windows. In both situations, provided that source spectra are non proportional (colored sources) or source energy varies differently (non stationary sources), the additional matrices add information (in fact, equations) sufficient for estimating all the parameters of the separating system. If the source is both colored and non stationary, one can use a set of both kinds of matrices, as we will illustrate.

The aforementioned popular JADE and SOBI algorithms are based on AJD and this is the case for many other BSS algorithms (for a review see [START_REF] Theis | On the use of joint diagonalization in blind signal processing[END_REF]. One advantage of AJD algorithms is that they execute fast and do not require setting parameters for convergence. In particular, the algorithms by Pham (2001 b) and by [START_REF] Ziehe | A Fast Algorithm for Joint Diagonalization with Non Orthogonal Transformations and its Application to Blind Source Separation[END_REF] enjoy sustained popularity because of their good performance and computational efficiency. Like ICA algorithms, the AJD approach allows extracting source components by groups, which appears to us an effective way to overcome the aforementioned limitation of assuming pair-wise spatial independence of all EEG source processes; source components may now be assumed independent between groups but not necessarily independent within each group. Mathematically, this amounts to require the products 1 ˆˆT BQ B , 2 ˆˆT BQ B , … be block-diagonal instead of diagonal. Such an approach has been foreseen by Cardoso (1998 b) and is nowadays referred to as independent subspace analysis (ISA). Block-AJD (B-AJD) allows seeking brain networks (groups of dependent source processes) instead of just several disjoint "hot spot", which is in line with current trends in brain neurophysiology (e.g., [START_REF] Mantini | Electrophisiological signatures of resting state networks in the human brain[END_REF]. As per today B-AJD is limited in practice by the necessity of specifying a-priori the numerosity and composition of the groups [START_REF] Theis | Blind signal separation into groups of dependent signals using joint block diagonalization[END_REF][START_REF] Févotte | Pivot selection strategies in jacobi joint block-diagonalization[END_REF]. Research on AJD algorithms is currently flourishing. Recent trends include pursuing decomposition by blocks and seeking optimal weighting (e.g., [START_REF] Tichavskí | A Fast Approximate Joint Diagonalization Algorithm using a Criterion with a Block Diagonal Matrix[END_REF]. We believe that the resulting improvements hold promise for the BSS field and its applications to human electroencephalogram.

SOS BSS methods solved by approximate joint diagonalization

The first proposed SOS method [START_REF] Féty | New Methods for Signal Separation[END_REF]Tong et al., 1991 b) exploited signal coloration. It consisted on joint diagonalization of two matrices, the covariance matrix and a lagged covariance matrix, allowing an exact solution via the well-known generalized eigenvalue-eigenvector decomposition [START_REF] Choi | Second Order Nonstationary Source Separation[END_REF]Parra and Saida, 2003). The corresponding procedure for exploiting energy time variation traces back to the work of [START_REF] Souloumiac | Blind Source Detection and separation using second order nonstationarity[END_REF]; if the energy of a source component changes in two successive time intervals, then the component can be estimated by joint diagonalization of the two covariance matrices estimated on those intervals. Importantly, if the source is active in one interval and inactive in the other the obtained filter is optimal [START_REF] Souloumiac | Blind Source Detection and separation using second order nonstationarity[END_REF].

Along these lines see the discussion on super-efficiency in [START_REF] Pham | Blind Separation of Instantaneous Mixtures of Non Stationary Sources[END_REF]. Although very simple and fast, these two-matrix joint diagonalization methods are very sensitive to estimation errors of those matrices. If the noise covariance structure is different in the two matrices then the joint diagonalization of the signal structure is severely distorted. A considerable improvement is obtained by AJD of a larger set of matrices. This idea, first applied in SOBI for colored source components [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF], has been then applied to non stationary source components [START_REF] Choi | Blind Separation of nonstationary sources in noisy mixtures[END_REF][START_REF] Pham | Blind Separation of Instantaneous Mixtures of Non Stationary Sources[END_REF] and finally extended to both colored and non stationary source components [START_REF] Belouchrani | Blind Source Separation Based on Time-Frequency Signal Representations[END_REF][START_REF] Pham | Exploiting source non stationary and coloration in blind source separation[END_REF].

Practically, in many SOS methods such as SOBI the data are first whitened and normalized (sometimes it is said they are sphered) as ( )

( ) t t = z
Hv , where

M N ⋅ ∈ H ℝ
is such that the covariance of ( ) t z is the identity. Then, the AJD of a set of delayed covariance matrices (several lags: SOBI), and/or a set of covariance matrices on several windows of z(t) is performed. It is known that the pre-whitening may jeopardize the separation performance due to the estimation error of the data covariance matrix, which is exactly diagonalized at the expense of the other matrices [START_REF] Cardoso | On the performance of orthogonal source separation algorithms[END_REF][START_REF] Yeredor | Blind Separation of Gaussian Sources via Second-Order Statistics with Asymptotically Optimal Weighting[END_REF]Pham, 2001a). Hence, a better procedure is obtained by using robust whitening [START_REF] Choi | Second Order Nonstationary Source Separation[END_REF] or obtaining the AJD of a set of covariance matrices directly on v(t), which amount to avoiding the pre-whitening step altogether [START_REF] Ziehe | TDSEP-an efficient algorithm for blind separation using time structure[END_REF], or by diagonalizing partial autocorrelation matrices [START_REF] Dégerine | Second-Order Blind Separation of Sources Based on Canonical partial Innovations[END_REF]. As compared to the two-matrix diagonalization the AJD approach is known to be more robust and efficient [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF][START_REF] Belouchrani | Blind Source Separation Based on Time-Frequency Signal Representations[END_REF][START_REF] Choi | Second Order Nonstationary Source Separation[END_REF]. One problem encountered by researchers with SOBI is how to choose an appropriate set of lags [START_REF] Tang | Validation of SOBI components from High-density EEG[END_REF][START_REF] Tang | Recovery of correlated neuronal sources from EEG: The good and bad ways of using SOBI[END_REF]. For colored Gaussian auto-regressive (AR) processes the asymptotically optimal set of lags includes as many lags as necessary to describe the maximal order of the process (see for example [START_REF] Doron | Asymptotically Optimal Blind Separation of Parametric Gaussian Sources[END_REF]. The AR order/number of lags depending on several factors (e.g., sampling rate, number of peaks in the source power spectrum etc.), one should estimate it on data at hand. A simpler solution to this problem is treated in appendix (B). It arises after shifting the AJD problem into the time-frequency domain, framework that we now delineate.

Time-frequency expansions

Source separation methods can be applied in different representation spaces. In fact, applying to (1.0) any invertible and linearity-preserving transform T leads to

[ ] [ ] ( ) ( ) t t = v A s T T ,
which preserves the mixing model. Then, solving source separation in the transformed space still provides estimation of the matrix A or of its inverse B, which can be used directly in Eq. (1.1) for recovering the source s(t) in the initial space. For example, the transform T may be a discrete Fourier transform, a time-frequency transform such as the Wigner-Ville transform or a wavelet transform.

AJD-based SOS methods such as SOBI can be easily and conveniently transposed in the frequency domain, thence in the time-frequency domain, whether we perform the frequency expansion for several time segments. Such approach is currently attracting much interest in the BSS community, especially for audio and speech applications [START_REF] Belouchrani | Blind Source Separation Based on Time-Frequency Signal Representations[END_REF][START_REF] Choi | Second Order Nonstationary Source Separation[END_REF][START_REF] Bousbia-Salah | A one step time-frequency blind identification[END_REF][START_REF] Deville | Temporal and time-frequency correlation-based blind source separation methods[END_REF][START_REF] Zhang | Blind separation of nonstationary sources based on spatial time-frequency distributions[END_REF]Aïssa-El-Bay et al., 2007). There exist several time-frequency expansions. For its simplicity in this study we consider the short Fourier transform, from which Fourier cospectral matrices are readily estimated2 . We will compute Fourier cospectral matrices ( )

N N fi ⋅ ∈ C ℝ
for a range f : 1…F of discrete frequencies and for a range i:1…I of temporal windows. Temporal windows should be short enough to capture the energy variations over time and wide enough to allow satisfactory estimations of cospectral matrices for each of them separately. For each temporal window i the f th cospectral matrix ( ) fi C holds the portion of the sensor covariance matrix corresponding to the f th frequency. Its diagonal elements hold the power (autospectra) of each measurement channel while its off-diagonal elements hold the terms describing the in-phase SOS dependency for that time window and frequency. As we have seen those off-diagonal terms are canceled by AJD in order to recover uncorrelated source components. Clearly, cospectral matrices are affine to the delayed covariance matrices used by SOBI, since they are a linear transformation of each other (e.g., Bloomfield, 2000, p. 12;Pham, 2001 a). Nonetheless, working in the frequency domain is advantageous for several reasons: first, covariance statistical estimations in the time domain are distorted for temporally correlated processes like EEG [START_REF] Beran | Statistics for Long-Memory processes[END_REF]. Second, estimating cospectral matrices in the frequency domain is computationally more efficient than estimating delayed covariance matrices in the time domain 3 . Finally, the AJD of cospectra has been connected to the Gaussian mutual information criterion (Pham 2001[START_REF] Pham | Exploiting source non stationary and coloration in blind source separation[END_REF]. This places the ensuing method at the hearth of the BSS theory and steers toward the Cramér-Rao bound (Pham, 2001 a;[START_REF] Pham | Blind Separation of Instantaneous Mixtures of Non Stationary Sources[END_REF]. We are aware of only one study comparing the AJD of delayed covariance matrices (SOBI) to the AJD of cospectral matrices [START_REF] Doron | Asymptotically Optimal Blind Separation of Parametric Gaussian Sources[END_REF]) and it clearly showed the better performance of the latter.

Approximate joint diagonalization of cospectra (AJDC): an extended time-frequency approach

Without loss of generality, the AJDC solution to the BSS problem (1.1) can be written compactly such as ( ) If only one source component was involved the shading could be directly interpreted as color coded energy, but since in general many source components are considered we shall think at the shading as a coding for the dependency structure. The variations of the cospectral structure in the defined space along the dimensions are called signatures. We say that a source component has a characteristic signature if no other source component has the same signature. Successful separation of a source component is obtained if the diagonalization set describes a characteristic signature of it. In other words, the diagonalization set should include at least two matrices differing in the dependency structure of this source component (with respect to the others) and those changes must not be the same for any other source components. For example, if the source component is narrowband and its frequency range differs from the others (characteristic coloration), the cospectral structure of this as FFTW3 [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF]; in typical situations we may expect the computation complexity of Fourier cospectral matrices be 20 to 100 times smaller as compared to lagged covariance matrices.

ˆAJD = B C , ( 1
source component along the frequency dimension will change uniquely and this change will enable the identification of that component. Algorithms like SOBI seek those changes to recover source components having non proportional power spectrum. The advantage of the time-frequency approach is precisely that either coloration or non stationarity characteristic signature can be captured in the time-frequency plane and that either one suffices to achieve separation. Thus, the multidimensional approach is robust with respect to possible violations of each taken separately. An important aspect of data expansion is that it enhances the characterization of source signatures; while the noise power tends to spread uniformly in the time-frequency plane the source power will concentrate in characteristic regions, thus the method is more robust with respect to noise as well [START_REF] Belouchrani | Blind Source Separation Based on Time-Frequency Signal Representations[END_REF]. The same arguments can be strengthened profiting of further source diversities simultaneously, such as those of physiological and experimental origin. For instance, to separate the posterior Alpha rhythms from the Rolandic Mu rhythms one may use the fact that posterior Alpha rhythms, but not Mu rhythms, are blocked by eyes opening (Niedermeyer, 2005 a). , where the cospectra at F frequencies (f : 1…F ) are estimated for I time intervals (i : 1…I) and K experimental conditions (k : 1…K). With such a diagonalization set one would exploit the diversity of source energy between conditions in addition to generic coloration and time-varying energy; notice that in this case source components can be identified if their energy differs in at least two experimental conditions, regardless the uniqueness of their spectral and stationarity signatures (that is, even if the basic assumptions of the SOS BSS method do not hold), but also if their characteristic signature is across the frequency or time dimension but not across experimental conditions. Along the same line, we can exploit the reactivity of EEG oscillations as aforementioned discussed, the presence/absence of a steady-state sensory stimulation, the presence/absence of electrical or magnetic stimulation, etc.; one may add as many indexes as desired and always proceed by (1.6).

We have seen that adding dimensions for expanding the data increases the chance to uncover the characteristic signatures of source components and increase the robustness with respect to noise.

However the number of matrices in the diagonalization set cannot be increased indefinitely. The essence of AJD algorithms consists in approximating the "average eigen-structure" of the input matrices. In general, any set (1.7) can be exactly jointly diagonalized if the instantaneous linear model holds exactly (Hyvärinen et al., 2001, p. 344), in which case all matrices in the set share common eigenvectors and the two subspaces spanned by those eigenvectors and the columns of mixing matrix A are identical [START_REF] Belouchrani | Blind Source Separation Based on Time-Frequency Signal Representations[END_REF]. This in practice will not quiet happen because of sampling estimation errors and noise, and while the latter is reduced the former is increased by data expansion, making more difficult finding the average eigen-structure. Another drawback of multiple dimension data expansion is that the instantaneous linear model (1.0) may not hold for all dimensions.

Finally, an open question is how the time-frequency plane should be sampled (on the other hand sampling of experimental conditions is given by definition). We see that the proper definition of the diagonalization set is the very challenge of AJD-based algorithms. A useful tool to identify regions where the characteristic signatures reside is described in appendix (B).

Discussion

Blind source separation (BSS) is a widespread method used in a number of scientific and technical fields [START_REF] Hyvärinen | Independent Component Analysis[END_REF][START_REF] Cichocki | Adaptive Blind Signal and Image Processing[END_REF]. Its use in EEG literature is currently growing at a fast pace. When applied to EEG data BSS decomposes scalp signals in a number of components. These components may correspond to the activity of cortical dipole layers generating the observed EEG. Precisely, BSS implicitly estimates their orientation and explicitly estimates their waveform (out of a sign and energy arbitrariness) and mixing coefficients. From the latter the spatial location can be estimated using an inverse solution method (Lopes da Silva, 2004; [START_REF] Greenblatt | Local Linear Estimators for the Bioelectromagnetic Inverse Problem[END_REF]; see for example [START_REF] Van Der Loo | Correlation between Independent Components of scalp EEG and intra-cranial EEG (iEEG) time series[END_REF]. Environmental and physiological artifacts may be extracted as well, while effective reduction of background noise may require additional noise suppression procedures (e.g., [START_REF] Vorobyov | Blind Noise Reduction for Multisensory Signals using ICA and Subspace Filtering, with Applications to EEG Analysis[END_REF]. Typically, BSS estimations feature higher SNR and strong suppression of the interference generated by other dipoles as compared to raw EEG. When the assumptions hold, BSS provides optimal estimations, in that they do not depend on physical modeling of the head. When the assumptions do not hold, BSS provides another representation of the sensor measurement space, which may still be useful (e.g., artifact reduction by BSS filtering as per Eq. 1.4), but that may encourage misleading interpretations of source waveform and associated topographies. One should be careful in claiming that the extracted components and associated topographies correspond to actual EEG physiological sources. Since checking the assumptions of the chosen BSS model and method is virtually impossible, the most credible arguments are those founded upon neurophysiological knowledge.

The linear BSS instantaneous model (1.0) makes a number of restrictive assumptions that are rarely checked or investigated. One assumption is that the number of sources is not greater than the number of sensors. When this is not the case (undetermined case) it is not possible to solve the BSS problem unless other constrains on the sources are introduced (e.g., sparsity: [START_REF] Gribonval | A survey of Sparse Component Analysis for blind source separation: principles, perspectives, and new challenges[END_REF]. One may also wonder if during the analyzed time interval the number of active dipoles is stable [START_REF] Li | Blind estimation of channel parameters and source components for EEG signals: a sparse factorization approach[END_REF]. In practice, brain electrical "source components" are macroscopic electric dipole with relatively high SNR formed by the synchronous activity of pyramidal cells over large cortical areas [START_REF] Nunez | Electric Field of the Brain[END_REF]. For sufficiently small time intervals one may assume that such high-SNR layers are limited in number. Other concurrently active cortical columns may be ignored if their current is comparatively negligible and it does not matter if the dipoles are active throughout the time interval or intermittently (actually such non stationarity signature can be explicitly exploited). Henceforth, assuming at least as many sensors as relevant sources does not appear problematic if we consider a sufficiently small time interval. Still, no definitive solution exists to the problem of estimating the number of source components in the overdetermined case (more sensor than source components). Whereas correct dimensionality reduction (appendix E) allows exact determination, over-reduction must be avoided since in this case identifiability is lost and several generators are extracted mixed in one component. A safe strategy is to identify a few meaningful components and keep reducing the dimension until those components are not distorted.

Another assumption is that the mixing matrix A in (1.0) is full-column rank. The columns of A are scalp spatial pattern vectors of the source components and the more the electrodes are close to each other, the more those vectors will be collinear. Consequently, it is always better to space the electrodes as much as possible on the scalp4 . Several restrictive assumptions are made by model (1.0) also on the nature of brain electric fields. One may ask whether it is reasonable to assume that dipoles keep fixed orientation and location in the analyzed time interval. Each row vector of the matrix B can be conceived as a spatial filter extracting the electric field of a dipole with a given fixed spatial extension, location and orientation. For a fixed spatial sensor configuration with respect to the brain, which is the case of a single EEG recording session, the orientation and location of electric dipoles are fixed by the anatomy and physiology of the grey matter forming the dipole. However, the dipole approximation becomes untenable for sources distributed over large areas [START_REF] Malmivuo | Bioelectromagnetism. Principles and Applications of Bioelectric and Biomagnetic Fields[END_REF][START_REF] Nunez | Electric Field of the Brain[END_REF]. Also, there is convincing evidence of traveling waves phenomena in the brain; long wavelength waves originating in a region and propagating via cortico-cortical connections to other regions (Lopes da Silva and Van Rotterdam, 2005;[START_REF] Srinivasan | Source analysis of EEG oscillations using high-resolution EEG and MEG[END_REF][START_REF] Thorpe | Identification of wave-like spatial structure in the SSVEP: Comparison of simultaneous EEG and MEG[END_REF]. These phenomena cannot be modeled by an instantaneous model and become more equivocal with larger time intervals. Also, the longer the time interval under analysis the less tenable is the stationarity assumption, which is basic to SOS estimations (Hyvärinen et al., 2001, p. 49). At the same time one must care to retain enough data points for analysis in order to avoid overfitting [START_REF] Müller | Blind source separation techniques for decomposing eventrelated brain signals[END_REF]. [START_REF] Särelä | Overlearning in Marginal Distribution-Based ICA: Analysis and Solutions[END_REF] reported that using small time intervals the output may contain artefacts that are not present in the data. For HOS method such as FastICA artifacts takes the form of artificial spikes and bumps, whereas for SOS methods such as SOBI they take the form of artificial sinusoid waves. [START_REF] Meinecke | A resampling approach to estimate the stability of one-or multidimensional independent components[END_REF] and [START_REF] Müller | Blind source separation techniques for decomposing eventrelated brain signals[END_REF] addressed the problem of obtaining robust and reliable source estimates. They proposed a resampling-based methods consisting in running the algorithms on different time intervals and retain only the source processes that can be found consistently. In conclusion, although statistical estimations improve with the number of samples we advocate the use of multiple time intervals as short as possible (enough to avoid overfitting while justifying the BSS method assumptions), modeling appropriately the stationarity within intervals while exploiting explicitly the non stationarity between intervals. In this sense an efficient time-frequency approach appears a precious option. Although we have contended that SOS-AJD methods such as AJDC fit well EEG data in general, a safe strategy is to compare the output to at least one HOS methods with any real EEG data problem at hand. We also notice that in the EEG field the instantaneous model has been rarely challenged [START_REF] Anemüller | Complex independent component analysis of frequency-domain electroencephalographic data[END_REF]Dyrholm et al., 2005). It is unfortunate that throughout comparisons of linear instantaneous, time-varying and convolutive model are lacking since the latter two families of BSS models may admit moving dipoles and traveling waves.

In this study we have described a simple time-frequency approach based on the approximate joint diagonalization of Fourier cospectral matrices (AJDC). AJDC is an extension of popular AJDbased algorithms such as SOBI, which are derived as restricted instances (exploiting source coloration or source non stationarity only), yet it is efficient statistically and computationally. Computationally, the AJDC equivalent of SOBI is several tens of times faster than SOBI. In turn, SOBI is known to be faster than, in the order, JADE, FastICA and InfoMax, the latter being the slowest [START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF]. Although those authors do not quantify precisely the complexity of each algorithm, we can safely say that AJDC is tens of times faster than SOBI and hundreds to thousands times faster than FastICA and InfoMax. Moreover AJDC (as all AJD-based algorithms) does not require parameter tuning for convergence. However, it requires an appropriate definition of the diagonalization set to correctly identifying the potential diversities in the data set. Instead of understanding this as a nuisance, we have contended that it amounts to correctly identifying the relevant aspects of the data variance at hand. Such an "informed" approach is somehow in between the completely blind setting, in which no a-priori knowledge on the source is assumed and the semi-blind approach, where temporal, spatial, spectral or other constraints are introduced in the cost function [START_REF] Roberts | Independent component analysis: source assessment and separation, a Bayesian approach[END_REF][START_REF] Ille | Artifact correction of the ongoing EEG using spatial filters based on artifact and brain signal topographies[END_REF]James and[START_REF] James | Temporally constrained ICA: an application to artifact rejection in electromagnetic brain signal analysis[END_REF][START_REF] James | On the use of Spectrally Constrained ICA applied to single-channel Ictal EEG recordings within a Dynamical Embedding Framework[END_REF]Lu and Rajapakse, 2005;[START_REF] Hesse | On semi-blind source separation using spatial constraints with applications in EEG analysis[END_REF][START_REF] Barbati | Functional source separation from magnetoencephalographic signals[END_REF][START_REF] Wang | Extracting Rhythmic Brain Activity for Brain-Computer Interfacing through Constrained Independent Component Analysis[END_REF][START_REF] Barbati | Functional source separation applied to induced visual gamma activity[END_REF][START_REF] Zhang | Morphologically constrained ICA for extracting weak temporally correlated signals[END_REF]. The basic time-frequency approach exploits the temporal dependency and energy variation over time of EEG. The diagonalization scheme can be defined so as to maximize the chance of separating dipole layers responsible for brain functions studied by experimental manipulation.

Assessing the difference in two or more experimental conditions is customary in cognitive and clinical studies using either continuous recording or evoked potentials paradigms. In this sense, AJDC may be an ideal companion for a very wide range of EEG experimental research. In the appendix we have collected several useful details about the effective use of AJDC for EEG data, which cannot be found elsewhere. Those details may be valuable to the reader interested in implementing AJDC or other time-frequency BSS algorithms. Code for AJD algorithms is commonly publicly available. An executable application performing BSS by AJDC can be obtained upon request to the corresponding author.

Appendix

A. Fourier cospectral matrices

For real data, if L is the length of the time window and S the sampling rate expressed in Hz, there are F=L/2+1 Fourier frequencies with resolution r=S/L equally spacing the range from 0Hz

(DC-level) to the folding frequency (0Hz, 1r Hz, 2r Hz,…,L/2r Hz). Typically, we take both S and L as a power of two and cospectral estimates may be averaged within arbitrary time intervals by sliding overlapping windows. The latter strategy allows arbitrary time intervals length. The Fourier cospectra and quadrature spectra are defined as the real and imaginary part of the Fourier cross-spectra [START_REF] Bloomfield | Fourier Analysis of Time Series[END_REF]. They are estimations of, respectively, the in-phase (or with a half cycle phase shift, i.e., opposite sign) and out-of-phase (a quarter cycle in either direction) covariance structure at frequency f. The discrete Fourier transform of sampled time-series x (t) over an epoch of length L is given by ( ) [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF]. Here below is the formula for computing the cospectral matrix at frequency f for time-series x (t) and y (t) :

1 2 1 0 ( ) L ift x L t d f x t e π - - = = ∑ . Let (
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) x x x x x y x y f y x y x y y y y d f d f d f d f d f d f d f d f d f d f d f d f d f d f d f d f ℜ ℜ ℑ ℑ ℜ ℜ ℑ ℑ ℜ ℜ ℑ ℑ ℜ ℜ ℑ ℑ   + + =     + +   C .
The formula readily extends to any N-dimensional input time-series to obtain its N-dimensional cospectral matrix. Notice that the cospectral matrix is symmetric and that the diagonal elements are the auto-spectra, better knows as power spectra. For an arbitrary long EEG segment we typically obtain an estimate of the cospectral matrix averaging C (f) over overlapping epochs of length L [START_REF] Welch | The Use of Fast Fourier Transform for the Estimaton of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms[END_REF]. Such estimates may then be summed across adjacent frequencies to obtain estimates within band-pass regions of interest. Summing all of them yields the covariance matrix as per Parseval's theorem. Pham (2001 a) shows that its Gaussian mutual information approximate joint diagonalization (AJD) criterion is the same using (complex) cross-spectral matrices or their real part (cospectral matrices). More in general, estimation of the separating matrix B can be obtained by AJD of cospectral matrix, quadrature spectra matrices or cross-spectral matrices if the mixing matrix is real [START_REF] Theis | Uniqueness of Complex and Multidimensional Independent Component Analysis[END_REF], which we always assume as per (1.0) due to the absence of capacitive effects in the brain (see introduction).

B. Defining diagonalization sets

Some AJD algorithms allow weighting the diagonalization effort across the input matrices (e.g., Pham, 2001 b). Here we seek an adaptive solution to the weighting problem. Let us consider a set w of non-negative weights to be associated to each matrix of the diagonalization set. For consistency, we take w so that the average of its elements equals 1.0. As a general strategy, we may encourage the mask to be sparse (many zero entries), so to enable data expansion in multiple dimensions while keeping the diagonalization set of reasonable size. Because of volume conduction dipole fields result in covariance structures with many non-null off-diagonal terms (see for example simulated forward solutions in [START_REF] Congedo | Subspace Projection Filters for Real-Time Brain Electromagnetic Imaging[END_REF]. On the other hand, higher frequency cospectra tend to be near diagonal because EEG energy decreases with frequency while spatially uncorrelated noise features a diagonal covariance structure. [START_REF] Whitham | Scalp electrical recording during paralysis: quantitative evidence that EEG frequencies above 20 Hz are contaminated by EMG[END_REF] even suggest that scalp EEG recording above 20Hz contains mainly electromyographic (EMG) activity. EMG exhibits near-diagonal covariance structures because it is spatially focal and does not propagate easily (through skull) to other leads. Consequently, matrices close to diagonal form should be down-weighted. In general, down-weighting low signal-to-noise ratio matrices amounts to effective noise suppression for the estimation of the separating matrix (Pham, 2001 a). Therefore, let us define for any cospectra

C (υ) 2 ( ) ( ) 2 ( ) 1 ( ) 1 rc r c rc r c C N C υ υ υ δ ≠ = = - ∑ ∑ C , (B.1)
where ( )rc C υ is the entry of matrix ( ) υ C at row r and column c and N is the size of the matrix (number of channels). For a positive definite matrix, measure (B.1) is bounded inferiorly by zero, for a diagonal matrix, and superiorly by 1.0, for a uniform matrix. Equation (B.1) provides a suitable nondiagonality weighting function: the higher the non-diagonality the higher the weight. Sparsification (noise-suppression) may be promoted by zeroing the weights above a cut-off frequency. According to our experience, such a weighting function generally allows satisfactory source estimation with continuously recorded EEG. We have observed that the non-diagonality function (B.1) is highly correlated with overall energy (trace of the cospectral matrices), but is not as much influenced by the dominant occipital rhythms (8-13Hz). In this fashion using a non-diagonality weighting function is in line with previous works in time-frequency BSS where the diagonalization effort has been concentrated on high-energy time-frequency regions [START_REF] Belouchrani | Blind Source Separation Based on Time-Frequency Signal Representations[END_REF].

C. Removing the DC-level (assuming zero-mean processes)

Typically, BSS models assume zero-mean processes. Thus, for DC EEG amplifiers, the DClevel needs to be removed. Simply, the first cospectrum (0Hz) is not considered in the diagonalization set. Notice that FFT estimates at positive frequencies are not affected by the DC level (Bloomfield, 2000, p. 90), hence there is no need to remove the mean, detrend or band-pass the signal before computing the FFT (the same is not true for lagged covariance matrices).

D. Evaluating the explained variance of source components

In the introduction we have suggested that the energy of the output source can be evaluated in spite of their sign arbitrariness. For simplicity, we illustrate the method for diagonalization sets of the

∑

. Notice that the explained variance or its relative portion VAR m /VAR TOT can be evaluated for any discrete frequency using cospectral matrix C (f) instead of V in computing both VAR TOT and VAR m above. In the same way, one may evaluate the explained variance for any frequency band pass region using instead the sum of cospectra within the region. This turns useful when we need to evaluate the energy of several components describing brain oscillations in a specific frequency band; for example, the several rhythms usually observed in the Alpha (8-12 Hz) range.

E. Subspace reduction and pre-whitening

When using many EEG sensors it may be useful to estimate fewer source components than sensors (M<N). Reducing the dimension of the input matrices makes them better conditioned, which enhance the performance of the separation [START_REF] Pham | Blind Separation of Instantaneous Mixtures of Non Stationary Sources[END_REF] 6 . The subspace reduction may follow different strategies. For instance, one may use model-driven beamforming to attenuate the signal originating outside the region of interest [START_REF] Rodríguez-Rivera | MEG and EEG Source Localization in Beamspace[END_REF][START_REF] Congedo | Subspace Projection Filters for Real-Time Brain Electromagnetic Imaging[END_REF]. Here we show how to perform subspace reduction to estimate the M most energetic source components, which is an extension of the common pre-whitening step. Let Note that it is not necessary to constrain the AJD matrix E to be orthogonal, effectively circumventing the aforementioned drawback of pre-whitening the data. Finally, note that working with AJDC we do not need to compute the covariance matrix of the data at all.

  most common artifacts are electric signals produced by the eyes and muscles of the face, jaw and neck. Artifacts have characteristic spatial, coloration and non stationary signatures, well distinct from EEG spontaneous activity[START_REF] Ille | Artifact correction of the ongoing EEG using spatial filters based on artifact and brain signal topographies[END_REF][START_REF] Da | Event Related Potentials: Methodology and Quantification[END_REF]. Still, their separation and removal has remained a difficult task. Eye blinks are attributed to change of conductance due to the moving eyelid on the cornea. They generate a potential peak which amplitude can be one order of magnitude stronger than the EEG. Eye movements, both volitional movements and saccades, are attributed to the cornea-retina dipolar field. Typically, vertical and horizontal movements are monitored by means of electrodes positioned above and next to the eye, providing an orthogonal space for movements in all possible direction. On the other hand, contamination by facial, jaw and neck muscles typically manifests as persistent low-voltage high-frequency signals (>20Hz:[START_REF] Whitham | Scalp electrical recording during paralysis: quantitative evidence that EEG frequencies above 20 Hz are contaminated by EMG[END_REF] with focused spatial distribution and rapidly decaying autocorrelation function.

  … are as diagonal as possible (subscript "T" indicates matrix transposition) . Given an appropriate choice of the diagonalization set {Q 1, Q 2 ,… } such matrix B is indeed an estimation of the separating matrix in (1.1) and one obtain an estimate of the mixing matrix as ˆˆ+ = A B . Matrices

  .6) where C : {C 1 , C 2 ,…} is the diagonalization set, i.e., a set of estimated Fourier cospectral matrices to be simultaneously diagonalized. The rational behind AJDC is expressed schematically in Fig 1. Each cube of the parallelepiped in the figure represents abstractly a cospectral matrix. The grid of cubes represents the sampling of some source property unfolding along two continuous dimensions (time and frequency) and one discrete dimension (experimental conditions). The different pattern of shading in each cube represents the different cospectral structure of each sampled region of the defined space.

  Two time intervals separated by eyes opening should then be considered. To exploit possible source energy diversity in several experimental conditions it suffices to average cospectral matrices separately within each condition, as indicated schematically in Fig. 1. This allows much flexibility, for an arbitrary number of cospectra computed on short time intervals can be averaged for each condition. To visualize the comprehensive nature of the method one may imagine the parallelepiped in Fig 1 in any number of dimensions susceptible to describe variations in some source statistical property. Putting all this in mathematical formalism turns out simple and elegant. Without loss of generality we shall always proceed by (1.6) after defining just a container for an arbitrary number of indexes and where each index indicates the sampling along a dimension. For example, the diagonalization set of Fig. 1 is obtained by defining fik υ ≜

Figure 1 .

 1 Figure 1. Schematic representation of the extended time-frequency blind source separation approach based on approximate joint diagonalization. See text for details.

  efficiently estimated by fast Fourier Transform (FFT: Cooley and Tukey, 1965; Frigo

  , when only source coloration is exploited. The method readily extends to any number of indices. First, let us scale the rows of estimated separating matrix B so that they all have unit L2 norm. Because of the energy arbitrariness this operation does not alter the output of the BSS. Let b m T and a m be, respectively, the normalized m th row of B (separating vector) and the m th column of ˆˆ+ = A B (spatial pattern) associated with the m th component. Let matrix of the raw EEG data. Its diagonal elements V nn hold the variance (energy) of the n th EEG channel. Using (1.1) and (1.0) and ignoring the noise term in the latter, the total explained variance of the source components is given by if M=N, since then ˆˆ= AB I . Similarly, the explained variance of the m th component alone is such as

  the original (unreduced) diagonalization set. As in Eq. (1.7), ν is a holder for a number of indeces. Now find a matrix (noise subspace). Note that for the casef υ ≜ , TOT Cis the sum of cospectral matrices at several frequencies and, if all Fourier frequencies are included in the diagonalization set, then matrix H is the well-known whitening matrix. For ... ijk υ ≜ our definition of TOT C is the natural extension to obtain a "global" whitening matrix. Let us now factorize the separating matrix such as ˆ= B EF , with ˆM M ⋅ ∈ E ℝ . We obtain a new diagonalization set D by applying the reduction to all cospectral matrices such as

  is called the demixing or separating matrix. Hereafter the caret indicates a statistical estimation. Although this is the classical BSS model we need a few clarifications for the EEG case: first, by ( ) t η we model instrumental noise only. In the following we drop the ( ) t η term because the instrumental (and quantization) noise of modern EEG equipment is typically low (<1µV

		v	( ) t	=	As	( ) ( ) t t + η ,	(1.0)
	where ( ) t ∈ v	N ℝ is the sensor measurement vector,	N M ⋅ A ℝ is a time-invariant full column rank ∈
	mixing matrix, ( ) t ∈ s	M ℝ holds the time-course of the source components and ( ) t ∈ η	N ℝ is additive
	noise, temporally white, possibly uncorrelated to ( ) t s and with spatially uncorrelated components.
	Our source estimation is given by
		s	t	)	=	( ( ) t Bv ,	(1.1)
	where	∈ B ℝ	M N ⋅
	above, i.e.,

Such processes are called colored, in opposition to iid processes, which are called white.

See appendix (A) for details on Fourier co-spectral matrices.

We have analyzed the computational complexity of estimating the former and latter matrices. Fourier cospectra estimations may take advantage of efficient split-radix fast Fourier transform (FFT) algorithms such

This suggests that placing many electrodes closely spaced above the brain region of interest, as it is sometimes done, is not a convenient strategy if multivariate statistical methods are to be employed.

For the first (0Hz) and last (L/2κHz) Fourier frequency the coefficients are real.

When the analyzed time interval is short the cospectral matrices may be non positive definite, a requirement of Pham's AJD algorithm. In this case the subspace reduction is necessary to obtain convergence. The algorithm of[START_REF] Ziehe | A Fast Algorithm for Joint Diagonalization with Non Orthogonal Transformations and its Application to Blind Source Separation[END_REF] does not impose this restriction but does not allow explicit weighting.

This Research has been partially supported by the French National Research Agency (ANR) within the National Network for Software Technologies (RNTL), project Open-ViBE ("Open Platform for Virtual Brain Environments"), grant # ANR05RNTL01601, by the European COST Action B27 "Electric Neuronal Oscillations and Cognition". During the period of the research the first author has been partially supported by Nova Tech EEG, Inc., Knoxville, TN and the second author by the French Ministry of Defense (DGA).