PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Biomed Eng. Author manuscript; available in PMC 2010 November 10.
Published in final edited form as:
PMCID: PMC2978031
NIHMSID: NIHMS238349

Integrated Approach for Fetal QRS Detection

Abstract

Fetal magnetocardiography provides reliable signals of the fetal heart dynamics with high temporal resolution that can be used in a clinical setting. We present a robust Hilbert transform method for extraction of the fetal heart rate. Our method may be applied to signals derived from a single channel or an array of channels. In the case of multichannel data, the channels can be combined to improve signal-to-noise ratio for the extraction of fetal heart data. The method is inherently insensitive to fetal position or movement and, in addition, can be automated. We demonstrate that the determination of R-wave timing is relatively insensitive to waveform morphology. The method can also be applied if the data were preprocessed by independent component analysis (ICA). We compared the Hilbert method, ICA, ICA + Hilbert, and raw signals and found that the Hilbert method gave the best overall performance. We demonstrated that there were approximately 171 errors in 46 789 fetal heart beats.

Index Terms: Hilbert transform, independent component analysis (ICA), magnetocardiogram (MCG), quantron resonance system (QRS)

I. INTRODUCTION

The detection and analysis of fetal cardiac electromagnetic signals has many uses including fetal monitoring during labor [1], identification of fetal state [2], [3], monitoring of fetal arrhythmia [4], and automatic detection of fetal movement [5]. Kariniemi et al. reported the first fetal magnetocardiogram (fMCG) using a super conducting quantum interference device (SQUID) [6]. SQUID systems are noninvasive and have excellent signal-to-noise characteristics compared to fetal electrocardiogram (fECG) recordings taken from the mother’s abdomen, especially after the 24th week of gestation when the fECG is not easily detectable, and SQUID systems also have high temporal resolution (<1 ms) [7], [8]. Muli-sensor SQUID arrays have been used for fMCG detection [9]. Compared to a single-channel SQUID system, an array has the advantage of covering a large area simultaneously and offers the possibility of improved signal-to-noise ratio by combining the fMCG detected from several channels. One difficulty for the researcher interested in studying the FHR using SQUID technology is that the healthy fetus may move during the recording period that, in turn, may result in significant changes in the amplitude or morphology of the fMCG waveform [10]. Because a healthy fetus may become active, the algorithms developed to extract the timing of R-waves from the typical adult ECG may have high error rates when applied to either fMCG or fECG data. This paper presents a robust Hilbert-transform-based algorithm for improved measurement of fetal R-wave timing from a SQUID array.

FHR variability and fetal cardiac time interval analysis has several potential clinical applications [11]. Many indicators of FHR variability vary with gestational age reflecting the maturation of the fetus [12]–[14]. A long low-intensity speech sound has been shown to induce cardiac deceleration in the fetus that is dependent upon the fetal sleep state [15]. For FHR variability studies, the peak of each R-wave has to be detected in time and the R-to-R (RR) intervals are determined. The detailed examination of the fetal cardiac waveforms and beat-to-beat timing can only be accomplished if a high-quality signal is available.

Various groups have reported excellent quality fMCG recordings using single- and multichannel SQUID systems [16], [17] in which the fetal P-wave, quantron resonance system (QRS) complex, and T-wave can be determined even without signal averaging. The fMCG can easily reveal arrhythmias and ectopic beats. By means of independent component analysis (ICA), the fMCG of twins have been separated [18]. In the later stages of gestation, when the fetal head is usually engaged into the birth canal and physically constrained, the fetal trunk may still flex or rotate during breathing movements or during gross body movements. Thus, the fMCG waveforms may show considerable variation in polarity, R-wave morphology, and signal strength across the maternal abdomen.

Several single-channel algorithms have been developed to extract RR intervals from adult ECG [19], [20]. Many algorithms follow a two step process; first, the ECG signal is passed through a digital filter or postprocessing step designed to enhance the QRS complex with respect to the P- and T-waves or simply to remove noise. In the second step, a template detection method or a threshold detection scheme is used to identify the R-wave positions and the fiducial point is usually taken as the highest point on the R-wave. Peak detection may be further refined to adjust for amplitude variation and to search for missing R-waves. If the morphology (mono-polar or bipolar nature) of the waveform remains homogenous throughout the entire duration of the recording time, the aforementioned methods can be reliably applied to find the fiducial points. The main objective of this paper is to present a new postprocessing method based on the Hilbert transform that allows the combination of multichannel SQUID data to enhance signal-to-noise ratio and better delineate the fMCG QRS complex, thus simplifying R-wave location.

The Hilbert transform is widely used in signal processing and has been applied to the problem of R-wave extraction [21], [22]. In the method presented by Bolton and Westphal [22], the cardiogram is first filtered by the method of first differences (numerical differentiation), and then, the Hilbert transform is performed on the differentiated signal. Thus, the two-step operation yields a high-pass filtered cardiogram without phase distortion. Since the QRS complex is of higher frequency than the P- and T- waves, the R-waves are significantly enhanced and the baseline is removed. Commonly, the peaks of the Hilbert transformed cardiogram are identified and used as guides to find the peaks in the unprocessed data. To deal with the dynamical changes in the polarity of the fMCG waveform and fetal repositioning, we propose a modified high throughput (HT) approach extending the Hilbert transform methods developed for ECG analysis. Nygards and Sornmo [23] have shown that the magnitude of the Hilbert analytic function, which is a complex signal formed by combining the signal (real component) with its Hilbert transform (imaginary component). The magnitude of the analytic function is always positive, and thus, independent of signal polarity, and to a large extent, is independent of QRS waveform morphology. These two properties allow the magnitude of the analytic function from multiple sensors to be combined linearly: 1) to improve the signal-to-noise ratio; 2) to eliminate transient signals that present as false R-waves; and 3) to maintain consistent fMCG signal amplitude during the entire recording.

In parallel to HT, ICA is used for extraction of fetal QRS complexes, especially in noisy recordings [24]–[26]. ICA is capable of separating multichannel recordings into the various sources creating the observed signals and has recently been used to extract an estimate of fetal heart rate in real time [27]. The basic requirement for ICA to perform properly is statistical independence. This basic requirement is fulfilled in fMCG recordings because maternal and fetal cardiac signals are generated by different sources that for practical purposes act independently of each other. We explore the idea that the QRS complex may be further improved by applying HT methods to the fMCG components found by ICA.

II. METHODOLOGY

A. Data Collection

Data were collected using the SQUID array for reproductive assessment (SARA). In the SARA instrument, the patient interface is a 151 channel array of SQUID gradiometers with diameter of 2 cm and baseline of 8 cm. [16]. SARA has a concave surface with sensors spaced approximately 3 cm apart over an area of 850 cm2 spanning the maternal abdomen longitudinally from the symphysis pubis to the uterine fundus and a similar distance laterally. The pregnant woman can comfortably lean forward with her abdomen pressing gently on the concave surface thereby positioning the fetus in proximity to several of the 151 SQUID sensors regardless of the fetal position within the womb.

Fifty five recordings of 6-min duration were collected from fetuses ranging in gestational age from 28 to 39 weeks using the SARA instrument. Each SQUID channel was sampled at a rate of 312.5 Hz using a bandwidth set for 0–100 Hz. The data were then processed to remove the baseline and high-pass filtered using a fourth-order high-pass Butterworth digital filter with 1 Hz-3 db frequency. The filter was applied in both the forward and reverse directions yielding an eighth-order high-pass response with zero group delay (symmetric noncausal impulse response).

B. HT Approach

1) Final Submission: For a real-time function x(t) the Hilbert transform is defined as

h(t)=1πP.V.x(τ)1tτdτ
(1)

where P.V. denotes Cauchy principal value. For a channel of interest, we compute the complex helical sequence, rn, such that the nth (t = n × sampling frequency) complex point is given by

rn=xn+ihn.
(2)

In practice, the Hilbert() function in the Matlab (The Mathworks, Inc., Natick, MA) was used to compute the sequence rn. Next, the distance Rn between successive points in the series rn that we define as rate of change of Hilbert amplitude (RHA) is computed by

Rn=(xn+1xn)2+(hn+1hn)2.
(3)

The series Rn contains R-wave-like features that are easy to identify by simple methods. We must mention that Rn is equivalent to the envelope function discussed in [20] that was used as a guide to peak selection. The method presented here is computationally efficient and is designed to emphasize the characteristics of the R-wave, thus making them much easier to distinguish from other cardiac components and the base line. There are three important properties that make this method especially useful in fMCG analysis.

  1. Rn is always positive, so that the polarity of the fMCG is not important.
  2. A typical bipolar R-wave has a unipolar response with only one peak.
  3. For a (stationary) sinusoidal signal, the response is given by Rn=|rn|22cosπυ, where υ is the normalized frequency range 0–1.

Property 1 is useful because signals from multiple channels or components can be integrated (summed) without regard to polarity, thus improving the signal-to-noise ratio of the QRS complex over a single channel. In this regard, the algorithm is very robust in that the signals from many SQUID channels may be combined and the fetal R-wave is still the dominant signal even though the fetal signal is found in only a few channels. Property 2 simplifies the identification of R-waves that present as a bipolar signal. This eliminates some problems in fetal MCG analysis where the waveform may change polarity due to fetal movement, or only bipolar R-waves are recorded. Property 3 gives insight into the frequency-domain characteristics of the method. Although (3) indicates that the method is nonlinear, nonetheless, this transformation behaves much like a high-pass filter in that the lower frequency P- and T-waves are suppressed, and in the low frequency range, the magnitude of a sinusoidal signal is approximated by Rn ≈ |rn | πυ.

To graphically demonstrate the versatility of the RHA algorithm, positive, negative, and bipolar prototype R-waves were computed and presented in Fig. 1 along with the output of the RHA algorithm and the output of the Hilbert transform algorithm described by Bolton and Westphal [22]. This Hilbert method acts as a kind of high-pass filter that does not time shift the peak, whereas the RHA is based on the Hilbert amplitude that is always positive (see Fig. 1). For the prototype R-waves, the RHA always produced a monopolar response with a single peak.

Fig. 1
Three prototype R-waves and corresponding results after being processed by the RHA algorithm and the method of Bolton and Westphal. (A) Gaussian pulses representing monopolar R-waves. (D) Derivative of the Gausian shown in (A). (G) Combination of positive ...

We advocate summing the RHA from multiple channels in order to gain better signal quality, but for this method to be effective, it is desirable that the maximum from the multiple channels be coincident, regardless of waveform morphology. In Fig. 1(B), (E), and (H), the fiducial points are simply located at the peaks. For the Gaussian example (Fig. 1), the fiducial points found by the RHA lie at the same point in time for all waveforms. However, other analytic functions and actual recorded waveforms may not be so well behaved. Given an arbitrary continuous function of time f(t) and its Hilbert transform h(t), we now define the amplitude square function as

M(t)=(df(t)dt)2+(dt(t)dt)2.
(4)

Note that the RHA is proportional to the square root of M(t). The actual maxima of an arbitrary M(t) can easily be found by setting the time derivative of M(t) equal to zero resulting in

ff+hh=0.
(5)

Since the Hilbert transform of a derivative is equal to the derivative of the Hilbert transform, the peak of the RHA is defined totally by the velocity and acceleration of the signal f(t). In general, if f(t) is either purely odd or purely even, then the RHA will be maximum at t = 0. However, if f(t) is neither odd nor even, the situation is not as simple to analyze. Using the even function f(t) = 1/(1 + t2) as a model of the R-wave, we define Y(t) = αf + (1−α)f′ where α is allowed to vary from 0 to 1. Thus, Y(t) is a linear mix of an odd function and an even function as determined by the parameter α. We numerically studied the position of the maxima of the RHA of Y(t) compared to the maxima of Y(t) as a function of α. The results are shown in Fig. 2. When α = 0, then Y(t) is odd and is a bipolar signal and the maximum occurs well before t = 0. This time is used to normalize the two curves presented in Fig. 2. As expected from the arguments made earlier, the maximum of the RHA of Y(t) occurs at t = 0 when α = 0 (odd case) and α = 1 (even case). However, the position of the maximum of the RHA(Y) is nonzero for all other values of α. We use the term dispersion to represent the spread in time induced by the RHA as a function of waveform morphology. Though not zero, the dispersion of the RHA of the model function is significantly less than the dispersion inherent in the signal. The only practical way to understand the dispersion from the fMCG is to compare the RHA from several MCG channels simultaneously.

Fig. 2
Normalized dispersion of Y(t) and RHA(Y) as a function of the parameter alpha such that X(t) = 1/(1 + t2) and Y(t) = αX(t) + (1−α)X′(t). The dispersion is normalized by the fiducial point set by the peak of Y(t) when α ...

C. ICA

Assuming an n-dimensional vector x of observables generated by a mixture of n independent sources s, we can formulate the system as

x=As.
(6)

The mixture matrix containing the elements aij will be denoted by A or more realistically

x=As+N
(7)

to account for additive noise N from unknown sources. ICA uses the observed signals denoted by x to estimate A and s with the assumptions of independence and a non-Gaussian distribution [28]. Once the matrix A has been estimated, the inverse W can be computed to obtain the independent components

s=Wx.
(8)

The particular ICA algorithm used for separation of the fMCG components is the FastICA, developed by Hyvarinen and Oja [29], [30]. FastICA finds the columns of the separating matrix by maximizing the absolute value of kurtosis, which works particularly well for QRS extraction due to high value of kurtosis for QRS complexes.

D. Example Waveforms

Fig. 3(A)–(E) shows the varying morphology of single-channel fetal cardiac activities of five different randomly chosen fetuses and Fig. 3(F)–(J) is the corresponding RHA. In Fig. 3(A), the R-wave is monopolar with a clear P-wave, no obvious T-wave, and some low-frequency baseline noise. The corresponding single-channel RHA is shown in Fig. 3(F). Fig. 3(B) displays beat-to-beat change in waveform morphology changing from a bipolar presentation to a monopolar presentation in only two fetal cardiac cycles. Note that the RHA shown in Fig. 3(G) is consistent both in amplitude and morphology despite the changes observed in Fig. 3(B). Fig. 3(C) is a purely bipolar presentation. In such cases, the dominance of the positive or negative peak may change quickly. The corresponding RHA, Fig. 3(H), distinguishes the R-wave position clearly from the baseline. Fig. 3(D) is dominantly negative polarity and shows considerable baseline instability while the resulting RHA is relatively clean and positive. Fig. 3(E) is a negative monopolar R-wave with variable amplitude. Note that the amplitude variation is captured in Fig. 3(J). Fig. 3(K)–(O) shows the ICA components corresponding to the fMCG signals shown in Fig. 3(A)–(E). All the ICA components show a constant polarity including the ICA component found for the purely bipolar representation seen in Fig. 3(C). It should be noted that the RHA presented in Fig. 3(F)–(J) is for a single channel and is not the summed waveform ultimately used to extract R-waves. The RHA may be applied to a single ICA component, or if more than one R-wave like ICA component is found, then the RHA can be applied to combine the ICA components.

Fig. 3
One-second duration of fMCG of five fetuses of different gestation ages and their QRS complexes identified by RHA compared to the ICA component. (A) Positive monopolar waveform. (B) Combination of a bipolar and positive waveforms. (C) Bipolar waveform. ...

E. Postprocessing

In this paper, the summed RHA method, FastICA, RHA applied to two or more FastICA components, and simple filtering are considered postprocessing steps designed to simplify the extraction of the fetal R-waves. Since the maternal cardiac signal is dominant over the fetal signal, the maternal MCG was first removed from all SQUID channels by a signal space projection algorithm [31] before applying any of the methods tested. To further reduce the influence of the residual maternal MCG, only data from the lower 82 channels of the SARA system, far from the maternal heart but nearer to the fetal heart, were used in our analysis. To be clear, the FastICA was not used to separate the fetal and maternal MCGs in this study, but rather it was only used to separate the fMCG from background signal. After removal of the maternal MCG, each dataset was postprocessed in four ways: 1) compute the RHA of each channel and sum (Hilbert); 2) apply the FastICA and select the dominant fetal component (ICA); 3) compute the RHA of all FastICA components containing fMCG and sum (ICA+Hilbert); and 4) manually select the channel with the largest fMCG signal (filtered). The label “filtered” refers to the 1–100 Hz bandwidth imposed by the initial filtering operations. Ideally, after applying any of the earlier postprocessing methods, the fiducial points may be identified by simply extracting all of the maxima and testing which maxima fall above a selected threshold.

To illustrate this process, the results from a single recording are presented in Fig. 4. Here, the local maxima from each method were extracted and normalized by the largest positive value, and then, used to build the histograms shown in Fig. 4(A)–(D). For clarity, the histograms were scaled by taking the log(n+1) of each histogram bin. Fig. 4(E)–(H) shows the resulting RR intervals obtained after selecting a threshold between 0 and 1 and then extracting the remaining maxima. For this dataset, the ICA and Hilbert method (RHA) separate the maxima into two nonoverlapping distributions, as shown in Fig. 4(B) and (D), where the histograms indicate a completely resolved distribution of R-wave maxima separated from the distribution of lower amplitude maxima arising from noise and/or other fMCG components. Fig. 4(F) and (H) are the graphs of the corresponding RR intervals for any value of threshold falling between the two resolved distributions. In contrast, Fig. 4(A) and (C) shows overlapping distributions. In these two cases, a single threshold cannot be selected that cleanly separates the maxima of R-waves from maxima originating from other sources. As an example, Fig. 4(E) shows an upward spike in the graph of RR intervals near the beginning of the dataset indicating that at least one true R-wave maxima fell below the threshold. We refer to this case as a “missed” beat. Fig. 4(G) shows several downward spikes indicating that a few maxima were falsely classified as R-waves. We refer to this case as “extra” beat.

Fig. 4
Histograms of peaks found in data and resulting RR intervals. (A)–(D) Histograms from the four methods described in the text. Note the bimodal nature of the histograms, with the distribution of high-amplitude R-waves lying to the right in each ...

F. Performance Measures

If the distribution of the local maxima, as shown in Fig. 4(B) and (D), is clearly separated, then the selection of the optimum threshold is trivial. In many cases, the optimum choice is not clear as seen in Fig. 4(C). In this case, there will always be errors regardless of the threshold selection. There are two types of errors to be considered: first, a “missed” beat, as seen in Fig. 4(E), results in an upward spike when the maximum of a true R-wave falls below the selected threshold. Second, an “extra” beat, as seen in Fig. 4(G), results in a downward spike when the maximum of a spurious signal exceed the selected threshold. Since Fig. 4(F) and (H) do not show either type of spike, it is clear that all the spikes shown in Fig. 4(E) and (G) are due to misclassified maxima. Further, ectopic beats or arrhythmias should be seen equally in the results of all four postprocessing methods. Two algorithms were designed to estimate the number of missed and extra beats so that the four postprocessing methods could be objectively compared. The estimators were used to automatically select an optimal threshold by finding the threshold that yields the minimum total error.

Examination of our data revealed that the fetal heart rate fell within the range of 105–190 beats per minute (0.5714 and 0.3158 s). These limits were obtained by constructing a histogram of RR intervals found by using combined results from all four methods after selecting a discrimination threshold for each dataset and method that minimized the sum of missed and extra beats. The resulting distribution was nearly Gaussian and was centered at 145 beats per minute (bpm). The range of 105–190 bpm includes the tails of the Gaussian, and therefore, represented the range of normal heartbeats in our data. For purposes of analysis, a “missed” beat is defined by an RR interval longer than 0.5714 s and an “extra” beat is defined by an RR interval shorter than 0.3158 s. A set of RR intervals is obtained by first normalizing the amplitude of the postprocessed dataset, specifying a threshold in the range 0–1.0, and then, extracting the time value of all maxima that exceed the threshold. Once a set of RR intervals are obtained, they may be analyzed to estimate the number of missing and extra beats by using the estimators presented later. Both estimators extract and use the median RR interval as estimate of the expected duration of a single heart beat. However, if the threshold is a poor choice, the median RR interval may fall outside the range of the normal heart beat, in which case the median RR is assigned a value of 0.4068 that is midway between the upper and lower limits. For any RR interval that exceeds 0.5714 s, the number of missed beats may be estimated by dividing that RR interval by the median RR interval and then subtracting one. Applying this process to all RR intervals that exceed the upper limit leads to the formula

missed=RRiRRmM
(9)

where missed is the estimated number of missed beats, RRi is an RR interval that exceed 0.5714, RRm is the median RR interval, and M is the total number of intervals that exceed 0.5714. The missed estimate is then rounded to the nearest integer. In a similar fashion, the number of extra beats is estimated by totaling the number of short duration RR intervals and subtracting an estimate of the number of intervals associated with a true R-wave. The extra beat estimator is then expressed as

extra=ERReRRm
(10)

where extra is the estimated number of extra beats, RRe is an RR interval less than 0.3158, RRm is the median RR interval, and E is the total number of intervals less than 0.3158. The extra estimate is then rounded to the nearest integer.

By summing the estimates of the extra and missed beats, the total number of errors is obtained for each threshold. By definition, the best threshold is the one that yields the lowest total errors. If the error estimate is zero, then all RR intervals fall within the upper and lower bounds set for the estimators and one is assured that the threshold cleanly separates the maxima of the RR intervals from all other maxima.

III. RESULTS AND DISCUSSION

Fifty five datasets were processed using the four methods described earlier and the RR intervals were computed with thresholds ranging from 10% to 80% of maximum in 2% steps. The number of missed and extra beats was computed for each threshold and the threshold that yielded the lowest total error was selected as the “best” threshold. The values from the two estimators associated with the “best” threshold were recorded for analysis.

A. Validation of Estimators

Since the missed and extra beat estimators are used to judge the postprocessing methods examined in this paper, it is essential to demonstrate their accuracy. As stated, if the error estimate is zero, then one is assured that the two distributions of the local maxima have been cleanly separated. All results were scanned for cases where at least one of the methods yielded zero for both estimators. In these instances, the number of RR intervals is the exact number of heart beats for that fMCG recording. Using only fMCG recordings where the exact number of heart beats was known, the exact count was subtracted from the number of observed RR intervals, yielding the net error count for the other postprocessing methods. The observed values of the net error count ranged from −228 to 260. To test the performance of our (“missed” and “extra” beat) estimators, the net error count was estimated by subtracting the “missed” beat estimate from the “extra” beat estimate both of which are obtained using (9) and (10), respectively. There was a good degree of correlation (R = 0.999) between the observed net error and the estimated net error indicating the robustness of our estimators.

B. Dispersion

Fig. 5(A)–(D) presents various example waveforms that may be seen simultaneously with a SQUID array; Fig. 5(A) is the classical monopolar R-wave, Fig. 5(B) is bipolar, Fig. 5(C) appears to be a mix of the first two, and Fig. 5(D) is the worst case example with respect to dispersion. Each waveform is normalized with respect to the maximum value of Fig. 5(A) and time zero is defined by the peak value of the summed RHA (Hilbert) for the arbitrarily chosen R-wave. Fig. 5(E)–(H) are histograms of each channel’s dispersion, respectively, computed by taking the time difference between the fiducial points found using the Hilbert method and the peaks of the RHA applied to the selected channels. The time is expressed in sample intervals. Fig. 5(E)–(G) shows only plus or minus one count difference from the Hilbert method indicating that the timing of the three waveforms shown in Fig. 5(A)–(C) as defined by their RHAs is consistent with respect to the summed RHAs (Hilbert) and that the dispersion is small. Fig. 5(H) is the worst case dispersion for this dataset. The dispersion histogram for each channel was computed and the mean dispersion for each channel plotted in Fig. 5(I). The majority of all RHA peaks fall within plus or minus one count of the peak found by the Hilbert method. The channels numbers in the range of 25–35 show considerable dispersion, but they are of generally low amplitude, and therefore, contribute little to the determination of the position of the peak found by summing the RHAs. Channel 57 is the worst case and is the waveform shown in Fig. 5(D). Future implementations of the Hilbert method could test for dispersion and discard problematic channels.

Fig. 5
(A)–(H) Time is expressed as sample intervals and t = 0 is taken as the summed RHA fiducial point for an arbitrarily chosen fetal heart beat. (A)–(D) Selected MEG channels showing the variation of fetal R-waves. All amplitudes were normalized ...

C. Comparison of Methods

The “missed” beat and “extra” beat estimators make it easy to tabulate the performance of the four postprocessing methods studied in this paper. The summed results for 55 datasets are presented in Table I. The first three rows summarize the results from the “extra” and “missed” beat estimators. In the fourth row, the actual number of beats was estimated for each method and the totals differed only in the third decimal place. Researchers have defined efficiency as the ratio of the total number of beats minus errors to the total number of beats. Thus defined, the efficiency is presented in the fifth row expressed as a percentage. The sixth row is labeled total best case. Here, the error estimate of each method was compared with the other three methods. If the error was equal to or less than all other methods, the example was classified as “best.” The value indicated is the total number of “best” classifications. For the next three rows, the value of the “extra,” “missed,” and “error” estimates were examined and marked if found to be zero. The total number of marked zeros is listed in the table for each category.

TABLE I
Summary of various measures of performance (see text for details)

The Hilbert method scored the best in all categories. Surprisingly, the ICA results were mixed showing a very high number of errors and yet ranking second in categories defined by zero error, zero extra beats, and best case. It was observed that the ICA either yielded excellent results or rather poor results. To quantify this, the standard deviation of the error, extra, and missed values were computed, excluding the cases where a method gave a zero value for the particular estimate of performance. Referring to the last three rows of Table I, the standard deviations for the ICA are highest. One interpretation of this result is that the ICA does not present consistent amplitude over the time course of each study and that limiting the RR interval extraction to a single threshold causes significant error even though the signal-to-noise ratio of the signal may be rather good. This idea is partially supported by noting that the number of ICA missed beats is rather high indicating that the signal frequently dropped below optimum threshold, but only in a limited number of recordings. Since the ICA is known to extract signals from noise, it seems unlikely that noise explains the behavior seen in this study. Rather, one might consider that the property (selective gain) that allows the ICA to select out the signal might also make it occasionally sensitive to fetal or maternal movement. An example that clarifies this idea is presented in Fig. 6 that shows simultaneous waveforms produced by the four postprocessing methods: (A) Hilbert, (B) ICA, (C) ICA + Hilbert, and (D) filtered. In Fig. 6(B), the ICA component briefly reverses polarity while the other three methods consistently yielded monopolar R-waves. For the ICA component shown, a single threshold can not be found that will properly identify the R-wave peaks that results in a high number of missed beats for the ICA. Examination of the heights of the R-waves in Fig. 6(A) shows that the Hilbert method responded to an overall increase in summed signal during the brief polarity reversal of the ICA component, indicating that one or more channels momentarily presented increased signal amplitude changes that were problematic for the ICA. The combination of ICA and Hilbert improved the performance of the ICA by reducing the number of errors, especially the missed beats. Fig. 6(C) demonstrates this point. However, the combination of ICA and Hilbert produced the largest number of extra beats. The FastICA tends to select for signals that have high kurtosis (R-waves or spikes) and the RHA acts similar to a high-pass filter, thus the ICA+Hilbert combination may exaggerate transients giving rise to the observed increase in extra beats compared to the other methods.

Fig. 6
Plots of 8-s duration showing the output from the four postprocessing methods. (A) Hilbert. (B) ICA. (C) Hilbert + ICA. (D) Filtered. In (B), the ICA component shows a momentary polarity inversion that results in “missed” beats. In (C), ...

It should be noted that the 171 errors found by the Hilbert method may not have been errors at all, but rather properly captured tachycardic or bradycardic beats. In a given dataset, if the majority of fetal RR intervals fall within the 105–190 bpm range, then threshold selection based on minimizing the total estimated error count finds the optimal R-wave amplitude threshold. We anticipate that the occasional irregular beat will have approximately the same amplitude R-wave as any other beat, and if so, will be properly classified. For now, RR intervals that fall outside the 105–190 bpm window must be manually verified. Finally, if tachycardic or bradycardic beats are in the data, then all four methods should indicate these intervals as either missed or extra. Therefore, there should be no bias for one method over another in the comparisons presented in Table I.

D. Automation

Except for the removal of the maternal heart signal using the orthogonal projection algorithm, the Hilbert method demonstrated here could be completely automated. The summed RHA, RR interval extraction, error estimation, and threshold selection are all tested algorithms that do not require user input. If implemented as a single automated approach, the result would be interpreted as 35 of 55 recordings (63.6%) were processed free of error with reliable confirmation. The remaining 20 recordings had only 171 possible errors (out of approximately 46 600 total beats) that could be presented to an operator for validation or editing. Considering the results presented in Fig. 5, it may be possible to improve the Hilbert algorithm by automatically excluding problematic channels based upon poor dispersion or low signal-to-noise ratio.

IV. CONCLUSION

FastICA and the rate of change of the Hilbert transform amplitude, as a postprocessing step for RR interval extraction, have been studied in this paper. The principle method for comparison was the use of a single threshold discriminator and “missed” beat and “extra” beat estimators. Our premise is that if the distribution of the maxima of the R-waves is completely distinct from the distribution of all other maxima, then the method clearly shows good signal-to-noise characteristics. Using this technique, our analysis indicates that the Hilbert transform method is particularly useful to combine SQUID array data for the detection of R-waves and consistently outperformed the other three methods examined. If the ICA is used to extract the fMCG, then the Hilbert method may be applied to one or more ICA components as needed to improve the signal quality. Additional benefits of the Hilbert method are that it is computationally efficient and easy to implement and requires little user intervention. For the case where the fetal heart rate falls within the range of 105–190 bpm, we demonstrate that the Hilbert method can be used to automate the extraction of the fetal heart rate, achieving 99.6% efficiency.

Acknowledgments

This work was supported by the National Institutes of Health under Grant 5R01-NS-36277-05A1 and Grant 5R33-EB-00978-02.

Footnotes

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Contributor Information

James D. Wilson, Graduate Institute of Technology, University of Arkansas, Little Rock, AR 72204 USA.

R. B. Govindan, Department of Obstetrics and Gynecology, University of Arkansas for Medical Sciences, Little Rock, AR 72205 USA.

Jeff O. Hatton, Department of Obstetrics and Gynecology, University of Arkansas for Medical Sciences, Little Rock, AR 72205 USA.

Curtis L. Lowery, Department of Obstetrics and Gynecology, University of Arkansas for Medical Sciences, Little Rock, AR 72205 USA.

Hubert Preissl, Magnetoencephalography (MEG) Center, University of Tübingen, Tübingen 72076, Germany, and also with the Department of Obstetrics and Gynecology, University of Arkansas for Medical Sciences, Little Rock, AR 72205 USA.

REFERENCES

1. Thaler I, Timor IE, Goldberg I. Interpretation of the fetal ECG during labor: The effect of uterine contractions. J. Perinat. Med. 1988;vol. 16:373–379. [PubMed]
2. Davidson SR, Rankin JH, Martin CB, Reid DL. Fetal heart rate variability and behavioral state: Analysis by power spectrum. Amer. J. Obstet. Gynecol. 1992 Sep.vol. 167(no. 3):717–722. [PubMed]
3. Nijhuis JG. Fetal Behavior: Developmental and Perinatal Aspects. Oxford, U.K.: Oxford Univ. Press; 1992. Ed.
4. Wakai RT, Strasburger JF, Li Z, Deal BJ, Gotteiner NL. Magnetocardiographic rhythm patterns at initiation and termination of fetal supraventricaular tachycardia. Circulation. 2003;vol. 107:307–312. [PubMed]
5. Wakai RT, Wang M, Leuthold AC, Martin CB. Foetal magnetocardiogram amplitude oscillations associated with respiratory sinus arrhythmia. Physiol. Meas. 1995;vol. 16:49–54. [PubMed]
6. Kariniemi V, Ahopelto J, Karp P, Katila TE. The fetal magnetocardiogram. J. Perinat. Med. 1974;vol. 2:214–216. [PubMed]
7. Carter MC, Gunn P, Beard RW. Fetal heart rate monitoring using the abdominal fetal electrocardiogram. Br. J. Obstet. Gynaecol. 1980 May;vol. 87(no. 5):396–401. [PubMed]
8. Jafari MG, Chambers JA. Fetal electrocardiogram extraction by sequential source separation in the wavelet domain. IEEE Trans. Biomed. Eng. 2005 Mar.vol. 52(no. 3):390–400. [PubMed]
9. Kandori A, Tsukada K, Horigome H, Asaka M, Shiqemitsu S, Takahashi M, Terada Y, Kubo T, Matsui A, Mitsui T. Multichannel SQUID system for detecting fetal magnetocardiograms. presented at the 2nd Int. Conf. Bioelectromagn.; Melbourne, Australia. 1998. Feb.,
10. Zhao H, Wakai RT. Simultaneity of foetal heart rate acceleration and foetal trunk movement determined by foetal magnetocardiogram actocardiography. Phys. Med. Biol. 2002;vol. 47:839–846. [PubMed]
11. Leeuwen PV. Fetal magnetocardiography: Time intervals and heart rate variablility. Neurol. Clin. Neurophysiol. 2004 Nov.vol. 46:1–7. [PubMed]
12. Leeuwen PV, Lange S, Bettermann H, Gronemeyer D, Hatzmann W. Fetal heart rate variability and complexity in the course of pregnancy. Early Human Dev. 1999 Apr.vol. 54(no. 3):259–269. [PubMed]
13. Echeverria JC, Woolfson MS, Crowe JA, Hayes-Gill BR, Pieri JF, Spencer CJ, James DK. Does fractality in heart rate variability indicate the development of fetal neural processes? Phys. Lett. A. 2004;vol. 331:225–230.
14. Govindan RB, Wilson JD, Preissl H, Eswaran H, Campbell JQ, Lowery CL. Detrended fluctuation analysis of short datasets: An application to fetal cardiac data. Phys. D. 2007;vol. 226:23–31.
15. Groome LJ, Mooney DM, Holland SB, Smith LA, Atterbury JL, Dykman RA. Behavior state affects heart rate response to low-intensity sound in human fetuses. Early Human Dev. 1999;vol. 54:39–54. [PubMed]
16. Lowery L, Campbell JQ, Wilson JD, Murphy P, Preissl H, Malak S, Eswaran H. Noninvasive antepartum recording of fetal S-T segment with a newly developed 151-channel magnetic sensor system. Amer. J. Obstet. Gynecol. 2003;vol. 188:1491–1497. [PubMed]
17. Stinstra J, Golbach E, Leeuwen PV, Lange S, Menendez T, Moshage W, Schleussner E, Kaehler C, Horiqome H, Shiqemitsu S, Peters MJ. Multicentre study of fetal cardiac time intervals using magnetocardiography. BJOG: Int. J. Obstet. Gynecol. 2002;vol. 109(no. 11):1235–1243. [PubMed]
18. Burghoff M, Leeuwen PV. Separation of fetal and maternal magnetocardiographic signals in twin pregnancy using independent component analysis (ICA) Neurol. Clin. Neurophysiol. 2004 Nov.vol. 39:1–4. [PubMed]
19. Hamilton PS, Tompkins WJ. Quantitative investigation of QRS detection rules using the MIT/BIH arrhythmia database. IEEE Trans. Biomed. Eng. 1986 Dec.vol. BME-22(no. 12):1157–1165. [PubMed]
20. Benitez S, Gaydecki PA, Zaidi A, Fitzpatrick AP. The use of the Hilbert transform in ECG signal analysis. Comput. Biol. Med. 2001;vol. 31:399–406. [PubMed]
21. Benitez S, Gaydecki PA, Zaidi A, Fitzpatrick AP. A new QRS detection algorithm based on the Hilbert transform. Comput. Cardiol. 2000;vol. 27:379–382.
22. Bolton RJ, Westphal LC. ECG display and QRS detection using the Hilbert transform. Comput. Cardiol. 1985:463–466.
23. Nygards ME, Sornmo L. Delineation of the QRS complex using the envelope of the e.c.g. Med. Biol. Eng. Comput. 1983;vol. 21:538–547. [PubMed]
24. Campbell JQ, Eswaran H, Wilson JD, Murphy P, Lowery CL, Preissl H. Fetal magnetocardiographic source separation: Independent component analysis techniques and signal-space projection. Int. J. Bioelectromagn. 2005;vol. 7:329–333.
25. Mantini D, Comani S, Pennesi P, Lagatta A, Cancellieri G. Fetal tailoring of the independent component analysis to multi-channel fMCG recordings for an optimal reconstruction of the fetal cardiac signal. Biomed. Tech. 2004;vol. 48:186–188.
26. Comani S, Mantini D, Pennesi P, Lagatta A, Cancellieri G. Independent component analysis: Fetal signal reconstruction from magnetocardiographic recordings. Comput. Methods Programs Biomed. 2004;vol. 75:163–177. [PubMed]
27. Waldert S, Bensch M, Bogdan M, Rosenstiel W, Scholkopf B, Lowery CL, Eswaran H, Preissl H. Real-time fetal heart monitoring in biomagnetic measurements using adaptive real-time ICA. IEEE Trans. Biomed. Eng. 2007 Oct.vol. 54(no. 10):1867–1874. [PubMed]
28. Comon P. Independent component analysis—A new concept? Signal Process. 1994;vol. 36:287–314.
29. Hyvarinen A, Oja E. A fast fixed point algorithm for independent component analysis. Neural. Comput. 1997;vol. 9:283–292.
30. Hyvarinen A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans. Neural. Netw. 1999 May;vol. 10(no. 3):626–634. [PubMed]
31. Vrba J, Robinson SE, McCubbin J, Lowery CL, Eswaran H, Wilson JD, Murphy P, Preissl H. Fetal MEG redistribution by projection operators. IEEE Trans. Biomed. Eng. 2004 Jul.vol. 51(no. 7):1204–1218. [PubMed]