Noninvasive fetal ECG extraction using doubly constrained block-term decomposition

Fetal electrocardiogram (fECG) monitoring is a beneficial method for assessing fetal health and diagnosing the fetal cardiac condition during pregnancy. In this study, an algorithm is proposed to extract fECG from maternal abdominal signals based on doubly constrained block-term (DoCoBT) tensor decomposition. This tensor decomposition method is constrained by quasiperiodicity constraints of fetal and maternal ECG signals. Tensor decompositions are more powerful tools than matrix decomposition, due to employing more information for source separation. Tensorizing abdominal signals and using periodicity constraints of fetal and maternal ECG, appropriately separates subspaces of the mother, the fetus(es) and noise. The quantitative and qualitative results of the proposed method show improved performance of DoCoBT decomposition versus other tensor and matrix decomposition methods in noisy conditions.


Introduction
The abdominal fetal electrocardiogram (fECG) is the electrical activity of the fetal heart that is recorded by abdominal electrodes. It is a very important tool for monitoring and evaluation of fetal cardiac activity. This signal is generated from a small heart, so the signal amplitude is low and it is morphologically quite similar to the adult ECG, with a higher heart rate. Non-invasive fECG extraction means extracting a clean fECG signal from the abdominal mixed signals, which is an interesting issue for several reasons [1]. The fECG provides useful information about fetal health for clinicians and helps them in better diagnosis. It also provides more precise data on fetal heart rate compared to ultrasonic Doppler techniques [2]. Nevertheless, there are limitations in non-invasive fECG extraction, due to the presence of interference signals such as the much stronger maternal ECG as the main interference source, muscle noise, motion artifacts, and etc. Moreover, since fECG and maternal ECG have temporal and spectral overlap, fECG separation from the abdominal mixture is complicated [3]. These limitations could be considered as the major difficulties in fECG extraction, which have turned it to one of the major challenges in biomedical signal processing.
In the field of non-invasive fECG extraction, there are many mathematical based methods such as Kalman Filter [4][5][6], wavelet transform [7][8][9][10][11], blind source separation [12][13][14][15][16][17], and matrix decomposition [18,19]. These methods attempt to detect R-peaks or extract the fECG from either abdominal or abdominal and thoracic recordings. In recent years, an approach has been developed in which, tensor decompositions have been employed for signal processing and fECG extraction, whereby promising results have been achieved [20]. An effective solution for the extraction and separation of sources from tensorized observations could be obtained by constrained tensor decomposition and modified cost function. In addition, a correct selection of decomposition type and the method of constructing tensors from observation data are very important in achieving significant results.
In [21], a tensor decomposition-based algorithm for fECG extraction is introduced by adding fetal periodicity constraint to Tucker decomposition. Because of the low-amplitude of the fetal signal compared to the high-amplitude of the maternal signal, the drawback of this method, especially in Low SNRs, is the use of only the fetal periodicity constraint.
Block-term decomposition (BTD) decomposes a tensor into tensors with a rank greater than one. In the application of fECG extraction, by changing the cost function and applying constraints on it, the structure of the BTD can be modified to decompose the tensor into three blocks containing maternal, fetal, and noise subspaces. In this paper, by using some source separation methods (πCA or ICA), the fECG peaks are first extracted and then, by using the extracted subspace from the observation matrix (abdominal mixed signals), the observation tensor is formed, and finally by applying periodicity constraints on BTD, the fECG signal could be extracted. In fact, by utilizing the useful information of the ECG signal, the semi-periodicity property, is the objective of constraining the BT decomposition. Regarding the semi-periodicity of both the maternal ECG (mECG) and fECG signals, we apply two constraints on the BTD, simultaneously.
The remainder of this paper is organized as follows. Section 2 presents the concepts and relations of BTD. Section 3 involves how to construct the observations tensor and the proposed constrained BTD. Finally, the results and conclusion are presented in sections 4 and 5, respectively.

Block term decomposition
A tensor is a multidimensional array, which is the higher-order generalization of vector and matrix [22]. Decompositions of tensors have applications in many fields such as signal processing and image processing. CP (Canonical Polyadic) decomposition approximates a tensor by a sum of Rank-1 tensors. In many applications, e.g. blind source separation (BSS), a given tensor should be decomposed into tensors of rank higher than one. BTD has been introduced for this purpose. In contrast to the CPD, BTD approximates a tensor by a sum of higher-rank terms ( Figure 1). In this section, definitions and characteristics of BTD are represented. A block term decomposition of a tensor is a decomposition of of the form where the main tensor is decomposed to a sum of Rank-( , , ) terms; A , , and C are the factor matrices where , , , and are the core tensors. A visual representation of a third-order tensor is shown in Figure 2. Another type of BT decomposition writes a third-order tensor as a sum of R Rank-( , ,1) terms. Each of these terms can be written as the outer product of a matrix and a vector ( The matricized form of BT decomposition based on factor matrices could be expressed as follows , (E ) (E ) (E )-C where A, and C are more general matrices and considered as A ,A A A -, , -, and C , -. Figure 3 shows this decomposition. terms [24].
The optimal estimate of the factor matrices A , and C could be obtained by solving the optimization problem (4) by using the alternating least square method [23].
Other tensor decomposition methods and more details of BTD are presented in [22,23,25,26].

Proposed method: DoCoBT decomposition
In this section, the proposed method for fECG extraction based on tensor decomposition is explained in two parts: The doubly periodicity constraints and the DOCoBT decomposition. Figure 4 shows the block diagram of the proposed method and contains tensor construction by using maternal abdominal signals, constrained BT decomposition by using maternal and fetal periodicity constraints, detection, and selection of fetal sources and reconstruction of fECG in the sensor space.

Tensor construction
The tensor used in this study is a third-order one with channel time subspace structure (L = number of channels, N = number of time indecis, K = number of subspaces).
In order to create this tensor, the abdominal signals and fECG subspace obtained by using a source separation methods (ICA or CA) are placed in the first and second slice of the third mode of the tensor, respectively. Figure 5 shows the structure of these tensors.

Fetal R-peak detection
To eliminate user interference and to speed up processing for fECG extraction, it is important to automatically identify the type of sources (maternal, noise and fetal sources) obtained from source separation methods. The type of sources can be determined by employing the maternal and fetal R-peaks and using correlation analysis. In order to extract fetal R-peaks, a πCA-based algorithm is proposed in this section ( Figure 6). Due to the dominance of mECG signals, at first, the maternal Rpeaks are extracted using one of the common peak detection methods [27,28]. Then, the observation m trix is de omposed into its sour es using πCA nd m tern l R-pe ks In the πCA lgorithm, the extracted sources are ranked according to their synchronization degree (periodicity) with a vector of peaks impulse train. Therefore, since the maternal R-pe ks re used s the impulse tr in in the πCA method, the extracted sources are ordered based on their similarity to the mECG. So, the most similar sources to the mECG could be observed in the source matrix from the first to the last extracted component. Consequently, the last signals have the least similarity to the mECG, and accordingly, the fetal components could be observed in the last rows. Since the average beat-rate for mother and fetus are different significantly, after extracting the R-peaks of the last signals, we calculate average RR interval and compare it with a threshold T = 0.5 sec (the average RR interval for mother and fetus are about 0.9 sec and 0.45 sec respectively). The first semi-periodic source with average RR interval less than 0.5 sec is recognized as fetal signal and its R-peaks are utilized as fetal R-peaks.

Periodicity constraint
Biomedical signals usually contain information that could be employed for separation and extraction of desired sources. The semi-periodicity of ECG and fECG signals could be considered as the mentioned information.
In [29], the cost function (5) has been introduced as a criterion of periodicity for semi-periodic signals. Thus, the maximization of the following cost function could be used as a means of extracting the semi-periodic signal ( ) ( 5 ) ε( ) * ( ) ( )+ * ( ) + where *+ represents averaging over time index , and is the time-varying period of ( ). This maximization problem could be reformed as a minimization problem: In [21], the penalty (5) has been used in order to force Tucker decomposition (called it Tucker) to extract fetal components as periodic signals. In this method, the mother's ECG as a strong periodic interference has been ignored.

DoCoBT decomposition
In order to extract the fECG from abdominal signals, we first construct an observation tensor and then decompose it into three blocks or tensors using BTD. The first and second blocks are constrained by fetal and maternal ECG periodicity constraints, respectively, and the third one is considered as a noise block (Figure 7).
Consider the tensor has the structure of channel time subspace. Therefore, , , -( ) and C , -( ) contain channel, temporal and subspace components, respectively. Thus, the penalty (6) could be written for B 1 and B 2 factor matrices that are initialized with R 1 fetal and R 2 maternal sources, extracted from a BSS algorithm such as ICA where and are the fECG and mECG time-varying periods, respectively. ( , ) and ( , ), as th rows of the factor matrices B 1 and B 2 , contain the th time index of fetal and maternal components, respectively. Moreover, B 3 is initialized with R 3 noise components extracted from the BSS algorithm. Finally, to estimate matrix B (sources matrix), the following optimization problem with doubly periodicity constraints could be obtained Since the optimization problem (13) does not have a closed-form solution, one can use the gradient descent algorithm to update the factor matrices B 1 , B 2 and B 3 , alternatively. So the derivatives , , and are calculated as follows: Using Eq 15 we introduce Algorithm 1 to estimate matrices A, , C and to extract the fetal ECG using doubly constrained block term decomposition. We named this algorithm as -doubly constrained block term (DoCo )‖

Results
In this section, we present the datasets and evaluation measures at first. Then, the proposed method using both synthetic and real data is evaluated.

Synthetic dataset
In order to evaluate the proposed method, we use a synthetic dataset that resembles real fECG signals. These data must contain mECG as the main interference and other noises such electrode movement and clear fECG to compute signal to noise ratio. The open-source electrophysiological toolbox (OSET) contains the required functions for generating synthetic maternal-fetal ECG datasets, using the following data model for modeling maternal abdominal signals [30].
is the maternal ECG as interference and is the fetal ECG. and represent colored and white noises respectively. Therefore, parameter controls the signal-to-interference ratio ( ) and the parameters and control the signal-to-noise ratio ( ). In order to evaluate the quality of the extracted fECG, the signal-to-interference and noise ratio ( ) is calculated for the extracted fECG (output) signal and the is calculated as a measure of fECG denoising performance [31]: -. The synthetic data that are generated and used in this paper have the dimensionality of 8 10000 (8 channel and 10000 time indices and thus the size of the tensor is 8 10000 2).

Real dataset
In addition to the synthetic signals, the real dataset Abdominal and Direct Fetal Electrocardiogram Database (ADFED) is used to evaluate the generalization of the proposed method. ADFED consists of five channels (four channels of abdominal ECG and one channel of fetal scalp ECG) with 1 kHz sampling rate and length of 10 seconds [32].

Evaluation of the proposed method
As for real data, the clean fECG is unknown. Therefore, the SINR and the performance of the proposed method can not be calculated quantitatively. Therefore, the proposed approach is evaluated on simulated data with different noise levels. To attain this goal, 10,000 samples of eight lead synthetic ECG mixtures with SINRinput = {-5,-10,-15,-20,-25} and 1000 Hz sampling rate are generated 10 times per SINR. The performance of the proposed method is compared with some related source separation methods such as ICA (using the JADE algorithm) [33], πCA [29] and π u ker [21] In πCA nd ICA methods, first, the components of the mother and fetus are extracted. Then, the fetal subspace is reconstructed using the fetal components. Regarding the π u ker method, the fetal sources are extra ted from the o serv tion tensor y using πCA or ICA hen, the mode 2 factor matrix including the fetal sources is updated, iteratively. Finally, the main tensor is reconstructed again via the updated mode 2 factor matrix.
In the proposed method, we need to have an estimation of the fetal and maternal signals as initial sources. Therefore, the observation tensor has to be matricized to extract initial sources using source separation methods. For automatic detection of mECG and fECG sources, the correlation criterion is applied between the R-peaks of these initial sources and the R-peaks of maternal/fetal signals derived from the maternal/fetal R-peak detection algorithm (section 3.2). Finally, the determined fetal and maternal sources are placed in the B 1 and B 2 factor matrices as the initial estimation of the tensor decomposition. The sources that do not pass the correlation criterion are considered as noise and are placed in the factor matrix B 3 .
Here, BTD decomposes the observations tensor into three tensors. According to the proposed method, the factor matrices B 1 , B 2 and B 3 are updated using the gradient descent algorithm, alternatively. However, the matrix B 3 could be calculated by solving , considering that no constraint is imposed on it. Eventually, the first, second and third blocks are reconstructed via updated factor matrices. The first slices of the first block and the second block represent fECG and mECG subspaces, respectively. Moreover, to obtain parameters , and , the data are divided into two parts. A portion of data is utilized to obtain the parameters experimentally, while the remaining part is employed for the algorithm validation. In this algorithm, , and are chosen to be , ( ) and 10 -6 , respectively. The quantitative results of the evaluation of the proposed methods in comparison with the other mentioned methods are summarized in Table 1 and a visual outline of SINR improvement of all methods can be observed in Figure 8. DoCoBT represents the proposed method. DoCoBT (ICA) corresponds to the case, where the initial sources of the factor matrices are the components of the ICA decomposition and the second slice of the input tensor is the fECG subspace estimated by employing ICA. In DoCoBT (πCA), the se ond sli e of the tensor nd the initi l sour es re estimated through the πCA lgorithm As could e seen, the proposed method performs etter th n ICA, πCA, nd π u ker methods for a large range of SINR in values. It can also be concluded that using the results of other methods (ICA nd πCA) will improve the ur y of fECG extraction. In fact, one of the advantages of tensor decomposition is that it uses additional information alongside the main observations to obtain more accurate results. The most important conclusion is that the simultaneous usage of maternal and fetal periodicity constraints can result in the extraction of fECG with higher SINR values. The extracted fECG sources using the mentioned methods are shown in Figure 9. In order to demonstrate the performance of the proposed method for subspaces reconstruction, several channels of reconstructed subspaces are illustrated in Figures 10 and 11. In this paper, the real data are used to evaluate the generalization of the methods. The reconstructed fetal components in sensor space are represented in Figure 12. The ability of the proposed method in fetal subspaces reconstruction is obvious in Figures 10 and 12.

Conclusion
In this study, we proposed a constrained tensor decomposition method based on BT tensor decomposition to extract fECG from maternal abdominal recordings. First, the observation tensor is constructed using abdominal recordings and the fECG subspace (obtained from an existing approach, such as ICA). Then, the observation tensor is decomposed into three tensors by using the BT decomposition constrained with two constraints related to the period of the fetus and mother's ECG. In this method, the estimated fetal and maternal sources are placed in the factor matrices as the initial values of the tensor decomposition and updated in an alternative algorithm. Experimental results show that decomposition of the main tensor to three tensors and simultaneous usage of maternal and fetal quasi-periodicity constraints lead to an increase in SINR improvement in noisy conditions compared to ICA and πCA methods. In addition, the limitations of the proposed method are: (1) Slow convergence of the gradient descent algorithm, and (2) the requirement for choosing (tuning) the values of the regularization parameters and .
In order to achieve superior results, a post-processing step such as Kalman filtering can be considered to denoise the sources of the factor matrices. Also, the use of other tensor decompositions such as CP and rewriting their cost functions based on periodicity constraints could be considered as another future work.