{'fraction_non_alphanumeric': 0.03479083455959638, 'fraction_numerical': 0.0059386167752785365, 'mean_word_length': 4.937285491419657, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Hypothesis Testing Approach to Nonstationary Source Separation

The extraction of nonstationary signals from blind and semi-blind multivariate observations is a recurrent problem. Numerous algorithms have been developed for this problem, which are based on the exact or approximate joint diagonalization of second or higher order cumulant matrices/tensors of multichannel data. While a great body of research has been dedicated to joint diagonalization algorithms, the selection of the diagonalized matrix/tensor set remains highly problem-specific. Herein, various methods for nonstationarity identification are reviewed and a new general framework based on hypothesis testing is proposed, which results in a classification/clustering perspective to semi-blind source separation of nonstationary components. The proposed method is applied to noninvasive fetal ECG extraction, as case study.


INTRODUCTION
A broad class of blind and semi-blind source separation (BSS) algorithms are based on the joint (approximate) diagonalization of a set of matrices/tensors containing second or higher order statistics of multichannel data. While a great body of research has focused on efficient methods and algorithms for approximate joint diagonalization (AJD) of a set of matrices [1,Chs. 5 and 7], the original selection of these matrices remains a highly problem-specific and nontrivial issue. The problem is more critical for mixtures of nonstationary or sparse signals, for which the desired signals of interest are difficult to detect and are considered as outliers that do not necessarily influence the average statistics captured by the diagonalized matrices/tensors. In fact, many BSS algorithms perform source separation based on the average statistics from the entire dataset and are incapable to identify local variations in the signal statistics that do not have significant effects on the global statistics.
Herein, we review some of the existing and new methods for supervised and unsupervised nonstationary sample identification (and clustering) during the formation of the matrices used in (approximate) joint diagonalization. As a result, the diagonalized matrices are calculated more targeted (and prior-based), making them more appropriate for nonstationary source signals. The idea is shown to be of great interest for unsupervised and semi-supervised source separation, and closely related to the problem of clustering and classification of samples in imbalanced and nonstationary scenarios.
In Section 2, we review methods for nonstationarity detection. These methods are later used in Section  how nonstationarity can be used as a means of source separation. As proof of concept, in Section 4, the proposed method is applied to the problem of noninvasive fetal electrocardiogram extraction.

NONSTATIONARITY IDENTIFICATION
Nonstationarity generally refers to signal attribute variations in time, space, frequency, time-frequency, stochastic properties, etc. In the sequel, we review the concept of nonstationary in different domains and how it can be quantitatively identified. Our data model is where k denotes time, x k ∈ R n are multivariate observations, A k is a linear invertible matrix and s k ∈ R n is the source vector.

Density Function Variations
The probability density function (PDF) of the observations in (1) is which is a function of the source vector density function and the (generally time-variant) mixing matrix. Apparently, the temporal nonstationarity of the source signals or the mixing matrix result in local variations of fX (·; k). Now suppose that regardless of the local variations due to nonstationarity one fits an empirical PDFfX (x) over the observed samples. Noting thatfX (x), or equivalently the empirical cumulative distribution function (CDF)FX (x), are no longer functions of time (k), one can conduct a sample-wise hypothesis test for nonstationarity detection, to identify the "odd events" that locally differ from the global density functionfX (x): Various statistical tests exist for evaluating whether or not a given data sample is drawn from a given PDF. Considering that the observations are commonly band-limited signals that satisfy the Nyquist rate, a nonstationary event is not an isolated point, but rather consists of a series of successive samples. Now suppose that we want to assess whether a segment of p successive data points {x k−p+1 , . . . , x k } corresponds to a nonstationary event, i.e., the null hypothesis in (3) holds for these samples. In an Anderson-Darling test the following statistics is calculated: where {z1, . . . , zp} is the sorted version of {x k−p+1 , . . . , x k }, in non-decreasing order. Depending on the segment length p, A 2 k can be checked versus tabulated values to accept or reject the null hypothesis with a desired confidence (p-value) [2,Ch. 4]. As a continuous index, A 2 k can be updated over time and compared with desired thresholds to identify the nonstationary points.

Reference-Based Nonstationary Identification
A popular approach for nonstationarity identification is to use a reference channel to solve the hypothesis test in (3). Various forms of this approach have been used in the literature, mostly without explicit reference to a hypothesis test.
In [4], the R-peaks of electrocardiogram (ECG) signals were detected from a reference channel and used as a means of deriving variable periods of the ECG (called the cardiac phase signal), which were used for extracting the ECG from noise by periodic component analysis. In [5], a reference lead was used to identify the artifacts due to a magnetic resonance device that obscured simultaneously acquired EEG signals. In [6], electrooculogram (EOG) artifacts of the electroencephalogram (EEG) were identified by thresholding the energy envelope of an EOG reference lead, which was closely placed to the eye. In [7], thoracic electrodes were used to identify the R-peaks of maternal ECG from the abdominal leads used for fetal ECG extraction. In [8,9], a reference-based method was used for defining periods including interictal epileptiform discharges for identification of brain regions involved in a reference state from intracerebral EEG.

Observations Moments/Cumulants Tracking
The moment generating function of x k defined as [10, Ch. 2]: uniquely characterizes the distribution of the random vector x k . Using the Taylor expansion of exp(·), variations of the PDF can be associated to and tracked by monitoring the first (mean), second (covariance), and higher-order moments and cumulants of the observations. This gives an alternative solution for solving the hypothesis test proposed in Section 2.1. Accordingly, depending on the source of nonstationarity, one can estimate and track the first few cumulants of the observations over sliding windows, to identify the nonstationary events.

Adaptive Source Separation
Adaptive source separation algorithms are well-known in the literature. As proof of concept, let us consider the equivariant adaptive source separation via independence (EASI) [11]. In this method, the separating matrix, here denoted by W k , is adaptively updated using an equivariant serial update: where y k = W k x k are adaptive estimates of the independent sources and H(·) is a nonlinear function of the estimated sources cumulants. For time-varying mixtures, the algorithm seeks the demixing matrix such that in steady state W k A k approaches identity and H(y k ) → 0. Now suppose that the mixing matrix A k = A is constant, but the sources are nonstationary. As a result, the function H(·) deviates from zero as the recursion reaches the nonstationary epochs.

Spectral Nonstationarity
Spectral variations are another common type of nonstationarity [12], which can be identified by using time-frequency signatures of the nonstationary events versus the background activity. This is a wellknown technique in the literature (see [13,Ch. 4.3] for a hypothesis test approach to time-frequency detector design). An alternative approach is to track the variations of the auto-correlation function, estimated over sliding windows of the observations. If one adopts a model-based approach to auto-correlation estimation-e.g. by fitting a time-varying auto-regressive (AR) model over sliding windows of the data-temporal tracking of the model coefficients (or equivalently the root loci of the fitted model) is an efficient means for tracking the temporal nonstationarities [14].

Temporally Dynamic Sources and Innovations Monitoring
Nonstationarity identification has also be performed in the time domain, by state-space modeling ( [15,7]). Suppose that the underlying multivariate signal of interest s k is modeled by a dynamic system: where f k (·, ·) and g k (·, ·) are generally nonlinear, differentiable and smooth functions considered as process and observation models, respectively. A very popular class of estimation schemes for such nonlinear models (including the Kalman or extended Kalman filter family) is in the form of a recursive estimator: whereŝ − k andŝ + k denote the state vector estimates before and after using the kth observation, respectively; i k is known as the innovation process; and K k is a gain matrix to be determined based on the sample-wise estimation qualities (covariance matrices). For a well-functioning filter, the innovation process is spectrally white (its auto-correlation is an impulse). In fact, the whiteness of the innovation process spectra is an indicator of "optimal filtering", showing that all the predictable parts of the observations have been retrieved, leaving only white noise behind. Now suppose that during the recursive time updates, the innovation process starts to deviate from white noise (becomes correlated in time). This behavior can be associated to changes in the underlying system's dynamics due to nonstationarity in the mixed signals or model parameters. This property was recently used in [7] to set up a hypothesis test for the noninvasive detection and extraction of fetal ECG.

Insufficient Information
Nonstationary events can suffer from the problem of insufficient information. In fact, for short nonstationary events, the sample estimates of the second or higher-order matrices/tensors are inaccurate. Therefore, model-based approaches are required to improve the statistics estimated from limited nonstationary points. Recent developments in random matrix theory can also be used in this context, where methods have been developed that guarantee the accuracy of eigenvalue decomposition, despite the inaccuracies of sample covariance matrices [16,17]. The extension of these techniques to AJD is currently an open problem.

Rate of Change and Nonstationarity Transitions
Another point of reservation is the rate of signal attribute variations and the margins of nonstationary time windows. This problem has been addressed in the theory of abrupt change detection [18], which is highly influential for the hereby studied problem.

NONSTATIONARY COMPONENT ANALYSIS (NSCA)
To this point, we reviewed various methods for detecting nonstationary events. In this section, we study how the detected nonstationarities can be used during source separation, assuming a constant mixing matrix A k = A.

Single Hypothesis Tests and Generalized Eigenvalue Decomposition
Consider multivariate signals x k ∈ R n (k ∈ T ), where T denotes the entire set of available samples and P ⊂ T is a subset of the sample points considered as nonstationary or odd events, with statistical properties that differ from the other samples. For this problem, the sample-wise hypothesis test in (3) can be restated as follows: H0 : k / ∈ P, H1 : k ∈ P The covariance matrices of the data over the samples of each hypothesis and the total data covariance matrix are defined as follows: where E{·} denotes averaging over the indicated subset of samples, and m0 the class-wise and total sample means, respectively. Now, defining y k = W T x k (11) where W ∈ R n×n , generalized eigenvalue decomposition (GEVD) of the matrix pair (C0, Cx) can be used to obtain W, such that: where Λ = diag(λ1, . . . , λn) is a diagonal matrix of generalized eigenvalues corresponding to the eigenmatrix W = [w1, . . . , wn], with real eigenvalues sorted in ascending order on its diagonal. In addition, the first eigenvector w1, corresponding to the largest generalized eigenvalue λ1, maximizes the Rayleigh quotient [19]: By using the Rayleigh quotient, it is straightforward to show that for single hypothesis test in (9), the joint diagonalizer is exactly equivalent to the one obtained by GEVD of the matrix pair (C1, Cx), but this time sorting the eigenvalues in descending order. Another property of this transform is that W retrieves uncorrelated linear mixtures of x k and y k = w T 1 x k has maximal energy during the nonstationary time epochs. Therefore, W extracts the nonstationary events, corresponding to H1, from the other background signals.

Multiple Hypothesis Tests and Approximate Joint Diagonalization of Eigenmatrices
The idea of sample clustering before source separation, can be extended to multiple classes, using a multiple hypothesis test. This idea is very close to a well-known method for separation of nonstationary sources, developed in [20]. In this case, the total data sample set T can be partitioned into K classes/clusters: T = P1 ∪ . . . ∪ PK , resulting in a multiple hypothesis test for the time samples: and the covariance matrices of each class is defined: where mi We may now seek the joint (approximate) diagonalizer W ∈ R n×n , such that are "as diagonal as possible." It is known that for K > 2, the diagonalization is only achieved approximately by using different variants of AJD. Depending on the application and diagonalization algorithm, the total covariance matrix Cx, defined in (10) may also be among the set of matrices to be diagonalized, to achieve decorrelated sources 1 .

Source Separation vs Clustering/Classification
The problem of (approximate) joint diagonalization of a set of matrices is closely related to the problem of clustering and classification, which has readily been used in common spatial patterns and similar algorithms [22,23,24]. In the two-class (two-cluster) case, the relationship is explicitly evident from the maximization of the Rayleigh quotient in (13), which can be seen as the linear separator between two classes represented by their covariance matrices. The relationship in the multiple hypothesis test case (detailed in Section 3.2), can be non-rigorously explained as follows: by excluding an arbitrary subset of the sample points Pj, one can set up a two class hypothesis test: H0 : k / ∈ Pj, H1 : k ∈ Pj (17) which is a test of one versus the fusion of all other classes. The covariance matrix corresponding to the null hypothesis is: where the differences of the cardinality of Pj and the significance of each matrix can be adjusted by the weights wi. Now, suppose that the separation matrix W satisfies (16), i.e. using some typical measure of joint diagonality, Λi are "as diagonal as possible" [1]. Therefore, the matrix C T \P j is also "approximately diagonalized" by W, and we have Therefore, defining the vector of eigenvalue ratios for each j, the eigenvector corresponding to the maximum entry of this vector, is "approximately" the separator of class samples j versus the union of the other sample labels. In other words, if m is the index of the maximum entry of the vector dj , then y k = w T m x k mostly resembles the samples of class Pj. Since this property holds for all j, the separation method based on AJD can be seen as a method that maximizes the total between-class distances. However, the relationship between m and j is not necessarily one-to-one. In other words, for any class j, represented by its covariance matrix Cj, there exists a separating vector wm that maximally discriminates class j from the union of others. But the same vector may happen to be the maximal separator for another class as well. Note also that the result of AJD is different from the solution obtained from the exact one versus the union of all other classes, which would be obtained by the GEVD of the matrix pair (Cj, C T \P j ). The latter guarantees the maximal distance between class j and the other classes; but does not attempt to maximize the inter-distance between the other classes.
The above property can be used for a "two-round targeted source separation", which first applies unsupervised source-separation algorithms such as JADE [25], to extract independent components (n clusters). In the second round, any desired channel among the estimated independent components can be refined by performing a twoclass NSCA between the desired channel versus the other components. From the clustering perspective, the first unsupervised ICA can be considered as n-group clustering, and the second round being a one vs (n − 1) classifier.

CASE STUDY
As a case study, we consider the problem of noninvasive fetal ECG (fECG) extraction from maternal abdominal recordings. In Fig. 1(a), a sample fECG signal from 5 maternal abdominal and 3 maternal chest leads is shown [26]. In Fig. 1(c), the following four indexes introduced in Section 2 are shown for proof of concept: innovations index after maternal ECG (mECG) cancellation (cf. [7] for details), Anderson-Darling test with a Gaussian kernel fit, energy envelope, and EASI nonlinear function norm. The signals and indexes have been normalized for visual comparison. For this example, the innovation process index, has been the most discriminative, which is due to the fact that the innovation process after mECG cancellation mainly consists of the fECG and background noise. This index was next thresholded at a level of 50% of its peak value to identify the nonstationary time epochs corresponding to the fECG. Fig. 2 demonstrates the results of NSCA by using this index, which shows the fECG in its first channels. The Matlab codes for generating these results are online available in the open-source electrophysiological toolbox (OSET) [27].

CONCLUSION
In this paper, various criteria which have been developed for the extraction of nonstationary components from multivariate data were reviewed from a unified hypothesis test perspective. It was shown how the problem is closely related to the joint (approximate) diagonalization of matrices/tensors containing second and higher-order statistics of the nonstationarty events and background stationary activities. The proposed approach was also shown to have a clustering/classification interpretation. Various aspects of this study can be  extended in future research, including: 1) integration with deflation algorithms such as [28], for low-rank and degenerate mixtures; 2) using developments in random matrix theory to overcome the problem of low number of samples during nonstationary epochs of the data [16,17]; 3) explicit integration of NSCA with classification and clustering algorithms.