Abstract:
This invention provides a novel technique separating two mixed signals. The method is particularly useful for fetal electrocardiogram estimation. A maternal electrocardiogram signal recorded at the thorax area is nonlinearly mapped onto an electrocardiogram signal recorded at the abdomen using polynomial networks. The fetal electrocardiogram component is then estimated by subtracting the nonlinearly mapped maternal electrocardiogram signal from the abdominal electrocardiogram signal.

Description:
BACKGROUND 
   a. Field of the Invention 
   This invention relates to a method of separating a signal from a mixture of two signals where a distorted version of one of the signal is available. It is particularly useful for the estimation of a fetal electrocardiogram (FECG) using signals recorded at the surface of the body. The FECG signal reflects the electrical activity of the fetal heart. It contains information on the health status of the fetus and therefore, an early diagnosis of any cardiac defects before delivery increases the effectiveness of the appropriate treatment. There are several technical problems associated with the extraction of FECG from signals recorded at the abdominal surface which result in a poor signal quality. The main sources of possible noise include the maternal ECG, the maternal electromyogram EMG, (50 or 60) Hz power line interference, baseline wandering and random electronic noise. To simplify the problem, it is assumed that state of the art low noise electronic amplifiers with high common mode rejection ratio are used to eliminate the power line interference and electronic random noise. The EMG noise can also be reduced, but not necessarily eliminated, with the use of classical low pass filtering techniques. Therefore, this invention provides a method of suppressing the maternal electrocardiogram (MECG) to extract the FECG. 
   b. Related Art 
   The problem of detecting and separating a desired biomedical signal corrupted by other periodic interference and random noise signals is of extreme importance in medicine. Examples include His Purkinje System Electrogram (HPSE) described by H. Al-Nashash, S. Kelly S. and D. Taylor in “Beat-to-Beat detection of His-Purkinje system signals using adaptive filters”, Med &amp; Biol. Eng. &amp; Comp., vol. 26, pp. 117-125,1988; and Ventricular Late Potentials (VLP)° and the Diaphragmatic Electromyogram (EMGdi) described by H. Al-Nashash, S. Kelly S. and D. Taylor in “Noninvasive beat-to-beat detection of ventricular late potentials”, Med. &amp; Biol. Eng. &amp; Comp., vol. 27, pp. 64-68,1989. 
   There have been different methods proposed for such extraction including adaptive filters for example in B. Widrow et al., “Adaptive noise canceling: Principle and application,” Proc. IEEE, vol. 63, pp. 1692-1716, Dec. 1975 or E. R. Ferrara and B. Widrow, “Fetal electrocardiogram enhancement by time-sequenced adaptive filtering.” IEEE Trans Biomed Eng., vol. 29, pp. 451-460, 1982., correlation techniques, for example S. Abbond, A. Alaluf, S. Finav, and D. Sadeh, “Real time abdominal fetal ECG recording using hardware correlator,” Comput. Biol. Med., vol. 22, pp. 32-335, 1992, singular-value decomposition (SVD) for example in D. Callaerts, B. D. Moor, J. Vandewalle, and W. Sansen, “Comparison of SVD methods to extract the fetal electrocardiogram from coetaneous electrode signals,” Med Biol. Eng. Comput., vol. 28, pp. 217-224, 1990, blind source separation, for example, V. Zarzoso and A. Nandiin “Noninvasive fetal electrocardiogram extraction: Blind separation versus adaptive noise cancellation”, IEEE Trans. Biomed Eng., vol. 48, no. 1, pp. 12-18, 2001; or L. De Laflaauwer. B. De Moor and J. Vandewalle, “Fetal Electrocardiograph Extraction by Source Subspace Separation”, Proc. IEEE SP/ATHOS Workshop on HOS, Girona, Spain, pp. 134-138, 1995, and wavelet transform, for example, A. Khamene and S. Negahdaripour. “A New Method for the Extraction of Fetal ECG from the Composite Abdominal Signal”, IEEE Trans. Biomed. Eng., vol. 47, no. 4, pp. 507-516, 2000 to mention a few. Despite the reported successes of these methods, they are still not used clinically at a large scale. 
   Other prior art techniques include U.S. Pat. No. 5,042,499 which describes a fetal heart rate monitor which is used to monitor weak fetal electrocardiogram signals in the presence of strong interfering material ECG complexes, general random background muscle noise, and 60 Hz power line noise. This invention uses an electrocardiographic adaptive cancellation process to cancel the maternal ECG component from an abdominal ECG signal recorded from a pregnant subject&#39;s abdomen signal. 
   U.S. Pat. No 5,372,139 describes a method for suppressing a maternal electrocardiogram signal from a fetal electrocardiogram signal obtained with invasive and non-invasive techniques using an almost pure maternal electrocardiogram signal as a trigger. The MECG suppression is achieved by using a substantially pure MECG trace as a trigger signal to calculate an MECG complex template which is subtracted from the composite signal using the MECG trace to align the MECG complex template with each MECG complex (16) in the composite signal. 
   U.S. Pat. No 6,751,498 describes an apparatus and method for non-invasive, passive fetal heart monitoring. An apparatus and method for fetal heart and maternal heart and uterine monitoring is described which acquires waveforms indicative of the mother&#39;s heart beat from sensors located at or near the mother&#39;s chest, and waveforms indicative of the combined maternal and fetal heart beats from abdominal sensors located on the mother&#39;s abdomen, lower back, or both. The signals from the abdominal sensors are divided into a plurality of channels. An adaptive signal processing filter (ASPF) algorithm or other suitable algorithm is then used to cancel the estimated maternal waveform from each channel derived from the abdominal sensors. The system then selects from the resulting waveforms at least one waveform to serve as the reference fetal waveform. The reference waveform is then processed against the other abdominal waveforms preferably using the ASPF algorithm again to form an enhanced fetal signal that is a representation of the fetus&#39;s electrocardiogram. 
   U.S. Pat. No 5,123,420 describes a method and apparatus for processing heart rate traces in a fetal monitor. A fetal monitor is disclosed which is capable of recording the heart rate trace, preferably the beat-to-beat heart rate trace, of a fetus and a second heart rate trace of the mother or of a second fetus. Coincidence between the heart rate traces is detected by means of a direct or indirect comparison of the two traces and comparison of the difference with a predefined or adaptive limit and a warning signal is generated if coincidence is detected. 
   SUMMARY OF THE INVENTION 
   The invention provides an improved method of estimating a desired signal from a first signal and a second signal, where the first signal is a mixed signal considered to comprise the desired signal combined with a distorted version of the second signal, in particular where the desired signal is a fetal electrocardiogram signal, the first signal is an abdominal electrocardiogram signal and the second signal is a thoracic electrocardiogram signal. The abdominal electrocardiogram signal is considered to comprise the fetal electrocardiogram combined with a distorted version of a maternal electrocardiogram signal and the thoracic electrocardiogram signal is assumed to accurately represent the maternal electrocardiogram signal. 
   The nonlinear distortion which the maternal electrocardiogram signal undergoes is modelled using polynomial networks. An advantage of this approach is that the method is non iterative. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The invention will now be described, by way of example only with reference to the accompanying drawings: in which 
       FIG. 1  illustrates schematically an outline of the method of the present invention; 
       FIG. 2  is a graph illustrating a maternal electrocardiogram signal plotted against an abdominal maternal electrocardiogram signal; 
       FIG. 3  is a flow chart illustrating the method of the present invention: 
       FIG. 4  is a flow chart illustrating the transforming method step in more detail; 
       FIG. 5  is a graph illustrating a transformed maternal electrocardiogram signal plotted against an abdominal maternal electrocardiogram signal; and 
       FIG. 6  is a graph comparing the plots shown in  FIG. 4  and  FIG. 2 . 
       FIGS. 7 ,  8  and  9  illustrate results of applying the method of the invention to various pair of signals 
       FIG. 10  illustrates the result of applying a prior art technique to the signals of  FIG. 7 . 
   

   DETAILED DESCRIPTION 
   Referring now to  FIG. 1 , w(n)represents an abdominal ECG signal and x(n)represents a thoracic ECG signal which is assumed to represent accurately the maternal ECG. 
   The abdominal signal w(n) may be considered to be the sum of a distorted version of the maternal ECG {circumflex over (x)}(n) and a noisy version of the fetal ECG ŝ(n) such that:
 
 w ( n )= {circumflex over (x)} ( n )+ ŝ ( n )
 
 ŝ ( n )= s ( n )+η( n )
 
{circumflex over (x)}(n) is distorted due to the fact that the signal is measured far away from its source (the mother&#39;s heart), and consequently it encounters some nonlinear transformation as it travels to the abdominal area ie:
 
 {circumflex over (x)} ( n )= T ( x ( n ))
 
   The problem would be simple if the transformation T was linear (i.e. if all that x(n) encounters by traveling from the heart to the abdominal areas was time delay and attenuation). In that case, x(n)could be aligned with w(n) via correlation and the signal ŝ(n) could be extracted by simply subtracting the aligned x(n) from w(n). Similarly, the problem would be easily solved in the frequency domain by spectral subtraction if the spectra of {circumflex over (x)}(n) and ŝ(n) were non-overlapping. 
   Unfortunately, the transformation between x(n)and the MECG component in w(n) is highly nonlinear as shown in  FIG. 2 . It can be seen that there is both time warping and time varying scaling.  FIG. 2  shows a 360-sample record of normalized w(n) plotted against the corresponding normalized record of x(n). The large values at the extreme right and left of the figure correspond to the MECG peaks, and the scattered values around the horizontal line x(n)=0 and around the region |w(n)|ε[0.2 0.5] correspond to the FECG peaks. The nonlinearity of the mapping can easily be seen by the amount of deviation from the unity-slope line. The thoracic signal x(n) is predominantly maternal (MECG), and it is assumed that the fetal component FECG in it is negligible. It should be noted that a proper placement of the thoracic and abdomen electrodes would result in a clean estimate of the FECG such that ŝ(n)≈s(n). If the abdomen electrode is not placed low enough on the mothers abdominal area then the resulting noisy fetal ECG estimate can be cleaned by post filtering using wavelet denoising and/or nonlinear filtering such as median filtering. 
   In order to estimate {circumflex over (x)}(n)and extract s(n) it is necessary to approximate the transformation T. In other words, it is desirable to align x(n)with w(n)in a way where the MECG component is perfectly aligned with x(n). Then the MECG component can be removed from w(n) to leave ŝ(n). 
   The transformation T is generated using polynomial networks which are described in W. M. Campbell, K. T. Assaleh, and C. C. Broun, “Speaker recognition with polynomial classifiers”, IEEE Transactions on Speech and Audio Processing, vol. 10, no. 4, pp 205-212, 2002; K. T. Assaleh and W. M. Campbell, “Speaker Identification Using a Polynomial-based Classifier”, Proceedings of the fourth international Symposium on Signal processing and its Applications ISSPA &#39;99, Brisbane, Australia, August 1999; and W. M. Campbell and Khaled Assaleh, “Low-Complexity Small-Vocabulary Speech Recognition for Portable Devices”, Proceedings of the fourth International Symposium on Signal Processing and its Applications ISSPA &#39;99, Brisbane, Australia, August 1999. 
   A nonlinear mapping between x(n)and w(n)is determined by minimizing a mean-squared error criterion. This minimization occurs when the dominant component of w(n) (i.e. the MECG) is aligned with x(n). 
   The method for generating the transformation T is frame-based. The signals x(n)and w(n) are partitioned into contiguous frames of N-samples long. 
   The i th  frames of x(n) and w(n) are given by:
 
 x   i ( m )= x ( iN+m )
 
 w   i ( m )= w ( iN+m ), where 0 ≦m≦N−- 1
 
   A time derivative vector sequence is generated from the frame x,(m)and its J time derivatives. Such that (in matrix notation): 
   
     
       
         
           
             X 
             i 
           
           ⁡ 
           
             [ 
             
               
                 
                   
                     
                       x 
                       i 
                     
                     ⁡ 
                     
                       ( 
                       0 
                       ) 
                     
                   
                 
                 
                   
                     
                       
                         x 
                         . 
                       
                       i 
                     
                     ⁡ 
                     
                       ( 
                       0 
                       ) 
                     
                   
                 
                 
                   ⋯ 
                 
                 
                   
                     
                       x 
                       i 
                       
                         ( 
                         J 
                         ) 
                       
                     
                     ⁡ 
                     
                       ( 
                       0 
                       ) 
                     
                   
                 
               
               
                 
                   
                     
                       x 
                       i 
                     
                     ⁡ 
                     
                       ( 
                       1 
                       ) 
                     
                   
                 
                 
                   
                     
                       
                         x 
                         . 
                       
                       i 
                     
                     ⁡ 
                     
                       ( 
                       1 
                       ) 
                     
                   
                 
                 
                   ⋯ 
                 
                 
                   
                     
                       x 
                       i 
                       
                         ( 
                         J 
                         ) 
                       
                     
                     ⁡ 
                     
                       ( 
                       1 
                       ) 
                     
                   
                 
               
               
                 
                   ⋮ 
                 
                 
                   ⋮ 
                 
                 
                   ⋰ 
                 
                 
                   ⋮ 
                 
               
               
                 
                   
                     
                       x 
                       i 
                     
                     ⁡ 
                     
                       ( 
                       
                         N 
                         - 
                         1 
                       
                       ) 
                     
                   
                 
                 
                   
                     
                       
                         x 
                         . 
                       
                       i 
                     
                     ⁡ 
                     
                       ( 
                       
                         N 
                         - 
                         1 
                       
                       ) 
                     
                   
                 
                 
                   ⋯ 
                 
                 
                   
                     
                       x 
                       i 
                       
                         ( 
                         J 
                         ) 
                       
                     
                     ⁡ 
                     
                       ( 
                       
                         N 
                         - 
                         1 
                       
                       ) 
                     
                   
                 
               
             
             ] 
           
         
       
     
   
   The time derivative vector sequence X i  is then used to generate a Polynomial matrix P i  (comprising a vector sequence p(x i (m))) which is generated by polynomial expansion. 
   For example, for a two dimensional vector x i (m)=[x i (m) {dot over (x)} i (m)] and a second order polynomial, the vector is given by
 
 p ( x   i ( m )=└1 x   i ( m ) {dot over (x)}   i ( m ) x   i ( m ) {dot over (x)}   i ( m ) 2   x   i ( m ) {dot over (x)}   i ( m ) 2 ┘
 
   Therefore, by polynomial expansion higher order terms of the vector elements and all possible cross terms with their derivatives are introduced, which allow for highly nonlinear modeling of the problem in hand. 
   In general, for a K th  order polynomial, the polynomial basis terms are monomials of the form 
   
     
       
         
           
             
               ∏ 
               
                 J 
                 = 
                 0 
               
               J 
             
             ⁢ 
             
                 
             
             ⁢ 
             
               
                 ( 
                 
                   
                     x 
                     i 
                     
                       ( 
                       J 
                       ) 
                     
                   
                   ⁡ 
                   
                     ( 
                     m 
                     ) 
                   
                 
                 ) 
               
               
                 k 
                 i 
               
             
           
           , 
           
             
               where 
               ⁢ 
               
                   
               
               ⁢ 
               
                 
                   ∑ 
                   
                     j 
                     = 
                     0 
                   
                   J 
                 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   k 
                   j 
                 
               
             
             ≤ 
             K 
           
         
       
     
   
   The mapping of the input sequence x i (m) onto the desired output sequence w i (m) is achieved by using mean-squared error as the objective criterion such that: 
   
     
       
         
           
             c 
             i 
           
           = 
           
             
               
                 arg 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 min 
               
               
                 c 
                 i 
               
             
             ⁢ 
             
               
                  
                 
                   
                     
                       P 
                       i 
                     
                     ⁢ 
                     
                       c 
                       i 
                     
                   
                   - 
                   
                     w 
                     i 
                   
                 
                  
               
               2 
             
           
         
       
     
   
   Solving for c i  can be accomplished by applying the method of normal equations as described in B. Noble and 3. W. Daniel, “Applied Linear Algebra”, Third Ed., Prentice Hall, 1988 such that
 
P i   t P i c i =P i   t w i 
 
 c   i =( P   i   t   P   i ) −1   P   i   t   w   i 
 
   The vector c i  maps the input sequence to the output sequence such that the MECG component in w i (m) is optimally aligned with x i (m). This mapping is the optimal nonlinear approximation of the transformation T in the context of the polynomial networks such that:
 
 T ( x   i ( n ))= {circumflex over (x)}≈P   i   c   i 
 
   Consequently, the FECG component ŝ i =[ŝ i (0) ŝ i (1) . . . ŝ i (N−1]I the frame w i (m) can be extracted by subtracting the estimated MECG component from it as follows.
 
 ŝ   t   =w   t   −P   c   t 
 
     FIGS. 3 and 4  are flow charts illustrating the method of the present invention using the principles described above. 
   Referring to  FIG. 3  which shows an overview of the method of the present invention, at step  30  the abdominal electrocardiogram signal w is measured; at step  32  the thoracic electrocardiogram signal x is measured. At step  34  the thoracic electrocardiogram signal is transformed to provided a distorted version of the thoracic electrocardiogram signal {circumflex over (x)} such that the difference between the distorted version {circumflex over (x)} and the abdominal electrocardiogram signal w is minimised. At step  36  the distorted version of the thoracic electrocardiogram signal {circumflex over (x)} is subtracted from the abdominal electrocardiogram signal w to provide an estimate of the fetal electrocardiogram. 
   The transformation step  34  is described in more detail with reference to  FIG. 4 . At step  42  the abdominal electrocardiogram signal w is partitioned into a sequence of frames. At step  44  the thoracic electrocardiogram signal x is partitioned into a sequence of frames. Each frame w i  of the abdominal electrocardiogram signal corresponds in time to a frame x i  of the thoracic electrocardiogram signal. In the preferred embodiment of the invention each frame comprises 360 samples of the respective electrocardiogram signal 
   A matrix X i  is then generated from a frame of the thoracic electrocardiogram signal x i . In the preferred embodiment of the invention 3 time derivatives are used, thus the matrix X i  has is 360×4 values. 
   Polynomial expansion is than applied to the matrix X i  at step  47  to create a polynomial matrix P i . In the preferred embodiment of the invention a third order polynomial expansion is used and hence P i  has the dimensions 360×35. 
   At step  48  the polynomial matrix P i  is used together with the corresponding frame of the abdominal electrocardiogram signal w i  to generate a vector c i  in accordance with the equation c i =(P i   t P i ) −1 P i   t w i  as described above. It wil be appreciated that in the preferred embodiment of the Invention (P i   t P i ) −1  has the dimensions 35×35. 
   Vector c i  is then multiplied by the polynomial matrix P i  at step  49  to produce the required distorted version of the thoracic electrocardiogram signal {circumflex over (x)}. 
     FIGS. 5  shows a graph of w(n) against the separated MECG component {circumflex over (x)}(n). The larger values of both signals are perfectly aligned do it is possible to subtract the MECG component from w(n). 
     FIG. 6  show the graph of  FIG. 5  superimposed upon the graph of  FIG. 2  to illustrate alignment before and after the transformation. It will be appreciated that that the MECG component in w(n) is aligned while the FECG component is almost unchanged. 
     FIGS. 7 ,  8  and  9  show how the proposed technique extracts the FECG component from w(n) for three representative frames. 
   FIG.  7 . a  shows one frame of w(n)containing non-overlapping 3 FECG beats and two MECG beats. FIG.  7 . b  shows the output x(n)of the nonlinear mapping process as described above. FIG.  7 . c  depicts the extracted FECG signal which clearly shows the separated 3 beats while completely eliminating the presence of the MECG in w(n). 
   FIG.  8 . a, b  shows another frame of w(n) containing a partial overlap between one of the FECG beats with one MECG beat, and the corresponding frame of {circumflex over (x)}(n). FIG.  7 . c  depicts the extracted FECG signal which clearly shows the separated 3 beats while completely eliminating the presence of the MECG in w(n). 
   To demonstrate the power of the method of the invention in extracting the FECG, FIG.  9 . a  shows a third frame of w(n) containing full overlap between the first FECG and the MECG beats. This case represents the extreme case where the FECG is completely dominated by the MECG component to the point that the FECG beat is no longer visually distinguishable. FIG.  9 . c  shows that the algorithm was still successfully capable of extracting the FECG signal. 
   Independent component analysis (ICA) technique has been successfully used by several research groups for the extraction of FECG. See, for example, V. Zarzoso and A. Nandi, “Noninvasive fetal electrocardiogram extraction: Blind separation versus adaptive noise cancellation”, IEEE Trans. Biomed Eng., vol. 48, no, I, pp. 12-18, 2001. and B. Widrow et al., “Adaptive noise canceling: Principle and application,” Proc. IEEE, vol. 63, pp. 1692-1716, December 1975. 
   However, multiple leads were required for successful separation of the FECG. In practice and for monitoring purposes, it is highly desirable that the FECG is extracted using fewer leads. As illustrated in  FIGS. 7-9 , the method of the invention is capable of extracting FECG using two leads only. However, this is not possible with the ICA method described by Zaroso and Nandi. 
   FIG.  10 . a  shows one frame of w(n) containing non-overlapping three FECG beats and two MECG beats (same as in FIG.  7 . a ).  FIGS. 10   b , and  10   c  show the output of the ICA separated MECG and FECG respectively using two channels only. Clearly, the generated outputs are still composite, and a strong component of the MECG is still apparent in the extracted FECG. 
   It is to be recognized that various alterations, modifications, and/or additions may be introduced into the method described above without departing from the spirit or ambit of the present invention.