{'fraction_non_alphanumeric': 0.032431781942898405, 'fraction_numerical': 0.009025542284665603, 'mean_word_length': 4.505134150380921, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, '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}}
Multimodal Latent Variable Analysis

Consider a set of multiple, multimodal sensors capturing a complex system or a physical phenomenon of interest. Our primary goal is to distinguish the underlying sources of variability manifested in the measured data. The first step in our analysis is to find the common source of variability present in all sensor measurements. We base our work on a recent paper, which tackles this problem with alternating diffusion (AD). In this work, we suggest to further the analysis by extracting the sensor-specific variables in addition to the common source. We propose an algorithm, which we analyze theoretically, and then demonstrate on three different applications: a synthetic example, a toy problem, and the task of fetal ECG extraction.


Introduction
The analysis of a physical phenomenon or some complex system at hand can often be made easier through the use of several sensors instead of a single complex one. The hope is that each of the sensors captures a different part of the convoluted system, while the fusion of all the information captures the global picture. This line of thinking has led to the abundance of multimodal and multi-sensory data in recent years and to an increased demand for algorithms In order to facilitate the extraction of the different sources of variability, we divide them into two conceptual categories: (i) sources of variability common to all sensors; and (ii) variables unique to a specific sensor. In our work, we focus on a two step implementation where we first reveal the common variable.
Once it is found, we extract the remaining sources of variability, i.e, the sensorspecific ones. Intuitively, our approach marginalizes the common variable, which is found in the first step, and then continues to extract the sources of variability left in the filtered data. This simplifies our task, since we do not attempt to extract all the sources manifested in the data at once.
In this paper, we use an unsupervised manifold learning approach to address the problem. Various manifold learning algorithms were proposed in the literature over the years, e.g., [2,3,4,5]. However, most of these classical methods assume that the data is captured by a single sensor, rather than in the multimodal multi-sensory setting we consider here. We focus on a particular paradigm -the Diffusion Geometry, as presented in [6,7]. Using this framework, the alternating diffusion (AD) algorithm was recently proposed in [8,9] for the purpose of extracting the source of variability common to multiple sensors. AD follows a recent line of papers that propose to use multiplications and manipulations of kernels for the purpose of fusing data from different sen-2 sors, e.g., [10,11,12,13]. Similarly to recently presented nonlinear methods, e.g., [14,15], AD is shown to reveal only the common components among all processed sensors. Successful applications of AD to real measured data were demonstrated, e.g., in [16] for the task of sleep stage identification. Herein, we rely on AD and aim to extend it by further analyzing the measurements and finding the sensor-specific variables. Our main motivation is that in some applications the sensor specific variables are far more important than the common variable. Indeed, we show one real-life example of such an application -fetal Electrocardiography (ECG) extraction.
Our main contribution in this work is a novel algorithm, attempting to recover all the sources of variability manifested in a set of multi-sensory multimodal measurements. We justify our proposed scheme theoretically, showing that it is guaranteed to find the underlying parametrizations under certain prescribed conditions. In addition, we demonstrate its applicability in three different applications: a synthetic example, a toy problem and a real-life application. This paper is organized as follows. In Section 2 we introduce formally the problem we address, and in Section 3 we review the diffusion maps and AD algorithms. In Section 4 we present the proposed method and in Section 5 we analyze it theoretically. In Section 6 we test our method on a synthetic example, a toy problem and a real-life application -the extraction of fetal ECG. We conclude this paper in Section 7.

Problem Formulation
Consider three latent random variables X, Y and Z in R dx , R dy and R dz , respectively, which are jointly distributed according to some probability density function (PDF) denoted by P (X, Y, Z). Following the work in [8], we assume that the variables Y and Z are independent given X, i.e., the joint PDF can be written as follows: 3 where P (X) is the marginal PDF of X, and P (Y |X) and P (Z|X) are the conditional PDFs of Y and Z given X, respectively. When measuring a system of interest, a measurement instance is defined by the triplet (x i , y i , z i ), which is a realization sampled from P (X, Y, Z). We do not have access to the latent variables; instead, we have two sensors observing the system at hand through two unknown observation functions given by g(x i , y i ) and h(x i , z i ). We assume g and h are smooth and locally invertible bilipschitz functions. Let {s and {s (2) i } N i=1 denote two sets of N measurement samples, taken simultaneously from the two sensors, such that s are N realizations of the system's hidden variables. In other words, we have hidden realizations (x i , y i , z i ) of three underlying variables and two sensor observations s (1) i and s (2) i ; x i is the common latent variable between the two observations, whereas y i and z i are two sensor-specific variables.
Given the two sets of measurement samples, the work in [8] showed that a method based on AD operators extracts a parameterization of the common variable X. In this work, we aim to further the analysis and extract a parametrization of the variables Y and Z as well. Such a complementing capability enables us to fully parametrize all the hidden variables underlying the measurements of the system of interest.
Although the analysis and methods used in this paper will be carried out from a different standpoint, the factorization in (1) can be used to explain the main concept. Intuitively, the extraction of the common variable X in [8] can be viewed as a marginalization operator applied to the joint probability P (X, Y, Z) obtaining P (X). In this work, we devise another operator which uses P (X) to construct the conditional probabilities P (Y |X) and P (Z|X). Then, given P (Y |X) and P (Z|X), it marginalizes the variable X and obtains a parametrization of the sensor-specific variables Y and Z.

Diffusion Maps
Diffusion maps [6,7] is a data-driven nonlinear dimensionality reduction algorithm. Given a set of N measurements {u i } N i=1 in R m , the method commences by constructing an affinity matrix W of size N ×N , whose (i, j)-th entry is given by Intuitively, W can be interpreted as a weight matrix of a graph with N vertices, where the coefficient > 0 dictates the sparsity of the edges. If is small, most edges have a negligible, close to zero weight and the graph is effectively sparse, whereas if is large, most edges are assigned with non negligible weights and the graph is dense.
The constant is usually chosen according to the data at hand, and in this work we set it using the method suggested in [8]. Therein, the constant was chosen to be = √ i j , where i a scaling constant corresponding to the i-th vertex. In particular, i is chosen to be the mean squared distance from the i-th vertex to its k nearest neighbors. The next step is to normalize the affinity matrix W, which results in the matrix K. Various normalization procedures have been suggested in the literature [17,18], each having a different interpretation when analyzed theoretically.
In this work, K is constructed by dividing each column of W by its sum, yielding a column-stochastic matrix. As a result, K can be viewed as a transition probability matrix of a Markov chain on the graph.
Once the affinity matrix is constructed and normalized, a d-dimensional is formed according to the following nonlinear map: where φ j is the j-th left eigenvector of the matrix K and φ i j is its i-th entry, λ m j is the j-th eigenvalue (when the eigenvalues are denoted in descending order) raised to the power of m, and m > 0 is a constant. Typically, d is set to be much 5 smaller than min(d 1 , d 2 ), thereby attaining dimensionality reduction. In addition to providing compact representation, this nonlinear map attempts to reveal the essence of the data in few dimensions, accurately representing their underlying intrinsic variables. In the context of diffusion maps, special attention is given to the Euclidean distance between the embedded samplesû i . Specifically, the Euclidean distance between the embedded samplesû i approximates the Euclidean distance between the corresponding columns of K m . This distance is termed the diffusion distance, since it takes into account transition probabilities on the constructed graph consisting of m Markov chain steps. We note that diffusion distance plays a large role in the algorithm presented in this paper.
For more details, as well as the motivation behind this particular dimensionality reduction method, we refers the reader to [6].

Alternating Diffusion
Given two sets of measurement samples originating from two sensors, i.e., {s (1) i } N i=1 and {s (2) i } N i=1 , the first step in the AD algorithm is constructing two pairwise affinity matrices, W (1) and W (2) based on the Gaussian kernel The constants (1) and (2) have a similar interpretation to the one presented in Section 3.1. The algorithm proceeds by normalizing W (1) and W (2) to be column-stochastic, yielding two matrices K (1) and K (2) , where the sum of each of their columns equals one. As a result, each stochastic matrix can be interpreted as a transition probability matrix of a Markov chain on a graph whose vertices are the samples (as described in Section 3.1). In other words, the (i, j)th entry in K (1) or in K (2) represents the probability of transition to the i-th vertex from the j-th vertex in the graph. Importantly, by construction (4), K (1) describes a Markov chain that jumps in high probability from the j-th vertex to the i-th vertex if the underlying values of both X and Y are similar (namely, x i is similar to x j and y i is similar to y j ). Analogously, by construction (5), K (2) describes a Markov chain that jumps in high probability from the j-th vertex to the i-th vertex if the underlying realizations of both X and Z are similar.
Given the normalized matrices K (1) and K (2) , an AD kernel is then defined This corresponds to a transition probability matrix K consisting of two consecutive, alternating steps -the first step is employed according to K (1) and second according to K (2) . Next, by raising K to the power of m, we obtain a transition matrix K m that corresponds to 2m steps, where the odd steps correspond to K (1) and the even steps correspond to K (2) . Consequently, the odd steps jump (in high probability) to a vertex where both X and Y values are similar, while the even steps jump (in high probability) to a vertex where both X and Z are similar. As a result, after many (odd and even) steps, we maintain similarity only to the X value whereas the Y and Z may vary significantly.
In order to obtain an affinity matrix in terms of the common variable X between pairs of samples (s i ) and (s j ) (recalling that each pair of samples shares the same X value according to the model assumptions presented in Section 2), the method computes the 2 distance between the corresponding columns in the matrix K m . In [8] a rigorous analysis is provided justifying this statement. Moreover, it was suggested to use a refinement step consisting of an additional diffusion maps application where the columns of K m are the new graph vertices, resulting in a low dimensional embedding as defined in (3). In the AD setting, since the underlying variable of the affinity matrix K m is X, we shall denote the resulting embedding byx instead of the general notation of u, which was used in (3).

Proposed Method
The first step towards a full parametrization of all the latent variables underlying the measurements is finding the common latent variable X, as previously suggested, using the AD algorithm. Once the common variable is extracted, we proceed to analyzing the measurements from the first and second sensors separately. Hereafter, for the sake of brevity, we will focus on the analysis of the first sensor only, while the analysis of the second sensor is analogous. For be a neighborhood of samples consisting of samples j with a similar common variable. Formally, define where η i > 0 is a small tunable threshold. In practice, instead of fixing a threshold η i for every signal, we choose all the neighborhoods N i to be of the same size q. In other words, the η i are fixed implicitly such that the size of each neighborhood N i is equal to q. The key point in defining these neighborhoods relies on the assumption that the AD algorithm is able to successfully recover the common variable X and to suppress the sensor specific variable Y .
That is, the measurements in these neighborhoods, i.e., s For each such neighborhood, we propose to compute its sample mean and its sample covariance both in the domain of the measurements {s i ) can be seen as a Gaussian representation of the local variability of the sensor-specific variable Y around every sample, and hence, in light of the discussion above, we have a local representation of Y .
In order to get a global parametrization, we compare the local neighborhoods by means of "registration" of point clouds or Gaussian distributions. Consider 8 the following affinity kernel We have denoted by C the Moore-Penrose pseudoinverse, which is employed since the rank of the covariance matrix is lower than its dimension. This follows the underlying assumption that the dimension of the measurements d 1 is larger than the dimension of the sensor-specific variable d y . We note that the omission of the low eigenvalues and their corresponding eigenvectors, as done by the pseudoinverse, results both in denoising possible ambient noise and estimation inaccuracies, and also in the attenuation of the common variable X remainders, thereby enhancing the desired variation -that of the sensor-specific variable only.
The distance in the Gaussian kernel in (10) is a modified Mahalanobis distance between the signal samples s j , which was presented in [19], with the exception that the local covariance matrices are computed based on neighborhoods in the extracted common variable domain. In [20,21,22], this distance was used in the context of manifold learning and diffusion maps to determine an intrinsic representation of (single) sensor data, invariant to interferences and measurement modalities. Such a manipulation of the Mahalanobis distance via the neighborhood choice was suggested in [23] for the task of sea mine detection in sonar images, and in [24] of reduction of stochastic dynamical systems. There, by controlling the locality within a pre-defined training set, a new metric, which is invariant to perturbations in the appearance of the target, was presented. In our work, rather than building invariances, we use a similar approach by appropriately choosing the neighborhoods in a multi-sensor setting in order to obtain a full parametrization of all the underlying sources of variability.
Once the affinity kernel W (1) is constructed, given that it captures only the variability of the (desired) sensor-specific variable Y , we apply the standard 9 Algorithm 1: The proposed algorithm.
originating from both sensors. Output: parametrizations of the sensor specific variables Y and Z.
1. Compute the parametrizationx of the common variable X using the alternation-diffusion algorithm. (b) Compute the local mean, using Equation (8), and the local Covariance, using Equation (9). 3. Compute the affinity matrix using the Mahalanobis distance between the signals s

Theoretical Analysis
In this section, we provide a theoretical analysis, showing that indeed the proposed algorithm approximates the distance between two signal samples in terms of the sensor-specific variable. As above, without loss of generality, we focus on signal samples arising from the first sensor, and therefore, our goal is to extract the variable Y . For simplicity, in this section we omit the sensor index.
A similar derivation to the one presented in this section was done in [19] and in [23]. Here, we highlight the significant differences both in terms of the analysis and in terms of the underlying assumptions. Assumption 1. If the measurement sample s j = g(x j , y j ) belongs to the neigh- This assumption relies on the ability of AD to capture the common variable X, as was proven in [8]. By definition, if a signal sample s j is in the neighborhood of s i , then the distance between their extracted values of common variable X is small (which in practice, is controlled by the tunable threshold η i ). Here we further assume that, if a signal sample s j is in the neighborhood of s i , then the distance between their respective X values (which is small by definition) is also smaller than the distance between their associated Y values by at least one order of magnitude.
Assumption 2. Locally, for every signal sample s i , the empirical covariance matrix of the sensor-specific variable Y given the extracted common variable X is isotropic, i.e., it is given by where I is the identity matrix.
While Assumption 2 may seem to be artificial and restrictive, in Section 6, we present experimental results supporting it empirically. In addition, we note that it was used in slightly different contexts in [19,22,25] and successfully applied in many applications with real measured data. The following result follows Assumption 1 and Assumption 2.
Theorem 1. For any signal sample s i = g(x j , y j ), if Assumption 1 and Assumption 2 are satisfied, then where J y i is the Jacobian of the function g with respect to the variables Y , computed at the i-th sample s i .
Proof. Using Taylor expansion, we can linearly approximate the observation function g(x j , y j ) around the point (x i , y i ), obtaining where J x i and J y i are the Jacobians at the i-th sample, s i , with respect to the variables X and Y , respectively. The last term in (14) encapsulates all the higher order derivatives that do not appear in this linear approximation. Under Assumption 1, (14) can be rewritten as for any j ∈ N i . Using (15) and by Assumption 2, the empirical covariance around the sample s i is given by concluding our proof.
Theorem 1 shows that in order to estimate empirically the Gram matrix J y i J y i T of the Jacobian J y i , we can simply compute the empirical covariance matrix of the samples in the neighborhood of x i , where the neighborhood is defined by the AD metric. This is accomplished without the knowledge of the function g itself. In the next result, we use this Gram matrix for estimating the distance between a pair of signal samples in terms of the sensor-specific variable Y .
Theorem 2. For any two signal sample s i = g(x i , y i ) and s j = g(x j , y j ), the Euclidean distance between the corresponding realizations of the sensor-specific variable is given by Proof. In the proof of Theorem 1, we considered the Taylor expansion of the observation function g from the domain of the latent variables X and Y to the range of the measured signal. Similarly, consider the Taylor expansion of its inverse function g −1 (recalling that g is assumed bilipschitz), which is given by where Q x i and Q y i are the Jacobian matrices of g −1 with respect to the variables X and Y , respectively. Isolating the Y variable yields By applying the 2 norm to both sides, we obtain The Taylor expansions in (15) and (21) correspond to g and g −1 , respectively.
Thus, due to the inverse function theorem, we have Typically, the dimension of the measurements is larger than the sum of the dimensions of the common and sensor-specific variables, i.e., d 1 > d x + d y . As a result, the number of rows in J y i is larger than the number of columns, and hence, the Gram matrix J y i J y i T is not full-rank. Consequently, it is not invertible and one needs to employ a pseudo-inverse operator instead. By substituting (23) into (22), we obtain as required.

13
Two important notes are due at this point. One is that the above analysis is based on the Taylor expansion around the sample s i . If we repeat the derivations with the Taylor expansion around the sample s j as well, then, the mean of the two resulting expressions is given by Once such symmetrization is employed, further analysis presented in [19] improves the order of the error term in (25) to s j − s i 4 2 . Two is that the above analysis assumes (unrealistically) that every signal has a local mean equal to zero, i.e., µ i = 0. However, as presented in [22], a similar derivation can be done without such an assumption. Combining these two notes results in the expression presented in (10), without a rigorous proof.
To conclude, the analysis presented in this section coincides with the proposed method. Indeed, in Algorithm 1, we begin by seeking for signals close in terms of the common variable using the parametrization obtained from the AD algorithm. Once such a neighborhood is found, we compute its local empirical covariance matrix (9), which is then used in the modified Mahalanobis distance (10)  Final remark concerns the accuracy of the Euclidean distances approximation. Theorem 2 implies that when the distances are large, the error terms are large and the (local) approximation via the linear terms is poor. This problem is "automatically" alleviated by the standard usage of the Gaussian kernel in (10); due to its fast decay, large distances are implicitly attenuated. 14 6. Experimental Results

Synthetic Example
Consider three independent and identically distributed random variables, X, Y , Z, sampled uniformly in [0, 1]. We generate from these variables 3000 triplets of (x i , y i , z i ). Assume two sensors observing these hidden samples through the following nonlinear functions g and h and so that we obtain 3000 pairs of signal measurements (s i ). We set R = 10, r (1) = 4, r (2) = 2. Notice that g and h correspond to two tori, with major angle X, serving as their common hidden variable, and with minor angles, Y or Z, serving as their respective sensor-specific variables. We apply Algorithm 1 to these samples where the size of the neighborhoods in the common variable domain is q = 11. In Figure 1 we color both tori according to the extracted parametrization of the common variable and also according to the obtained parametrizations of the two sensor-specific variables. Indeed, we observe that our method accurately extracts the three hidden variables. The coloring of the tori according to the common variable X is highly coherent with the major angle, while the coloring with respect to the sensor-specific variables, Y and Z, are consistent with the minor angles.

Playing with Toys: Yoda, Bulldog and Rabbit
In this experiment, we consider the toy problem presented in [8]. The setting of the problem includes three objects: a figure of Yoda (green alien), a Bulldog,  In [8], it was shown that the AD algorithm attains a parametrization of the angle of the Bulldog X, namely, the common variable hidden in the sets of images. In this work we infer a parametrization of the angle of the sensor-specific Yoda Y . In Figure 3 we present the result of applying Algorithm 1 in this setup where the size of the neighborhoods in the common variable domain is q = 15. We scatter plot the first two coordinates in the obtained parametrization, y, and observe that the parametrization takes the shape of a circle, correctly representing a rotating angle. To show that this angle is indeed associated with the rotating angle of Yoda, we overlay several images corresponding to the embedded points 1 .
It can be seen that the orientation angle of Yoda corresponds to angles on the obtained circle.
For further evaluation, we compute the ground truth parametrization of the angle of Yoda, denoted byỹ. To this end, we crop all images captured by the first camera, discarding Bulldog and maintaining only Yoda, and then we apply diffusion maps to the set of cropped images. We emphasize that the information of how to appropriately crop the image is not available to the proposed algorithm, which does not use any prior knowledge on the experimental setting, and it was done for illustration and evaluation purposes only. In Figure   4, we present the discrete Fourier transform of both the first components ofỹ and y; the parametrizationỹ is obtained by diffusion maps applied to the cropped images, and the parametrizationŷ is obtained by Algorithm 1. We observe a sharp peak around the frequency 310, which according to [8] corresponds exactly to the rotation speed of Yoda 2 . Moreover, the curves are similar, implying on the successful recovery of the sensor-specific variable by the proposed algorithm.

Non-Invasive Fetal ECG
Fetal heart rate monitoring [26,27,28] is widely-used for the assessment of contains, in addition to the fetal's heart beats, the maternal ECG, masking the desired information. In order to suppress the maternal ECG and to extract the fetal ECG, another (reference) electrode is often placed on the mother's thorax for the purpose of measuring only the maternal ECG.
In practice, in addition to being occluded by the maternal ECG, the fetal ECG is also contaminated by noise. Power line disturbance and maternal muscle movements (electromyographic activity in the abdomen and uterus muscles of the mother) are only two typical examples of the possible interferences hindering the extraction of the fetal ECG [29].
As reported in [29], due to its time-varying statistical character, the ECG of the fetal is a highly nonstationary signal. Moreover, the relation between the measured abdomen signal and the fetal ECG is arguably nonlinear. As such, standard approaches, e.g., the adaptive least mean squares (LMS) algorithm [30], provide only coarse estimations in recovering the fetal ECG, and the solution for this problem is not trivial, and it is still considered an open 19 problem.
In [29], the authors suggest to tackle the fetal ECG extraction problem by first extracting the maternal ECG from the two measurements using adaptive noise cancelers. Then, given the result, that is the fetal ECG plus muscle noise, the authors suggest to employ an adaptive signal enhancer in order to extract the fetal ECG and attenuate the remaining noise. In particular, this adaptive signal enhancer relies on the alignment in time of an ensemble of similar pulses and the extraction of their statistics. The solution presented in that work is specificallytailored for fetal ECG extraction as it requires, for example, the detection of the peaks in the fetal ECG signal using some peak detector. In contrast, in the sequel, we show that our approach does not require any knowledge about the task at hand.
We use the fetal ECG extraction problem as a testbed for our algorithm not only to demonstrate its applicability, but also to show the relevance of the problem setting we present in this paper to real measured data. Let us now return to the problem formulation, as defined in Section 2. In our context, the variable common to both the abdomen and thorax signals is the maternal ECG, while the sensor-specific variable in the abdomen signal is the desired fetal ECG.
We demonstrate that our proposed method is capable of not only recovering the maternal ECG (common variable), but also factoring out the mother's pulse from the measured abdomen signal. This results in revealing the fetal ECG (sensor-specific variable), which is relatively weak when compared to the maternal signal. More specifically, we show that our method builds a parameterization of the fetal ECG, which, in turn, could aid in detecting fetal QRS complexes which are used in measuring the fetal's heart rate.
In our experiments we use the "Non-invasive Fetal ECG Database" from PhysioNet [31], and apply our algorithm to the raw data without any preprocessing. In Figure 5 we present a short 2 second interval from both the thorax and the abdomen signals. The length of the entire signals is 5 minutes and 20 seconds and in the following experiments we operate on a sample extracted from them that is of length 32.76 seconds. The sampling rate of the ECG signals is Given these two signals, we first apply the AD algorithm to extract the common variable and present the obtained parametrization in Figure 6. The algorithm is applied to time segments of length 256 samples (lag map) with 16 samples overlap, which are extracted from both the abdomen and thorax signals 3 . Since the sampling rate is 1 kHz, each segment is of duration 256 miliseconds. In the context of this paper, these segments are viewed as the sensor samples, and thus, are denoted by s (1) i and s (2) i . Each 2D point in the scatter plot in Figure 6, representing a pair of segments, is colored according to the angle created between the axis origin and the point itself (we choose this method of coloring due to the parametrization resembling a circle stemming from the signals periodicity). To emphasize the validity of our assumption that the common variable is indeed related to the maternal ECG, we present in Figure 7 the thorax signal and color its samples according to the extracted parametrization of the common variable. Clearly, the common variable coincides with the cardiac cycle of the mother. 3 The mean of every segment was subtracted.  Next, we proceed by extracting the sensor-specific variable from the abdomen signal using our proposed algorithm. In this experiment we set the size of the neighborhoods in the common variable domain to be q = 21. In Figure 8, we present the abdomen signal colored according to the obtained parametrization of the sensor-specific variable. The results imply that indeed the sensor-specific variable is related to the ECG of the fetal. Importantly, at t = 70.5 sec, t = 77.8 sec, and t = 88.4, we observe that the fetal's heart beat is detected, even in pathological cases where it is completely "buried" in the maternal heart beat.
To demonstrate the generality of our method, we repeated the experiment and applied our algorithm to a signal measured from another patient. The results of this experiment are presented in Figure 9, showing that the parametrization of the sensor-specific variable manages to capture the ECG of the fetus in this case as well. Similarly, at times t = 69.9 sec, t = 78.5 sec, and t = 92.7 sec the algorithm manages to capture the fetal ECG despite the significant occlusion by the maternal heart beat. This result also demonstrates cases in which the identification fails. For example, at times t = 78.2 sec and t = 93.4 sec a (possibly) redundant fetal heart peak is detected. One should note however that due to the lack of a ground truth, we can not be certain whether this is indeed a mis-identification. It might be the case that these are anomalies in the fetal heart rate and that the prediction is in fact correct.
To better support our results, we further apply our algorithm to two additional patients and depict the results in Figure 10. Here as well, it is observed that our algorithm manages to detect all the spikes corresponding to the QRS complexes of the fetus, exhibiting robustness to abnormal and uncharacteristic peaks, such as the one at time t = 43.8 sec.
We note that in order to establish baseline results, in addition to using the proposed algorithm, we have attempted to use simpler methods such ICA [32,33] in order to separate the maternal ECG from the fetal ECG. Broadly, these methods attempt to transform the data at hand (after some whitening) into components that are as statistically independent from each other as possible.
However, these methods did not obtain satisfactory results in separating the fetal and maternal ECG signals, since some nontrivial preprocessing (in addition to whitening and dimensionality reduction) must be employed in order to facilitate their employment. Conversely, we emphasize that our proposed method does not rely on any preprocessing of the data (except for the segment mean subtraction).

Conclusions
Given a set of measurements, originating from several sensors, the AD algorithm extracts a parametrization of a variable common to all sources. In 23 this work, leveraging on AD, we proposed a method which further analyzes the signals by extracting the sensor-specific variables. We provided a theoretical justification as well as various applications. A shortcoming of our method is the need to extract the intermediate common variable parametrization. Proposing a method that could skip this stage is a promising future direction.  . The signal is colored similarly as in Figure 8 and Figure 9. Notice the uncharacteristic peak at time t = 43.8 sec that is captured by our algorithm.