Fetal motion estimation from noninvasive cardiac signal recordings

The pre-processing consists of baseline wander removal and power-line cancellation (whenever applicable). We use general methods and tools previously developed in Sameni (2008) and (2014). For baseline wander removal, a high-pass zero-phase forward-backward infinite impulse response (IIR) filter is used with a  −3 dB cutoff frequency of 0.5 Hz. For datasets with a stationary power-line noise, a standard second-order IIR notch filter is used. For datasets contaminated by non-stationary power-line artifacts (varying power-line amplitude), the robust Kalman filtering scheme developed in Sameni (2012) and (2014) is used.

Herein, a variant of the deflation procedure proposed in Sameni et al (2010) is used for maternal ECG/MCG cancellation. This method can separate the fetal cardiac signal subspaces from the maternal cardiac signal and background noise, in low SNR and degenerate situations (few number of abdominal electrodes), without reducing the data dimensions.

The algorithm consists of three stages: multichannel decomposition, channel-wise denoising, and multichannel re-composition. In the first stage, periodic component analysis (πCA) is used (Sameni et al 2008), which linearly transforms multichannel signals into the same number of channels ranked in order of resemblance with the maternal heart rhythm. Next, the first few channels of the ranked signals are denoised using a single-channel denoiser (e.g. wavelet denoiser or Kalman filter (Sameni et al 2007b)). In the third stage, the denoised signals are back-projected to the original subspace using the inverse of the linear transform of the first stage. These three stages are repeated several times until all components corresponding to the maternal ECG/MCG are removed. The advantage of this scheme is that it overall acts as a nonlinear denoiser, which eliminates the maternal subspace stage-by-stage, without omitting the fetal subspace even in low-rank scenarios (Sameni et al 2010). Moreover, the resulting signals remain in the original sensor space, and not in the transform domain. Therefore, the method is even applicable to relatively low number of channels.

To this point, the maternal artifacts have been removed. The following stages of the algorithm can either be applied in the lead space (after maternal ECG/MCG cancellation), or over a linear transformed version of these signals. The objective of this (optional) transform is to enhance the signal quality by transforming (or even projecting) the data to the null-space of the multi-channel noise.

In a blind scenario, where no prior knowledge regarding the fetal ECG/MCG is available, blind source separation (BSS) methods such as independent component analysis (ICA) combined with principal component analysis (PCA) dimension reduction can be be used at this stage (Comon and Jutten 2010). In semi-blind scenarios, where prior knowledge such as the fetal R-peak locations are available, semi-blind source separation (SBSS) methods such as πCA can be used to linearly transform the fetal signals to a domain in which the channels are ranked in order of resemblance with the FHR (Sameni et al 2008). SBSS is apparently more effective than BSS; but requires fetal R-peak detection in advance. In either case (blind or semi-blind), in order to preserve the relative amplitudes and fetal rotation angles, this stage is performed using a global separation matrix over the entire dataset, and not using block-wise demixing matrices, which are commonly used in online source separation schemes.

After removing the maternal cardiac interferences, the most dominant pseudo-periodic peaks correspond to the fetal QRS. However, depending on the sensor and fetal positioning, the different abdominal channels do not equally show the fetal cardiac activity. In practice, only a subset of the channels contain traces of the fetal QRS, which should be selected by visual inspection or automatic channel ranking. In the automatic scheme, after removing the maternal cardiac interferences and quality enhancement, the channels are ranked in order of quality and resemblance with the fetal QRS peaks. Several automatic ECG quality assessment methods have been developed in previous studies (Niknam 2012, Behar et al 2013). We use a computationally efficient approach based on matched filtering with predefined templates. The method is based on the fact that due to the small fetal heart size and its distance from the abdominal leads, fetal cardiac signals are recorded at far-field (see Geselowitz (1989) for a discussion on cardiac signal propagation models). As a result, fetal PQRST complexes are not morphologically as diverse and complex as adult signals. In figure 2, the three most common fetal QRS morphologies are shown, which have been synthetically constructed by mathematical functions. The idea is inspired from similar methods previously used for ECG waveform modeling using Gaussian kernels (Sornmo et al 1981, Laguna et al 1996, McSharry et al 2003, Sameni et al 2007a). We use these samples as templates for matched filters, denoted by h_i(n) ($i=1,\cdots,K$), where K is the number of predefined templates.

Zoom In Zoom Out Reset image size

For real-time FHR and FM monitoring systems, the sample fetal beat templates may alternatively be selected manually or updated adaptively over time, which is beyond the scope of the current study.

By passing each data channel x_k(n) (after maternal ECG/MCG cancellation) through the matched filter h_i(n), the normalized cross-correlation and the vectorcardiogram amplitude (VA) are respectively calculated as follows:

$$\begin{eqnarray}&&{{r}_{ki}}(n)=\frac{{\sum\nolimits_{m}{{h}_{i}}(m){{x}_{k}}(n-m)}}{{{\left[\sum\nolimits_{m}{{h}_{i}}{{(m)}^{2}}\right]}^{1/2}}}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{r}_{k}}(n)={{\left[\underset{i=1}{\overset{K}{\sum}}{{r}_{ki}}{{(n)}^{2}}\right]}^{1/2}}\end{eqnarray} \tag{ 2 }$$

VA is a measure of the overall strength of the vectorcardiogram at the matched filter outputs (Sameni 2008, Schmidt et al 2014)¹.

The signal r_k(n) is next used as the reference signal for fetal R-peak detection by a local peak search algorithm over a sliding window, previously developed in Sameni (2008) and (2014). For the studied datasets, the sliding window length was selected to be S  =  375 ms, which enables the detection of the FHR up to 160 beats per minute (BPM), considered as the upper bound of normal FHR (American College of Obstetricians and Gynecologists and others 2009). After detecting the fetal R-peaks from r_k(n), the robust weighted average (RWA) (Leski 2002, Leski and Gacek 2004), of all x_k(n) channels are calculated. RWA uses the detected R-peaks to align the ECG beats of each channel k in a matrix ${{\mathbf{X}}_{k}}\in {{\mathbb{R}}^{B\times D}}$; where B is the total number of beats and D is the window length around each R-peak (D  =  40 ms in the later reported results). Next, the RWA vector signal ${{\mathbf{m}}_{k}}\in {{\mathbb{R}}^{1\times D}}$ and the sample-wise variance around the RWA $\sigma _{k}^{2}$ are calculated. Intuitively, a 'good quality' fetal ECG/MCG channel has a small variance around its RWA. Therefore, the better the signal quality, the smaller the variance should be around its RWA.

The overall fetal R-peak detection and channel ranking algorithm is summarized in algorithm 1.

Algorithm 1. Fetal R-peak detection and channel ranking algorithm.

Require: Multichannel signal xk(n) ($k=1,\cdots,M$), after maternal ECG/MCG cancellation

1: Calculate matched filter outputs rki(n) from (1)

2: Calculate VA rk(n) from (2)

3: Detect local peaks of rk(n) in a sliding window of length S to construct the HR sequence pi ($i=1,\cdots,B$)

4: Construct the R-peak time-aligned beat matrix ${{\mathbf{X}}_{k}}\in {{\mathbb{R}}^{B\times D}}$ for each channel k

5: Calculate the RWA ${{\mathbf{m}}_{k}}=\text{RWA}\left({{\mathbf{X}}_{k}}\right)\in {{\mathbb{R}}^{1\times D}}$

6: Calculate the sample-wise variance around the RWA $\sigma _{k}^{2}=E\{{{\left({{\mathbf{X}}_{k}}-\mathbf{1}{{\mathbf{m}}_{k}}\right)}^{2}}\}$, where ${{\mathbf{1}}^{B\times 1}}$ is a vector of ones, and $E\left\{\cdot \right\}$ denotes averaging over both dimensions.

7: Rank the channels in ascending order of $\sigma _{k}^{2}$.

The channel ranking can be applied either in the lead space or in the transform space (as discussed in section 2.3.). In order to avoid using the channels without any fetal ECG/MCG contribution, the first M channels with minimal RWA variance—as described in algorithm 1—are used for fetal QRS-loop alignment.

A spatio-temporal alignment of QRS-complex VCG loops was proposed in Astrom et al (2000) and Sörnmo (1998) to compensate motion-induced QRS variations, such as respiratory-related ECG variations. Herein, the same method is used for estimating the FM from fetal signals after maternal ECG/MCG cancellation and FHR detection. In the application of interest, the major source of fetal ECG/MCG morphological variations is the FM. Although the original VCG loop alignment method was applied to the Frank lead system, we can apply this method to different datasets, in which the leads are not necessarily orthogonal; but the sensor positions are fixed with respect to the maternal body surface.

Consider a multichannel QRS complex $\mathbf{Z}\in {{\mathbb{R}}^{M\times W}}$ obtained from a fixed window of length W (30 ms in the following results), around a typical fetal R-peak. We would like to align $\mathbf{Z}$ with a reference QRS complex ${{\mathbf{Z}}_{R}}\in {{\mathbb{R}}^{M\times \left(W+2 \Delta +1\right)}}$, with the integer $\Delta$ defined in the sequel. We assume that $\mathbf{Z}$ has been altered by several transformations over ${{\mathbf{Z}}_{R}}$, due to rotations, scaling and shifts of the fetal heart signals (due to the fetal rotations and movements). Therefore, a data model similar to the one proposed in Astrom et al (2000) and Sörnmo (1998) is adapted for FM estimation, as follows.

• Scaling: The movement of fetal body, fetal extremities and maternal respiration can lead into variations of the maternal body volume conductor, resulting in amplitude variations in the fetal ECG/MCG recorded on the body surface. We can model such amplitude variations with a positive-valued scalar parameter α, used for defining the actogram² and accounts for fetal QRS-complex amplitude variations (figure 3).
• Rotation: Fetal body rotations change the angles of the fetal body with respect to the maternal body axes and the surface leads (figure 3). Such rotations can be modeled by a rotation matrix, which maps ${{\mathbf{Z}}_{R}}$ to $\mathbf{Z}$. The rotational changes of the QRS-complex are modeled by an orthogonal matrix $\mathbf{Q}\in {{\mathbb{R}}^{M\times M}}$, which can be decomposed into M independent rotation angles using the Givens transform (Golub and van Loan 1996). However, since we are interested in gross fetal rotations due to the fetal body movements, and not the beat-wise rotations considered in Astrom et al (2000) and Sörnmo (1998), due to adult respiration and baseline wander, only the principal subspace angle between ${{\mathbf{Z}}_{R}}$ and $\mathbf{Z}$ is calculated and a smoothing procedure is later applied on the final FM results to obtain gross rotations.
• Time synchronization: It is initially assumed that each $\mathbf{Z}$ is well synchronized in time with ${{\mathbf{Z}}_{R}}$, i.e. the fetal R-peak detection was exact. However, in order to compensate minor R-peak mis-detections, a shift matrix ${{\mathbf{J}}_{\tau}}\in {{\mathbb{R}}^{\left(W+2 \Delta +1\right)\times W}}$ is also considered, which permits minor peak displacements over a finite number of integer delays $\tau \in \left[- \Delta, \Delta \right]$. The shift matrix ${{\mathbf{J}}_{\tau}}$ is a $W\times W$ identity matrix shifted up (down), depending on the sign of τ, and filled with zero rows from the bottom (top) to keep its size $\left(W+2 \Delta +1\right)\times W$ (see Astrom et al (2000)).

Zoom In Zoom Out Reset image size

The above transformations are merged together as follows:

$$\begin{eqnarray}&&\mathbf{Z}=\alpha \mathbf{Q}{{\mathbf{Z}}_{R}}{{\mathbf{J}}_{\tau}}+\mathbf{V}\end{eqnarray} \tag{ 3 }$$

where $\mathbf{V}\in {{\mathbb{R}}^{M\times W}}$ is additive white noise. With this data model, following Astrom et al (2000), we seek the least-squares estimation of the parameters α, $\mathbf{Q}$, and ${{\mathbf{J}}_{\tau}}$, for each fetal QRS window, i.e.

$$\begin{eqnarray}&&\left\{{{\alpha}^{\ast}},{{\mathbf{Q}}^{\ast}},\mathbf{J}_{\tau}^{\ast}\right\}=\underset{\alpha,\mathbf{Q},{{\mathbf{J}}_{\tau}}}{{\arg \min}}\Vert \mathbf{Z}-\alpha \mathbf{Q}{{\mathbf{Z}}_{R}}{{\mathbf{J}}_{\tau}}\Vert _{F}^{2}\end{eqnarray} \tag{ 4 }$$

This is an optimization problem previously solved in Astrom et al (2000) using iterations of singular value decomposition (SVD) of the following matrix, over fixed values of the shift τ:

$$\begin{eqnarray}&&\mathbf{D}\overset{\Delta}{{=}}\mathbf{ZJ}_{\tau}^{T}\mathbf{Z}_{R}^{T}\end{eqnarray} \tag{ 5 }$$

Denoting the SVD of $\mathbf{D}$ as $\mathbf{D}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}$, the rotation matrix ${{\mathbf{Q}}_{\tau}}$ and scale factor ${{\alpha}_{\tau}}$ are found as follows

$$\begin{eqnarray}&&{{\mathbf{Q}}_{\tau}}=\mathbf{UV}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\alpha}_{\tau}}=\frac{\text{tr}\left({{\mathbf{Z}}^{T}}{{\mathbf{Q}}_{\tau}}{{\mathbf{Z}}_{R}}{{\mathbf{J}}_{\tau}}\right)}{\text{tr}\left(\mathbf{J}_{\tau}^{T}\mathbf{Z}_{R}^{T}{{\mathbf{Z}}_{R}}{{\mathbf{J}}_{\tau}}\right)}\end{eqnarray} \tag{ 7 }$$

These calculations are repeated for all $\tau \in \left[- \Delta, \Delta \right]$ and the optimal set $\left({{\alpha}^{\ast}},{{\mathbf{Q}}^{\ast}},\mathbf{J}_{\tau}^{\ast}\right)$ is selected as the one, which minimize the argument of (4).

The reference signal ${{\mathbf{Z}}_{R}}$ can be any arbitrary clean QRS complex, which is selected by visual inspection or through an automatic procedure. For the automatic selection we propose to calculate the RWA of all $\mathbf{Z}$ using the detected fetal beats. Next, we calculate the error between the RWA and all $\mathbf{Z}$. The beat with the minimum error energy is chosen as the reference beat ${{\mathbf{Z}}_{R}}$. For highly noisy data in which the probability of missing or false fetal QRS complexes increases, the RWA calculation can be performed twice, over a subset of beats having the least error energy in the first run. This makes the selection of ${{\mathbf{Z}}_{R}}$ less susceptible to mis-detected QRS complexes.

The procedure detailed in section 2.5 is repeated for each of the detected fetal QRS complexes using the reference beat ${{\mathbf{Z}}_{R}}$, as detailed above. We define the beat-wise parameter ${{\alpha}^{\ast}}$ as the fetal actogram, accounting for the fetal extremities movement. We further coin the variations of the principal subspace angles of ${{\mathbf{Q}}^{\ast}}$, as the fetal rotatogram, which accounts for fetal movements that change the angle of the fetal body axes with respect to the maternal body. The notion of principal angle between subspaces (Golub and van Loan (1996), chapter 12), Knyazev and Argentati (2002), used for tracking the rotations of the fetus is summarized in the appendix. The procedure consists of tracking the principal angle between the beat-wise rotation matrix ${{\mathbf{Q}}^{\ast}}$ and an arbitrary reference orthogonal matrix over time. For simplicity, the reference matrix can be an identity matrix of the same size, which simplifies the principal angle calculations. A similar idea has been used in Schmidt et al (2014).

It should be noted that for M  =  1 (only one fetal related ECG/MCG channel) only ${{\alpha}^{\ast}}$ may be calculated and the rotation matrix ${{\mathbf{Q}}^{\ast}}$ no longer exists.

To this point, the fetal actogram and rotatogram have been calculated beat-by-beat, without considering any inter-beat dynamics. In reality, fetal movements are rather slow as compared with the fetal heart rate, which is typically above 110BPM (American College of Obstetricians and Gynecologists and others 2009). Therefore, true variations in the actogram and rotatogram (not the ones due to noise, mis-detection of fetal R-peaks, or casual abnormal beats), are expected to follow rather smooth dynamics in time. This motivates the application of post-processing for smoothing the FHR, actogram and rotatogram. Various methods can be used for this purpose, ranging from Kalman filters to Tikhonov regularization (Tarvainen et al 2002). Among the three parameters, the FHR is less susceptible to noise as it only requires a robust peak detection, described in previous sections. Moreover, for heart rate variability (HRV) analysis, it is better not to filter the FHR as much as possible, to avoid the manipulation of its temporal and spectral properties. Therefore, for the FHR, we use a variant of nonlinear trimmed median filter (Arce 2004). In this method, we calculate the median over a sliding window of length T and replace the FHR samples, which have an absolute difference greater than a predefined threshold th, from the window median, with the window median. The advantage of this filter is that it only acts over the outlier FHR samples and keeps the calculated FHR unchanged in normal epochs. This makes the algorithm more suitable for HRV analysis, which rely on minor fluctuations of the FHR. On the other hand, a moving average filter with the same sliding window length T is used for the actogram and rotatogram. All post processing filters are applied in a non-causal smoothing scheme, to guarantee zero-phase delay between the raw and denoised signals and avoid any phase delays, which might influence the causality analysis, which follows.

Besides ultrasound imaging, the standard procedure of electronic fetal monitoring (EFM) is currently based on the FHR and its variations (Vandenberghe 1993). In American College of Obstetricians and Gynecologists and others (2009), a three-tiered procedure was proposed for categorizing normal, abnormal and intermediate FHR. Accordingly, normal and abnormal epochs are identified from the existence or lack of baseline FHR, FHR accelerations and/or decelerations. Herein, we use the coherency between the baseline-removed actogram/ rotatogram and the FHR, to assess the performance of the proposed FM calculation methods.

The coherency between fetal movements and FHR accelerations has been reported in recent studies (Schmidt et al 2014). Additionally, in our studied datasets, it was observed that there is significant coherency between the onsets of FHR, actogram and rotatogram accelerations/decelerations; although the duration of the FHR accelerations/decelerations were commonly longer than the active epochs of the motion-based indexes (actogram/rotatogram). For illustration, a sample segment of the baseline-removed and the normalized actogram and rotatogram parameters are shown in figure 4.

Zoom In Zoom Out Reset image size

Besides the coherency, the study of causal relationships between the onsets of the proposed indexes is also of interest. We are interested to know whether the FHR accelerations/decelerations are initiated by the fetal motion, or vice-verse. The answer to this question can lead to the notion of motion-induced FHR acceleration/deceleration. The hypothesis is that as in adults, some of the variations of the FHR can be due to the fetal motion, rather than spontaneous behaviors of the fetal sympathetic and parasympathetic brain systems. In other contexts, methods such as the Granger causality test have been used (Granger 1988). In what follows, we propose a novel test for causal relation assessment between the proposed parameters.

To verify the causality hypothesis, first, epochs of considerable fetal motion and FHR variations are calculated by thresholding the energy envelopes of the FHR/actogram/rotatogram signals around their baseline (see American College of Obstetricians and Gynecologists and others (2009) and Macones et al (2008) for the definition of the FHR baseline). The local energy of a signal x(t) around the baseline b_x(t) in a sliding window of length w is defined as follows:

$$\begin{eqnarray}&&{{p}_{x}}(t)=\frac{1}{w}\underset{\tau =-w/2}{\overset{w/2}{\sum}}{{\left(x\left(t-\tau \right)-{{b}_{x}}\left(t-\tau \right)\right)}^{2}}\end{eqnarray} \tag{ 8 }$$

where the subscript $x\in \left\{\text{FHR},\text{ACT},\text{ROT}\right\}$, refers to the three possible types of signals: FHR, acogram (ACT), or rotatogram (ROT), and b_x(t) is the low frequency baseline of x(t). In the reported results we use w  =  25 samples, which was found to be suitable to detect even the steepest accelerations and decelerations of the fetal motions: 'baseline to peak in less than 30 s' (American College of Obstetricians and Gynecologists and others (2009), table 1).

Next, activity pulses of the FHR/actogram/rotatogram, are reported when p_x(t) exceeds some predefined threshold T_x (the algorithm later sweeps over different values of this threshold):

$$\begin{eqnarray}&&{{a}_{x}}(t)\overset{\Delta}{{=}}\left\{\begin{array}{*{35}{l}} 1\quad \text{if}~{{p}_{x}}(t)\geqslant {{T}_{x}} \\ 0\quad \text{otherwise} \end{array}\right. \end{eqnarray} \tag{ 9 }$$

The overlap, vicinity, leading or lagging of the activity pulses is next used to find coherency and causal relationships between the different parameters. For this, the activity onset sequence $\left\{t_{x}^{1},t_{x}^{2},\cdots,t_{x}^{p}\right\}$ and the activity offset sequence $\left\{u_{x}^{1},u_{x}^{2},\cdots,u_{x}^{p}\right\}$ are respectively defined as the rising and falling edges of the a_x(t) pulses.

For two signals x(t) and y(t) ($x,y\in \left\{\text{FHR},\text{ACT},\text{ROT}\right\}$), $\alpha _{xy}^{ij}\overset{\Delta}{{=}}t_{x}^{i}-t_{y}^{j}$ and $\beta _{xy}^{ij}\overset{\Delta}{{=}}t_{x}^{i}-u_{y}^{j}$ are respectively defined as the rise-to-rise and rise-to-fall times between the ith and jth pulses of a_x(t) and a_y(t). The activity onsets $t_{x}^{i}$ and $t_{y}^{j}$ are said to be coherent (or overlapping), if $|\alpha _{xy}^{ij}|\leqslant \epsilon$, where ε is the coherency window length. Moreover, if $0<\alpha _{xy}^{ij}\leqslant \epsilon$ the activity onset $t_{y}^{j}$ is leading $t_{x}^{i}$, and if $-\epsilon \leqslant \alpha _{xy}^{ij}<0$ the activity pulse $t_{y}^{j}$ is lagging $t_{x}^{i}$. In the former case (when x is leading y), a rise-to-rise counter $L_{xy}^{RR}$ is incremented and in the latter case (when y is leading x), the rise-to-rise counter is decremented. Similar definitions are proposed between the onset $t_{x}^{i}$ and the offset $u_{y}^{j}$ to consider the rise-to-fall instances counted by a rise-to-fall counter $L_{xy}^{RF}$.

If the two events x and y were independent, the probability of successive rising or falling edges of their activity pulses would be equal to 0.5 ($L_{xy}^{RF}\approx 0$). However, when the two events are correlated, $L_{xy}^{RF}$ tends to positive or negative values. At specific activity thresholds T_x and T_y, for a significant causal relationship between x(t) and y(t), we expect the rise-to-rise counts to be more frequent than their rise-to-fall counts. Therefore, the signal x(t) is assumed to be the cause of y(t) on average if $L_{xy}^{RR}>0$, i.e. the total number of leading coherencies of x(t) exceeds the total number of leading coherencies of y(t) at various thresholds T_x and T_y. In practice, the thresholds T_x and T_y can be swept over L (a finite number of) levels ranging from the median up to the maximum of p_x(t) and p_y(t), to make the algorithm scale invariant and more robust to irrelevant and minor fluctuations of p_x(t) and p_y(t) around their minima.

Finally, we use these measures to define a normalized causality index

$$\begin{eqnarray}&&{{C}_{x\to y}}=\frac{L_{xy}^{RR}-L_{xy}^{RF}}{L_{xy}^{RR}+L_{xy}^{RF}}\end{eqnarray} \tag{ 10 }$$

and an average x−y lead-time index:

$$\begin{eqnarray}&&{{\bar{\alpha}}_{xy}}=E\{\alpha _{xy}^{ij}\}\end{eqnarray} \tag{ 11 }$$

where $E\{\cdot \}$ is the averaging operator applied over all successive pulse pairs of x and y. The proposed causality index is $-1\leqslant {{C}_{x\to y}}\leqslant 1$. Positive values of ${{C}_{x\to y}}$ indicate that x(t) is the cause of y(t) and negative values indicate that y(t) is the cause of x(t).

For illustration, the proposed causality index has been tested on random pulses and sinusoidal signals with random leads and lags. Figures 5 and 6 show the typical waveforms used for this simulation and the causality test results of 100 trials with random leads and lags. The random pulses in figure 5(a) have been obtained by generating random spikes and convolving the spike sequence with bell-shaped waveforms of arbitrary height and width (the signal x), plus additive noise. The second sequence (y) was obtained by convolving a randomly shift version of the same spike sequence with bell-shaped waveforms of arbitrary height and width. The sine-wave samples in figure 5(b) consist of sinusoids of arbitrary amplitude and frequency with random phase lags between x and y and under additive white noise. For both sample cases, it is seen that the sign of ${{C}_{x\to y}}$ is always in accordance with the signal's lead or lag property, despite the different amplitudes and widths of the components of x and y. In various simulations performed on similar synthetic data, this coherency existed even in rather low SNR.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first database contains noninvasive fetal biomagnetic recordings of 39 pregnant woman using 156 primary channels, with a sampling rate of 610.3516 Hz during 6 minutes. The datasets were recorded with the SARA (Squid Array for Reproductive Assessment, VSM MedTech, Ltd.) installed at the University of Tübingen (Kiefer-Schmidt et al 2013). This database is referred to as TÜBINGEN in the sequel.

This database has many channels and a very high quality (as compared with typical fetal ECG recordings). Therefore, the maternal MCG was eliminated by a single stage of the deflation algorithm described in figure 1, followed by fetal channel sorting and extraction. After extracting the fetal MCG, the channels were ranked according to their resemblance with the fetal MCG. Next, the FHR, actogram, and rotatogram were calculated according to the proposed scheme over the entire 6 min. A typical case is shown in figure 7. Accordingly, there are very few outliers in the calculated FHR (across more than 900 beats) and by visual inspection, there is significant coherency between onset pulses of the extracted parameters.

Zoom In Zoom Out Reset image size

The second database is the Noninvasive Fetal ECG Database (NIFECG) with 55 subjects PhysioNet (2007), consisting of two thoracic and three to four abdominal channels sampled at 1 kHz, with a duration of 4.5–16.5 min. For proof of concept, the proposed method was applied to 16 subjects of this database. Due to the few number of abdominal leads, five iterations of the deflation algorithm was used for maternal ECG cancellation without any dimension reduction. The FastICA library (Gävert et al 2005), was then applied to the residual signals to extract the fetal ECG, followed by fetal ECG ranking to sort the channels according to their resemblance with the fetal ECG. The fetal R-peaks and the FHR were obtained from the first-rank channel, followed by actogram and rotatogram calculations. A typical sample of the first channel used for FHR calculation, and the obtained FHR, actogram and rotatogram are shown in figure 8. Again the FHR outliers are very few across more than 2500 beats, for this subject.

Zoom In Zoom Out Reset image size

The next database is the Abdominal and Direct Fetal ECG Database (ADFECG) PhysioNet (2012), recorded from five subjects during labor, consisting of four abdominal and a single fetal scalp lead sampled at 1 kHz, with a duration of 5 min. Four out of five subjects were used for this study. Since this database contains a dedicated fetal scalp lead, the fetal R-peaks can be directly calculated from the scalp lead. Therefore, the πCA developed in Sameni et al (2008), which uses both the maternal and fetal R-peaks to separate the fetal and maternal ECG can be used in the deflation algorithm. Hence, the maternal ECG was removed in a single stage of the deflation algorithm. The proceeding stages were similar to the ones used for the NIFECG database. Figure 9 shows typical FHR, actogram and rotatogram results obtained from the ADFECG database.

Zoom In Zoom Out Reset image size

The final sample is a twin fetal MCG recording with a duration of 30 min and 195 channels, sampled at 1025 Hz from the Jena Fetal Monitoring Database (Biomagnetic Center/Department of Neurology, Department of Obstetrics, Jena University Hospital). The recording details are given in Hoyer et al (2013) and Schmidt et al (2014). This dataset is referred to as JENA in the sequel.

Due to the high quality of the channels, the maternal cancellation procedure was similar to the TÜBINGEN database. Next, in order to separate the two fetuses, a single stage of ICA was used for BSS and quality enhancement, as shown in figure 1. The FastICA library was used at this stage (Gävert et al 2005), with a PCA dimension reduction to project the 195 channels to a 50 dimensional subspace (for reducing the processing load of the following stages). The channel ranking was again performed automatically using the method described in section 2.4. However, the association of the obtained MCG traces to each fetus were done by visual inspection (since only one subject was studied in this case). Sample segments of the twin MCG are shown in figure 10(a) (longer traces of the same twins were previously reported in (Sameni (2008), chapter 7). After obtaining the MCG signals of each of the twins, the procedure of FHR, actogram, and rotatogram extraction was performed for each of the twins separately. The results are shown in figure 10. The signal quality of one of the twins was better in this case (perhaps due to its position with respect to the surface SQUID leads), resulting to more robust FHR, actogram, and rotatogram parameters.

Zoom In Zoom Out Reset image size

The causal relationship between the FHR, actogram, and rotatogram was calculated for all four databases, using the proposed index. As mentioned before, the sign of the causality index can be used as an indication of lead or lag property between two parameters. According to figure 11, in 34 out of 60 subjects (56.6%), the actogram was ahead of the FHR, in 30 out of 60 subjects (50%), the rotatogram was ahead of the FHR and in 43 out of 60 subjects (71.7%), either the actogram or rotatogram (one of the motion-induced parameters) were ahead of the FHR, which approves the motion-induced FHR acceleration/deceleration hypothesis for the these subjects.

Zoom In Zoom Out Reset image size

The causality time between the FHR, actogram, and rotatogram were also calculated for all subjects of four databases, using the pulses defined in (9) for L  =  22 different thresholds ranging between the median and maximum value of each signal. A coherency window length of $\epsilon =35$ samples was found to be optimal over all datasets. The mean and per-subject lead-times are reported in table 1 and figure 12, respectively.

Zoom In Zoom Out Reset image size