Abstract
Fetal heart rate extraction from the abdominal ECG is of great importance due to the information that carries in assessing appropriately the fetus well-being during pregnancy. In this paper, we describe a method to suppress the maternal signal and noise contamination to discover the fetal signal in a single-lead fetal ECG recordings. We use a locally linear phase space projection technique which has been used for noise reduction in deterministically chaotic signals. Henceforth, this method is capable of extracting fetal signal even when noise and fetal component are of comparable amplitude. The result is much better if the noise is much smaller (P wave and T wave can be discovered).
1. Introduction
Since the early work of Cremer in 1906, various methods for fetal monitoring have been proposed to obtain information about the heart status. The cardiac electrical activity of a fetus can be recorded noninvasively from electrodes on the mother’s body surface. Such recordings of fetal electrocardiograms (ECGs) are complicated by the existence of the mother’s ECG and effectively random contaminations due to noncardiac sources. Furthermore, the fetal signal is rather small due to the size of the fetal heart and the intervening tissue. We have to separate the fetal ECG (FECG) from the maternal trace and from the other contaminations. The ECGs is the tool for the clinical diagnostic because the nonlinear chemilogical excitation of cardiac tissue and signal show both fluctuation and remarkable structure. Moreover, because the length of cardiac cycle which is measured by the distance of two successive QRS spike fluctuates with the predictable component, ECGs is deterministic chaotic. Besides, ECGs also comprise the excitation of the mother, the distortion of tissues on the transmission.
In general, linear filter is used to separate signals based on their difference in frequency domain which can be expressed in terms of Fourier spectrum. However, even the optimal linear filter, the Wiener filter, cannot be successful in this case because ECGs of the pregnant include both maternal and fetal signals which have the same spectral contents, and the noise coming from the electric equipments has a broad band and random [1].
Another solution is nonlinear filter which offers some superior features to linear projection in this case. However, because the method is based on theory of deterministic dynamical system, we must make sure that the observed data contents the typical properties of deterministic chaotic signals. We can find it true with ECGs by seeing that the maximal Lyapunov exponent is positive [2].
In nonlinear filter, though the fetal signal and maternal signal are similar in shape and spectral contents, we can separate the components by a very natural way: the magnitude and the heart beat. In fact, because the fetal heart is much smaller than the maternal heart, the fetal signal is much smaller than the maternal signal. Generally, the heart beat of the fetus is about one-third of the mother. In fact, our method is used for noise reduction, we just consider fetal signal as a contaminated noise.
Up to now, in order to deal with this problem, many works have been done and have given satisfactory results: “wavelet transform” was applied to extract wavelet-based features of fetal signal [3], “blind source separation (BSS)” was used to separate a set of source signal from numerous observed signals [4], “source extraction” is quite related to BSS except using additionally prior information about FECG [5], and some other method such as matched filtering [6], dynamic neural network [7], adaptive neurofuzzy inference systems [8], fuzzy logic [9], frequency tracking [10], and polynomial networks [11].
The technique we apply below is actually based on phase space reconstruction. The posttransient trajectory of the system is frequently confined to a set of points in state space called an “attractor” [12]. By using the delay coordinates, attractor is then empirically found to be constrained to a low-dimension manifold [13]. Hence, by estimating the attractor, noise can be reduced by projecting onto it. Whenever a multidimensional reconstruction of a signal can be approximated by a low-dimensional surface (or attractor), a projection onto this surface can improve the signal-to-noise ratio. In the present application, the fetal component is first treated as a contamination of the maternal ECG, whence noise reduction techniques are suitable for signal separation.
On the other hand, the extracted FECG recorded in the form of the nonstress test (NST) by using cardiotocography (CTG) was analyzed by wavelets to monitor fetal well-being [14].
2. Method
It has been proved that if a system is controled by an attractor, we can find out the dynamics of the full system just by single variables (theorem of Takens). In this paper, we use the delay coordinates of dimension. The geometry of a state space trajectory or a shape of attractor can be obtained by using delay coordinates to construct vectors valued time from a single-channel observation :
Here time is measured in the sampling intervals, is called delay or lag, and is the embedded dimension. All of the vectors are then inserted into a matrix, in this way we can draw delay plots: a plot pf versus a delay itself is able to reveal the characteristic of the attractor. According to Taken’s theorem, under general conditions, if the embedded dimension is large enough, the local topology of the attractor is preserved. This process is called a delay reconstruction in dimension (see Figure 1).
(a)
(b)
For noise reduction, according to Schreiber’s algorithm [12], we have to move all points of the phase space to the attractor manifold. First, we have to find the nearest points to , called neighbors, within the radius and the set of these points called . In addition, the number of neighbor should not be lower than a specified number which is usually 50 [15]. From here, we compute the mean:
Then we compute the covariance matrix:
is a weight matrix which is chosen to be diagonal with and large (about 1000) and all other diagonal entries . This makes the two largest eigenvalue of lying in subspace spanned by the first and the last coordinates of embedded space and prevents the correction vector from having any component in these direction.
Particularly, when using MATLAB, instead of using the for loop to compute covariance matrix, the method we use is a little bit different. First, we create a “deviation” matrix whose each column is the “deviation” vector . Then, the covariance matrix in equation is computed simply by . Empirically, this makes the result more accurate and much faster than using for loop.
In the next step, we determine the orthonormal eigenvector and the eigenvalue of . The eigenvector of represent the semiaxes of the ellipsoid best approximating the cloud of neighbor points . Ideally, the largest eigenvalues of the covariance matrix span the attractor manifold and the lower span the others. Projecting vector onto the subspace spanned by the largest eigenvalue will move it closer to the attractor manifold thereby creating a more accurate approximation of the true dynamics of the system, because the contaminating noise and the fetal signal span in another subspace. (where is the number of dimensions of the manifold that will be locally approximate by eigenvector corresponding to the largest eigenvalues).
When the projection is finished for all the points, we will get corrected vector because each element in scalar time series exists in vector. Therefore, we just average them all; this will not project the vector exactly to the manifold but will still move it closer to the manifold.
In this algorithm, there are three important parameters: the embedded window which is used to select components by time scale, the radius of cloud of neighborhood points which is used to select components by magnitude, and the number of dimension of the manifold .
In order to choose an optimal , many studies have been done. According to [2], should be large for two reasons. First, the larger is, the more deterministic signs are presented in the dataset. In fact, due to the fluctuating of body condition such as the respiration, the biological signal such as ECG becomes nonstationary and that violates the condition to apply the locally projective noise reduction. A solution is to increase . Instead of using in Taken’s theory, is used where is the number of nonstationary parameters [16]. Although, mathematically, must form a stationary chaotic system, this still works well if the nonstationary component varies slowly and rarely has any sudden change. The other reason is that with the large , the selectivity or appropriate neighbor is increased because we are using max norm to define the distance between neighbors. For more precise value of , the technique of Kennel and Abarbanel [17] is often used, which examines whether points that are near neighbors in this dimension are still near neighbor in the next dimension. Empirically, the larger is, the fewer false neighbors and the more fully the image unfolded.
In alliance with , choosing is a very important factor. A suitable should fulfill two criteria. First, must be large enough compared to the time scale or the two successive or they will be strongly correlated, that is, the value of the time series at must be significantly different form its value at . Therefore, we can gather enough information to successfully reconstruct all whole phase space with the reasonable choice of . Another criteria is that must not be too large than the time scale of the system, in which the two successive will be independent and uncorrelated [18, 19], and the system will lose memory of its initial state. According to [18], it has been showed that the optimal value of is typically around of the mean orbital period of the attractor. Often value of is around the correlation time.
Finally, should be larger than the noise amplitude and the fetal amplitude (when the fetal signal is considered as a contaminated noise), but still small enough not to average out the curvature radius of the time series to reserve the manifold shape. However, if we find the neighbors for all points in phase space by using fixed will not give a good result. In this paper, we use a graph between the number of neighbor versus distance to determine for each point in phase space (see Figure 2 for relation between ε and distance).
In the graph, each curve represent one point in phase space. Clearly, we can see that it is separated into 3 parts. Part (a) is those points which lay on the peaks of ECG, part (b) is those points which lay on the peaks of FECG, part (c) is other points. Thus, basing on the slope and the height of the curve, we can apply this simple rule for better filter: because our aim is to calculate the fetal heart beat, should be sufficiently large for the (b) part, and for avoiding distortion, should be fair small for the (a) and (c) parts. However, at the (b) part, there is, in fact, some points that does not lay on the FECG’s peak. The number of those points is few but the reason for this phenomenon is unknown. Fortunately, in some dataset, we find out that its slope is larger than the points lying on the FECG’s peak when the number of neighbor grow sufficiently large (about > 150). Note that this method is primarily used in the first filtering for we need to avoid the large distortion at this point (large error will severely affect the second filtering), while in the second filtering, we will come back with the fix to suppress the noise better. In further research, a better rule or algorithm for the neighbor finding procedure may solve this problem.
Besides, we approximate the attractor as a collection of locally linear manifolds. For example, the loop can be approximated by a collection of short line segments; in this case the approximating manifold is one dimensional. When the embedding dimension is larger than two, it can be appropriate to select locally linear manifold with dimension where . For instance , the manifold is locally planar.
In order to acquire FECG, in this paper, we will follow this strategy. First step, using Schreiber’s algorithm for the input data to extract to clear ECG, without noise and fetal signal. Second, subtracting the input data with the clear ECG to acquire the secondary input (including noise and fetal signal). The last step, using Schreiber’s algorithm again for the secondary input data to extract the clear FECG. The first and the last step should be repeated 2-3 times to get a better result. In sum, there are a lot of tradeoffs in choosing parameters and times of iteration, so that we have do a careful visual inspection of time series and few test runs to find the optimal results (see Figure 3).
3. Result
3.1. Artificial Signal
At first, we generate an artificial ECG by repeating one clear heart cycle of mother (the sample rate is 4 ms), plus the fake fetal signal created by scale artificial mother signal and then plus random noise following the Gaussian distribution (see Figure 4).
(a) Clean FECG RMS: Noise RMS = 1 : 1 (noise ratio = 1 : 1).
Following the following strategy, at first, to extract the clear ECG, we form an embedded window of 200 ms (equal 1/3 heart cycle of the mother) with and and using dynamic . around 2 is enough to reconstruct the phase space. The noise reduction ratio (see Figure 5).
After that, the secondary input (using dynamic ) is processed again with the embedded window . With of about 3, we will get the noise reduction ratio , and see that wave and wave are distorted. However, we just need to calculate the fetal heart beat so that the result is good enough (see Figure 6).
(a)
(b)
(b) Clean FECG RMS: Noise RMS = 1 : 5 (noise ratio = 1 : 5).
The first filtering (using dynamic ) is done with and , we will get the noise reduction equal to 2.58. However, the noise ratio is too large so that none of the characteristic of the FECG is reserved, thus the FECG extraction failed (see Figure 7).
(a)
(b)
(c)
3.2. In the Real World
The entire following sample is taken from MIT website [20].
(1) ECG 23rd Weeks
Because the fetus has grown a lot, the magnitude of the fetal signal is quite big so that its heart beat can be seen virtually.
The first filtering is quite good despite the distortion at the wave. For this result, we use , dynamic . The manifold dimension is 2 (see Figure 8).
For extracting the fetal signal, the next filtering is done with , using fix = 1.5. The manifold dimension is 15. The extracted FECG lacks many details of wave and wave, its shape is distorted a lot. Due to the nonstationary characteristic of the real ECGs, the neighbor versus distance graph becomes somehow unstable (some points do not lay on the FECG’s peak also have a slope as steep as the ones on the FECG’s peak), this makes it quite difficult to extract the FECG at the “nonstationary” section (using dynamic ) (see Figure 9).
For comparison, we add here a result of Martín-Clemente et al. [21] which uses fast ICA. Clearly, the applied sample is more stationary than ours and its amplitude is a little bit larger than either of them may be in the same weeks. Hence, the result is very clear while ours is still interfered by some unwanted harmonics as Figure 7. Though the main purpose is to measure only the fetal heart beat, our algorithm still meets the requirement with just one channel.
(2) ECG 25th Weeks
In this case, the fetal signal is very small; the noise ratio is large so that any effort to recover FECG fails because our algorithm will average both noise and FECG. This phenomenon happens because at this time, the fetus will create a membrane to cover itself, thus the fetal heart beat signal received at the probe will be greatly reduced (see Figure 10).
(a)
(b)
(3) ECG 38th Weeks
At this time, the fetus has grown a lot, its heart beat is bigger as well. Thus the noise ratio is lower, then we can calculate its heart beat once again.
The first filtering is done with , using dynamic . Then the second filtering is done with , using fix = 0.75 (see Figure 11).
(a)
(b)
4. Conclusion
According to all results we got, this algorithm can be mainly used for counting the heart beat because many details such as wave and wave are distorted, few peaks of FECG disappear sometimes. on the other hand, if we want a better result, the diagnosis will be done after the 38th week or before the 25th week. This is a physiological characteristic that was reconfirmed by our results.
However, these results still show the dominance over the linear method which is based on frequency domain to separate signal. In our method, we define the difference between ECG and FECG in a very natural way: amplitude and time scale. The natural language to implement such filtering procedures is in terms of the geometry in a reconstructed state space. Here, we have estimated the geometry of the maternal ECG using projections onto locally linear surfaces, which proves effective at filtering FECG from the measured signal.
In further works, some automatic information extraction method may be applied to test the efficiency of this method in reality as well as the development of a better algorithm or the research for the optimal parameters will be carried out to reserve as much information as possible for further diagnosis, manually or automatically.
Acknowledgments
This paper was partly supported research grants from Vietnam National University-Ho Chi Minh City, -Ho Chi Minh City International University, and Vietnam National Foundation for Science and Technology Development—NAFOSTED research Grant no. 106.99-2010.11.