IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 2, FEBRUARY 2011 227
Fast Technique for Noninvasive Fetal ECG Extraction
Ruben Martı´n-Clemente∗, Member, IEEE, Jose Luis Camargo-Olivares, Susana Hornillo-Mellado, Mar Elena,
and Isabel Roma´n
Abstract—This letter describes a fast and very simple algorithm
for estimating the fetal electrocardiogram (FECG). It is based on
independent component analysis, but we substitute its computa-
tionally demanding calculations for a much simpler procedure.
The resulting method consists of two steps: 1) a dimensionality re-
duction step and 2) a computationally light postprocessing stage
used to enhance the FECG signal.
Index Terms—Blind source separation, fetal electrocardiogram
extraction, independent component analysis.
I. INTRODUCTION
CURRENTLY, there is an intense research into noninvasivemethods for detecting the fetus at risk of damage or death
in the uterus. The fetal electrocardiogram (FECG) provides
useful information about the fetus’s condition: for example,
there is strong evidence that certain fetal heart rate patterns are
closely related to fetal acidemia; in addition, the duration of the
ST-segment is important in the diagnosis of fetal hypoxia, and
it has been also shown that both the QT interval and T-wave
changes are associated with fetal acidosis [1].
Given electrical signals from electrodes placed on a pregnant
woman’s abdomen, the technique of independent component
analysis (ICA) [2] can be used to separate out the FECG, which
is extremely low voltage, from the maternal ECG (MECG),
and from other unwanted background interferences, such as
the electrical activity produced by the uterus muscles. It has
been suggested that ICA outperforms most of the other signal
processing techniques [3]. However, ICA is computationally
demanding due to its use of higher-order statistics, and hence, it
is not well suited for implementation in real-time applications
nor in lightweight and portable monitors, as would be desirable
to provide mothers freedom of movement. To cope with this
problem, traditional approaches have been based on adaptive
algorithms or, more recently, on the use of a priori information
on the ECG signal [4]. The solution proposed in this letter
is one in which the computationally intensive step of ICA is
substituted for a much simpler procedure—specifically, it does
Manuscript received April 22, 2010; revised June 19, 2010; accepted July 5,
2010. Date of publication July 19, 2010; date of current version January 21,
2011. This work was supported by the “Junta de Andalucı´a” (Spain) under a
grant with reference P07-TIC-02865. Asterisk indicates corresponding author.
∗R. Martı´n-Clemente is with the Department of Signal Theory and Commu-
nications, University of Seville, Seville 41092, Spain (e-mail: ruben@us.es).
J. L. Camargo-Olivares and S. Hornillo-Mellado are with the Department of
Signal Theory and Communications, University of Seville, Seville 41092, Spain
(e-mail: jlcamargo@yahoo.es; susanah@us.es).
M. Elena is with the Department of Electronic Engineering, University of
Seville, Seville 41092, Spain (e-mail: mar@gte.us.es).
I. Roma´n is with the Department of Automatica and System Engineering,
University of Seville, Seville 41092, Spain (e-mail: isabel@trajano.us.es).
Digital Object Identifier 10.1109/TBME.2010.2059703
not require the estimation of statistics—without appreciable loss
of performance. We present experiments with real data showing
that the resulting algorithm is fast and effective.
II. PREPROCESSING
Let v1(t), . . . , vp(t) be zero-mean signals recorded from elec-
trodes placed on the mother’s body, where t ∈ Z is the discrete
time, and let v(t) be the vector whose ith component is vi(t).
The proposed preprocessing consists in reducing the number
of signals under consideration and it is intended to speed up
the estimation process. Our preprocessing is based on princi-
pal component analysis (PCA), which has the added advantage
of reducing appreciably the MECG interference [5]. Consider
that we are given N samples v(1), . . . ,v(N) of the vector sig-
nal v(t). Preprocessing starts from the data covariance matrix
Rv = (1/N)
∑N
t=1 v(t)v
T (t). This matrix can be always fac-
torized into Rv = QDQT , where
D = diag (λ1 ≥ λ2 ≥ · · · ≥ λp)
is the p× p diagonal matrix whose elements are the eigenval-
ues of Rv and Q is the matrix containing the corresponding
eigenvectors. If the MECG is strong enough, Callaerts et al. [5]
showed that the M largest eigenvalues in D are associated with
the MECG. Matrices D and Q can be then partitioned into two
groups:
D =
(
D1 0
0 D2
)
, Q = (Q1 Q2)
where D1 contains those M largest eigenvalues and the columns
of Q1 are the corresponding eigenvectors. The MECG can be
reduced significantly by projecting the data onto the subspace
spanned by Q2 [5]. Specifically, this can be written as
z(t) = D−1/22 Q
T
2 v(t) (1)
where z(t) is a (p−M)× 1 vector that contains little MECG
contribution. Matrix D−1/22 is a mere scale factor, which ensures
that the signals in z(t) are of unit variance. Of course, the
determination of M is an important problem. Some previous
works consider M = 3 [5]; however, it has been recently argued
that from M = 4 to M = 6 may be required in some cases [6].
In practice, experiments suggest finding M empirically from
the gap between the eigenvalues of the data covariance matrix.
Finally, observe that we have projected the data from a vector
space of p dimensions to a vector space of p−M dimensions.
The complete procedure can be accomplished in real time with
low computational cost by using, for example, the method in [7].
0018-9294/$26.00 © 2010 IEEE
228 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 2, FEBRUARY 2011
III. METHOD
This section presents our original contribution. After the
preprocessing, the SNR is still poor and further processing is
required. The method is based, like many ICA algorithms, on
maximizing the kurtosis (fourth-order cumulant) of the FECG
estimates. The optimization procedure is the main innovation,
and its interpretation is discussed in the context of the FECG-
estimation problem. The FECG-extraction method itself is de-
scribed in Section III-A. Mathematical justification is then pre-
sented in Section III-B.
A. Algorithm
Given z(1), . . . , z(N):
1) Define z(K) as the vector z(t) with the largest length.
2) Output y(t) = z
T (K)
‖z(K)‖ z(t).
Observe that, to speed up the process, z(K) can and should
be calculated at the same time that z(t) is calculated in the
preprocessing step.
B. Justification
Given an arbitrary vector w, it follows from the central limit
theorem that
y(t) = wT z(t) (2)
is more Gaussian when it is a sum of the FECG and the interfer-
ences than when it is equal to only one of them. In consequence,
to find w in such a way that the distribution of y(t) is as far
as possible from Gaussian seems to be a promising idea. This
general approach to the problem of “unmixing” mixed signals
is very common in ICA and is usually called maximization of
nonGaussianity [2].
The simplest measure of nonGaussianity is the kurtosis or
fourth-order autocumulant, defined by:
κy =
1
N
N∑
t=1
y4(t)− 3
N
N∑
t=1
y2(t). (3)
We maximize the kurtosis of y(t) under the unit-power
constraint
1
N
N∑
t=1
y2(t) = 1 (4)
which avoids the solution y(t) →∞. It is easily shown that this
is equivalent to constraining the norm of w to be the unity. Next,
let us deal with the maximization of κy . Traditional ICA algo-
rithms use standard procedures such as gradient cancellation.
However, as explained in Section I, we have devised a simpler
method without sacrificing much on accuracy. First, consider
the following theorem, which is stated without proof due to the
lack of space:
Theorem 1: Let {v(t), t = 1, . . . , N} be the samples of a
generic discrete-time signal. The kurtosis of v(t), defined by
κv =
1
N
N∑
t=1
v4(t)− 3
N
N∑
t=1
v2(t)
is maximized under the unit-power constraint 1N
∑N
t=1 v
2(t) =
1 by signals of the form
v∗(t) = ±
√
N ek (t)
where ek (t) is a discrete-time signal that equals 1 at t = k and
is 0 elsewhere.
To explore the vicinity of the maximum
√
N ek (t), where
k ∈ {1, . . . , N}, we perform a first-order Taylor expansion of
the kurtosis around this point:
κv ≈ N − 3 +
(
4
√
N − 6√
N
) [
v(k)−
√
N
] (5)
where N − 3 is the value of the kurtosis at the maximum. The
unit-power constraint implies that
∑N
t
=k v
2(t) = N − v2(k).
Adding [y(k)−√N ]2 , we obtain
[
v(k)−
√
N
]2 +
N∑
t
=k
v2(t) = −2
√
N
[
v(k)−
√
N
]
.
Substituting this formula into (5), we obtain after a little bit of
algebra that, for any practical value of N
κv ≈ N − 3− 2
N∑
t=1
(
v(t)−
√
N ek (t)
)2
. (6)
This formula approximates the kurtosis of the unit-power se-
quences that are in the neighborhood of
√
N ek (t). Since y(t)
is by definition of unit power, we also use this formula to ap-
proximate κy
κy ≈ N − 3− 2
N∑
t=1
(
y(t)−
√
N ek (t)
)2
. (7)
According to (7), the kurtosis κy is maximized when
N∑
t=1
(
y(t)−
√
N ek (t)
)2 (8)
is minimum—in other words, the optimum y(t) is the signal
that is as close as possible to
√
N ek (t). To determine the best
value for the index k, it is assumed that the accuracy of the
approximation (7) increases as (8) decreases. Consequently, we
minimize (8) among all possible values of k. Taking into account
that y(t) = wT z(t), a bit of algebra shows that the minimum
is obtained simply by setting
w∗ =
z(K)
‖z(K)‖ , where K = argmaxk ‖z(k)‖. (9)
By construction, y(t) is the signal that is as close as possible to
the impulse signal
√
N eK (t). Furthermore, if z(t) is periodic,
one can prove easily that y(t) is also the best approximation to an
impulse train having the same period and centered upon t = K.
Consider then the following additional interpretation: the ECG
resembles an impulse train, but the interferences degrade the
measurements and disturb this property. The algorithm restores
this property and, as a result, restores the signal itself. The
method may be then considered as a particular application of
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 2, FEBRUARY 2011 229
Fig. 1. Simulated data: (a) measured signal, (b) FECG dipole component.
the class of waveform-preserving methods for recovering ECG
signals.
Observe that this is actually an approximate optimization
procedure. Furthermore, to the best of our knowledge, this pro-
cedure is presented here for the first time. Further discussion is
illustrated with examples in the Section IV.
C. Sequential Estimation
To extract sequentially more signals, we can use the procedure
described in [2, ch. 4]. Basically, we remove y(t) from the
mixture by z′(t) = z(t)−w y(t). Then, PCA is applied again to
reduce the dimensionality in one unit. The algorithm is repeated
until all the desired signals are recovered.
IV. EXPERIMENTS AND RESULTS
A. Experiment and Results With Simulated Data
We have used the Open Source ECG Toolbox for generating
synthetic maternal and fetal ECG mixtures with realistic ECG
noises [8]. Eight simulated recordings of 5000 samples each
were generated (SNR = 10). Fig. 1 shows part of one of the sig-
nals recorded on the mother’s abdomen and one of the original
fetal cardiac dipole components [6], where both of them have
been normalized to unit-variance. In the preprocessing step, M
was set to 4. Fig. 2 represents the four elements of z(t)—observe
that z4(t) is a residue of the MECG. We have plotted a vertical
dotted line at t = K, which is the instant at which z(t) yields
its largest norm. There is a fetal R wave at t = K, as expected:
the ECG is maximum at the time instants in which there is an R
wave. We have z(K)/‖z(K)‖ = [−0.39, 0.44, 0.79,−0.15]T
and
y(t) = −0.39 z1(t) + 0.44 z2(t) + 0.79 z3(t)− 0.15 z4(t).
The larger the peak value of the R wave, the larger the contribu-
tion of zi(t) to y(t), i.e., the signals zi(t) with a more significant
presence of the FECG are more weighted. Signal y(t) is shown
in Fig. 3. Even though the preprocessing step has done “half of
the job,” in the sense that the location of the R waves is perfectly
visible, the SNR of y(t) is, at least, 2.54 dB better than that of
any of the signals zi(t), at a minimum computational cost. The
correlation coefficient between y(t) and the true FECG signal is
0.9046. To compare, FastICA, which is other ICA method [9],
produced a correlation coefficient equal to 0.9080. In both cases,
Fig. 2. Simulated data: Signals zi (t) obtained in the preprocessing step.
Fig. 3. Simulated data: Estimated fetal heartbeat signal.
Fig. 4. Real data: Abdominal recording from a pregnant woman.
before calculating the correlation, we applied a simple filtering
for the removal of the FECG baseline.
B. Experiment and Results With Real Data
Eight real cutaneous potential recordings of a pregnant
woman have been obtained from [10]. Data were recorded
at University Hospitals Leuven with the informed consent of
the patient [11]. The first five signals correspond to electrodes
placed on the woman’s abdominal region. The last three signals
correspond to the electrodes located on the mother’s thoracic re-
gion. The sampling rate is 250 Hz with 12-bit resolution. Fig. 4
shows a portion of one of the abdominal’s recordings. Even
though the FECG is much weaker than the MECG, it is slightly
visible. Parameter M was set to 4. Following the procedure of
Section III-C, we extracted sequentially the four signals shown
in Fig. 5. The first and the second one from the top correspond to
the FECG, since the heart rate is about twice the maternal one.
The first signal is also displayed in Fig. 6(a). According to the
most usual interpretation of ICA [3], these signals are estimates
of the cardiac dipole vector components [6]. To validate the
method, we repeated the experiment using FastICA. Fig. 6(b)
230 IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 58, NO. 2, FEBRUARY 2011
Fig. 5. Real data: Output signals generated sequentially by the algorithm.
Fig. 6. Real data: (a) FECG estimated by the proposed method, (b) FECG
estimated by FastICA, (c) signal (a) after postprocessing.
shows the FECG estimated by FastICA: there is a strong similar-
ity with the results obtained by our algorithm. To show that the
algorithm retains enough information on the FECG, in Fig. 6(c),
we show a simple postprocessing of the signal generated by our
method using the filter described in [12]. Observe that the loca-
tion of the P and T waves has been improved significantly.
We also estimated the FECG using other ICA meth-
ods: π-ICA [4], which is a powerful ICA algorithm specif-
ically devised for the FECG estimation problem, JADE
[2], and Pearson ICA [13]. Similarity between the esti-
mates was measured by the correlation coefficient ρi =∑N
t=1 (y(t) yi(t))/
√∑N
t=1 y
2(t)
∑N
t=1 y
2
i (t), where y(t) is
the best FECG estimate obtained by our algorithm, and yi(t)
represents the best FECG estimate obtained by: FastICA (i = 1),
π-CA (i = 2), JADE (i = 3), and Pearson ICA (i = 4). A mea-
sure of the execution time for generating all the estimates was
defined by λi = Ti/T , where T is the time required by our
TABLE I
COMPARISON BETWEEN FIGURES OF MERIT FOR EACH ICA METHOD
proposed method (0.01 s) and Ti represent the mean time re-
quired by the others, averaged over ten experiments (the actual
time may depend on the initialization point). Here, index i has
the same meaning as mentioned earlier. Results are shown in
Table I. In addition, to reduce further the baseline wander and the
high-frequency noise, we filtered all the estimated fetal heart-
beat signals with a simple bandpass filter with cutoff frequencies
at 1 and 100 Hz. After this postprocessing, the correlation coef-
ficients mentioned earlier increased to 0.9492, 0.9277, 0.9491,
and 0.9647, respectively.
V. CONCLUSION
The algorithm is simple and fast. Execution time should be
considered here as a measure of computational complexity. Even
though other methods are also able to operate in real time,
they carry out very complex calculations that result in increased
execution time. From this viewpoint, the proposed algorithm is
advantageous for the design of battery-powered devices, since
it is usually admitted that computational complexity is related
to energy consumption.
Results are comparable to those obtained by other more com-
plex ICA-based approaches. However, note that only R–R in-
tervals can be reliably determined in general at the current state
of the art.
REFERENCES
[1] M. Taylor, M. Thomas, M. Smith, S. Oseku, N. Fisk, A. Green, S. Paterson,
and H. Gardiner, “Non-invasive intrap. fetal ECG: Preliminary report,” Br.
J. Obstet. Gynaecol., vol. 112, pp. 1016–1021, 2005.
[2] A. Cichocki and S.-I. Amari, Adaptive Blind Signal and Image Processing:
Learning Algorithms and Applications. New York: Wiley, 2002.
[3] V. Zarzoso and A. K. Nandi, “Noninvasive fetal electrocardiogram extrac-
tion: Blind separation versus adaptive noise cancelation,” IEEE Trans.
Biomed. Eng., vol. 48, no. 1, pp. 12–18, Jan. 2001.
[4] R. Sameni, C. Jutten, and M. Shamsollahi, “Multichannel electrocardio-
gram decomposition using periodic component analysis,” IEEE Trans.
Biomed. Eng., vol. 55, no. 8, pp. 1935–1940, Aug. 2008.
[5] D. Callaerts, B. D. Moor, J. Vandewalle, and W. Sansen, “Comparison of
SVD methods to extract the foetal ECG from cutaneous electrode signals,”
Med. Biol. Eng. Comput., vol. 28, pp. 217–224, 1990.
[6] R. Sameni, G. Clifford, C. Jutten, and M. Shamsollahi, “Multichannel
ECG and noise modeling: Application to maternal and fetal ECG signals,”
EURASIP J. Adv. Signal Process., vol. 2007, art. 43407, 2007.
[7] D. Han, Y. N. Rao, J. C. Principe, and K. Gugel, “Real-time PCA (principal
component analysis) implementation on DSP,” in Proc. IEEE Int. Conf.
Neural Netw., Jul. 25–29, 2004, pp. 2159–2162.
[8] (2010). [Online]. Available: http://ee.sharif.edu/ecg/
[9] A. Hyva¨rinen, “Fast and robust fixed-point algorithms for independent
component analysis,” IEEE Trans. Neural Netw., vol. 10, no. 3, pp. 626–
634, May 1999.
[10] (2010). [Online]. Available: ftp://ftp.esat.kuleuven.ac.be/pub/SISTA/
delathauwer/data/foetal_ecg.dat
[11] L. de Lathauwer, private communication, Jul. 2010.
[12] R. Sameni, M. Shamsollahi, and C. Jutten, “Model-based bayesian filtering
of cardiac contaminants from biomedical recordings,” Physiol. Meas.,
vol. 29, pp. 595–613, 2008.
[13] A. Gupta, M. C. Srivastava, V. Khandelwal, and A. Gupta, “A novel
approach to fetal ECG extraction and enhancement using blind source
separation and adaptive fetal ECG enhancer,” in Proc. Int. Conf. Inform.,
Comm. Signal, 2007, pp. 1–4.