# Patent

## Title
Method and apparatus for predicting material hypertension during pregnancy using coherence analysis of material and fetal blood velocity waveforms

## Abstract

An embodiment of a method of predicting maternal hypertension during pregnancy is disclosed. The method may include measuring a maternal blood velocity waveform and measuring a fetal blood velocity waveform. The method may further include the calculation of a coherence value using a magnitude-squared coherence function between the maternal and fetal blood velocity waveforms and comparing the coherence value to a predetermined threshold value.

## Background
This non-provisional application claims priority to U.S. application No. 60/499,699 filed Sep. 2, 2003.


This application is funded from a Predoctoral Fellowship from the American Heart Association, grant number 0110005Y.


BACKGROUND OF THE INVENTION
1. Field of the Invention
This patent application is related to the field of signal processing. More particularly, the invention relates to the use of ultrasound Doppler measurements.
2. Related Art
Preeclampsia is a major cause of maternal mortality as well as perinatal morbidity and mortality. It is mainly characterized by hypertension and protienuria. Preeclampsia affects as many as one in ten of all pregnancies if its milder forms are counted. The incidence of preeclampsia is higher in women who are nulliparous, carry multiple fetuses, have chronic hypertension, diabetes, or if the woman has had previous preeclamptic pregnancies. Although clinical manifestations do not appear until the final three months, the foundations for preeclampsia are thought to be established in the first 10-12 weeks of gestation.
Maternal complications due to preeclampsia are associated with the vascular system. Mothers are at risk for intravascular coagulation, bleeding, and organ failure. Maternal blood volume is reduced and has altered distribution. The treatment of hypertension in pregnancy does not eliminate or reverse the process but does help prevent maternal cardiovascular complications, especially during labor and delivery.
Evidence exists to suggest that there is inadequate trophoblast secondary invasion in pregnancies that later developed preeclampsia. The adaptation to a limited supply of nutrients due to inadequate secondary trophoblast invasion in preeclampsia may lead to damages in the structure and the metabolism of a fetal life. Suboptimum supply of oxygen and nutrients may also lead to fetal intrauterine growth restriction. This, together with premature delivery, poses major threats to the fetus leading to various degrees of morbidity and even death. Long term follow up studies have demonstrated that babies who have suffered intrauterine growth retardation and those who are small and disproportionate at birth are more likely to develop hypertension, coronary heart disease and non-insulin-dependent diabetes in adult life. These observations suggest that fetuses are quite sensitive to the maternal placental function.
Various methods have been developed to detect preeclamptic pregnancies early in the gestational period. The use of ambulatory pressure monitoring for identifying hypertensive pregnancies has been used. Using this method, blood pressure monitoring can be performed in four-week intervals for 48 hours throughout gestation starting in the 16th week. Ambulatory pressure monitoring has shown significant results with specificity and sensitivity above 90%. However, the instrumentation required for pressure monitoring is expensive and reported lower acceptability and sleep disturbances among patients.
Doppler ultrasound blood velocity waveforms have also been employed for predicting preeclampsia. A uterine arterial pulsatility index (PI) has been studied as an indicator to the increased resistance to the placental flow during the 11th through 14th week and 23rd week of gestation in various studies using Doppler ultrasound measurements. Analysis of the Doppler ultrasound measurements has proven difficult, however. Current techniques to analyze the measurements have resulted in sensitivities of mean PI for preeclampsia at rates as low as 25-27% during 1-14 weeks and 41% at 23 weeks. Superior techniques for analyzing data from Doppler ultrasound measurements are greatly desired.
SUMMARY OF THE INVENTION
Embodiments of a method and apparatus for predicting maternal hypertension during pregnancy using coherence analysis of maternal and fetal blood velocity waveforms are disclosed.
An embodiment of a method for predicting maternal hypertension during pregnancy is disclosed. The method may include measuring a maternal blood velocity waveform and measuring a fetal blood velocity waveform. The method may further include the calculation of a coherence value using a magnitude-squared coherence function between the maternal and fetal blood velocity waveforms and comparing the coherence value to a predetermined threshold value.
An embodiment of an apparatus for predicting maternal hypertension during pregnancy is also disclosed. The apparatus may include a sensor for measuring a maternal blood velocity waveform and a sensor for measuring a fetal blood velocity waveform. The apparatus may further include a processor in communication with the sensors for calculating a coherence value using a magnitude-squared coherence function between the maternal and fetal blood velocity waveforms and comparing the coherence value to a predetermined threshold value.
Additional features and advantages of the invention will be apparent from the detailed description which follows, taken in conjunction with the accompanying drawings, which together illustrate, by way of example, features of embodiments of the present invention.



BRIEF DESCRIPTION OF THE DRAWINGS
The following drawings illustrate exemplary embodiments for carrying out the invention. Like reference numerals refer to like parts in different views or embodiments of the present invention in the drawings.
FIG. 1 is a flow diagram of a method for predicting maternal hypertension during pregnancy in accordance with an embodiment of the present invention.
FIG. 2 is a block diagram of an apparatus for predicting maternal hypertension during pregnancy in accordance with an embodiment of the present invention.
FIG. 3 shows the estimated spectrum of a segment of Doppler signal recorded from an 18 week old fetus.
FIG. 4 shows the integrated spectrum of the signal in FIG. 3.
FIG. 5 shows the distribution of the square-root of the MSC function between maternal and fetal arterial blood flow at the maternal heart rate with gestational age.
FIG. 6 shows the distribution of the threshold level with weighting constants.
FIG. 7 shows the distribution of the sensitivity and specificity with weighting constants.
FIG. 8 shows the threshold level and coherence distribution.
FIG. 9 shows a table containing test data.



Reference will now be made to the exemplary embodiments illustrated, and specific language will be used herein to describe the same. It will nevertheless be understood that no limitation or the scope of the invention is thereby intended.
DETAILED DESCRIPTION OF EXAMPLE EMBODIMENT(S)
FIG. 1 is a flow diagram of a method 100 of predicting maternal hypertension during pregnancy in accordance with an embodiment of the present invention. The method 100 may include measuring 102 a maternal blood velocity waveform and measuring 104 a fetal blood velocity waveform. The measuring steps 102, 104, may be performed sequentially in any order or in a presently preferred embodiment, simultaneously as shown in FIG. 1. This method may further include calculating 106 a coherence value using a magnitude-squared coherence (MSC) function between the maternal and fetal blood velocity waveforms. An operation in the method may further include comparing 108 the coherence value to a predetermined threshold value. Such a method may be contained in computer instructions suitable for implementation in computer hardware, software or firmware according to embodiments of the present invention.
FIG. 2 is a block diagram of an apparatus 200 for predicting maternal hypertension during pregnancy in accordance with an embodiment of the present invention. The apparatus 200 may include a maternal sensor 202 for measuring a maternal blood velocity waveform and a fetal sensor 204 for measuring a fetal blood velocity waveform. The sensors 202, 204 may be Doppler ultrasound sensors or any other sensor capable of detecting blood velocity waveforms. A processor 206 can be in communication with the sensors 202, 204 and be configured for calculating a coherence value using a magnitude-squared coherence function between the maternal and fetal blood velocity waveforms. The coherence value can be compared to a predetermined threshold value. Processor 206 may be a general-purpose microprocessor such as a central processing unit (CPU), a digital signal processor (DSP), a reconfigurable field programmable gate array (FPGA) such as a Xilinx programmable logic device, a general-purpose computer or any other computer hardware with or without software or firmware for performing the calculations specified herein.
Yet another embodiment of the present invention includes a processor 206 configured to receive maternal and fetal blood velocity waveforms as input and to calculate a coherence value using a magnitude-squared coherence function between the maternal and fetal blood velocity waveforms and comparing the coherence value to a predetermined threshold value. This embodiment could be configured to receive input from one or more Doppler ultrasound machines.
One embodiment of apparatus 200 contemplates the simultaneous measurement of maternal and fetal blood velocity waveforms. However, sequential and time delayed measurement of maternal and fetal blood velocity waveforms are also contemplated to be within the scope of the present invention and for use with the method and apparatus described above. Experimental results using sequential measurement of maternal and fetal blood velocity waveforms are described below.
The main advantages of this coherence analysis for the prediction of preeclampsia are: 1) it can be used to identify women at risk for hypertensive disorders at an earlier stage of pregnancy than currently possible, 2) it is non-invasive and simple to perform, and 3) it is highly sensitive.
Physiology Related to Development of Preeclampsia
Prior to implantation, the cells of the early embryo separate into two types, the inner cell mass and the outer cells. The inner cell mass becomes the fetus and its protective membranes, and the outer cells are collectively known as the trophoblast. The trophoblast is the only embryonic tissue that comes into contact with maternal tissues. Its main purpose is to form the placenta that juxtaposes the maternal and fetal circulations. A secondary function of the trophoblast is to engineer structural changes in the maternal blood vessels terminating in the placenta bed, so that they can carry a substantial volume of blood to the placenta. These vessels are known as spiral arteries due to their coiled shape.
The transformation of the spiral arteries takes place in two phases. During the primary trophoblast invasion within the first few weeks after implantation, the endometrial blood vessels are converted to the uteroplacental vessels. In the second phase, between 12 to 18 weeks of pregnancy, the trophoblast invades more deeply into the segments of the vessels within the uterine muscle and further increases their capacity. The secondary trophoblast invasion results in an increase in the vessel diameter and a decrease in the vessel wall thickness. When there is inadequate trophoblast invasion, the spiral artery remains narrow and its walls thicken. This results in reduced blood flow to the intervillos space. The inadequacy of secondary trophoblast invasion results in a spectrum of pregnancy complications including maternal hypertension and preeclampsia.
There is evidence that absence of or inadequate secondary trophoblast invasion can contribute to the development of hypertensive disorders later in human pregnancy. The functional relationship between the maternal and fetal blood flows may be different in normotensive and hypertensive pregnancies. Furthermore, the secondary trophoblast invasion in the placental function typically occurs by the 12th week of gestation. This can make it possible to detect hypertensive pregnancies through the analysis of maternal-fetal blood flows as early as the late first or the early second trimesters of pregnancy.
Doppler Ultrasound Techniques
Blood flow velocity waveforms can be measured from fetal arteries as early as the eighth week of gestation using transvaginal Doppler Ultrasound techniques. Ultrasonography can be a safe, noninvasive, and cost effective tool for monitoring the fetal cardiovascular system through imaging and blood velocity measurements. However, exposure to ultrasound beam for long durations of time may cause temperature increase and damage in fetal tissues, especially in the brain encased in the fetal skull. Consequently, long term monitoring of blood velocity waveforms is not recommended.
The blood velocity waveforms estimated by the use of the Doppler ultrasound techniques can be modeled as a single sinusoid embedded in white noise over short intervals of time. Such an inference can be justified when considering the autocovariance function of the signal, which matches the sinusoidal model. Therefore, the blood velocity waveforms can be modeled as a pure sinusoid embedded in white noise over short intervals.
In Doppler Ultrasonography, a beam of ultrasound is transmitted toward the blood flow whose velocity is being estimated, and it is backscattered mainly by the red blood cells of the flow. The frequency difference between the transmitted and backscattered waves is the Doppler shift, and is directly proportional to the relative velocity between the reflective interface and the receiver. The velocity is estimated using the relationship




  
    v
    =
    
      cf
      
        2
        ⁢
        
          f
          0
        
        ⁢
        cos
        ⁢
        
          
        
        ⁢
        θ
      
    
  
  ,




where c is the velocity of the ultrasound in the medium in meters/sec, f is the Doppler shift frequency in Hz, fo is the frequency of the incident ultrasound beam in Hz, and θ is the angle of incidence in radians.

Reconstruction of Maximum Blood Flow Velocity

Doppler insonation of the blood vessel produces a spectrum of frequencies. The maximum frequency of the Doppler spectrum, fmax is proportional to the maximum blood flow velocity. In current practice, the maximum frequency waveform is estimated as the largest frequency value at which the Doppler spectrum exceeds a predetermined threshold. The most common method for estimating the maximum frequency waveform takes the largest frequency value that exceeds a certain threshold level. There are four main variations of this method including (i) the percentile method, (ii) D'Alessio's threshold crossing method, (iii) the modified threshold crossing method, and (iv) the hybrid method.
The percentile method uses the integrated Doppler power spectrum defined as




  
    
      φ
      ⁡
      
        (
        f
        )
      
    
    =
    
      
        ∫
        0
        f
      
      ⁢
      
        
          S
          ⁡
          
            (
            λ
            )
          
        
        ⁢
        
          d
          λ
        
      
    
  
  ,




where S(f) is the power spectrum at frequency f. In this method the maximum frequency fmax is defined as the frequency at which





  
    
      φ
      ⁡
      
        (
        
          f
          max
        
        )
      
    
    =
    
      
        
          100
          -
          α
        
        100
      
      ⁢
      
        φ
        T
      
    
  
  ,




where ΦT is the total signal power and α(0<α<100) is a designed parameter whose optimum value depends on both the signal-to-noise ratio (SNR) and the signal bandwidth.

Both the D'Alessio's method and the modified threshold crossing method compare the estimated power spectrum in a sliding window of r successive frequency bins to a chosen threshold level. The highest frequency that exceeds the threshold is recorded for each window for which the threshold is exceeded in at least m consecutive bins. The maximum value of this sequence of recorded frequencies is selected as the maximum frequency fmax. The threshold level in the modified threshold crossing method is chosen as a multiple of the noise density estimated as the average value of the Doppler spectrum over high frequency values where the noise components predominate.
In the hybrid algorithm, the noise density is estimated as the gradient of the integrated spectrum in the high frequency range. The maximum frequency is then estimated as the frequency at which the straight line that passes through the maximum point of the integrated spectrum at a gradient greater than the estimated noise density intersects the integrated spectrum. These maximum frequency estimation methods depend on arbitrarily selected parameters, and the estimated root-mean-square (RMS) error varies significantly with both the SNR and the bandwidth of the signal and was greater than 4% of the fmax even at SNR values greater than 9 dB.
A geometric method, wherein maximum frequency estimation algorithms are based on the geometry of the integrated spectrum, worked better in experiments than both the percentile method and the modified threshold crossing method in experiments. The geometric method first estimates the mode frequency, the frequency at which the signal has maximum power. A reference line is then drawn connecting the points on the integrated spectrum at the estimated mode frequency and the highest frequency value (half the sampling frequency in digitized signals). The maximum frequency is then estimated as the frequency at which the Euclidean distance from this reference line to the spectrum is maximum. The geometric method has been shown to work well for Gaussian shaped spectra where the maximum frequency is relatively far away from the mode frequency, but it fails when the mode frequency is close to the maximum frequency or when the spectrum is relatively flat. It has also been shown that the estimated maximum frequency using the geometric method is susceptible to amplifier gain.
A modified geometrical method can also be implemented wherein the reference line is selected from the first point or origin of the integrated spectrum to its maximum point. The maximum frequency is estimated as the frequency corresponding to the maximum vertical distance from the reference line to the integrated spectrum. This method showed good performance in experiments. A proof of the optimality rotated integrated spectrum based methods is below.
Let

y(n)=x(n)+η(n)

represent the digitized Doppler signal at time n, where x(n) is the signal of interest and η(n) is a zero-mean and white noise process with variance ση2. Let fN denote the Nyquist frequency (half the sampling frequency) of the signal. Assuming that x(n) is band limited to the frequency fmax, where fmax<fN.

The integrated spectrum of the signal component x(n) is given by




  
    
      
        φ
        x
      
      ⁡
      
        (
        f
        )
      
    
    =
    
      
        ∫
        0
        f
      
      ⁢
      
        
          
            S
            x
          
          ⁡
          
            (
            λ
            )
          
        
        ⁢
        
          d
          λ
        
      
    
  
  ,




where Sx(f) is its power spectrum at frequency f. The integrated spectrum of y(n) is given by

φY(f)=φX(f)+ση2f. 

Since x(n) is band-limited to fmax, for frequencies f>fmax,

φY(f)=φX(fmax)+ση2f. 

An example of the estimated spectrum of a portion of a Doppler signal recorded from the umbilical artery of an 18 week old fetus and the corresponding integrated spectrum are depicted in FIGS. 3 and 4 respectively. The power spectrum was estimated as the periodogram of a Hann-windowed 512-sample segment of data. The FFT size was 1024 samples.

The Gradient of the integrated spectrum




  
    ∂
    
      
        φ
        Y
      
      ⁡
      
        (
        f
        )
      
    
  
  
    ∂
    f
  




in the frequency range fmax<f<fN corresponds to the instantaneous noise density. The average gradient of the integrated spectrum within the signal bandwidth is greater than the noise density. An estimate of the maximum Doppler frequency can be found using the following method. The integrated Doppler spectrum can be rotated over the frequency to obtain

φrot(f)=φY(f)cos β−f sin β,

where β is the angle of rotation. For 0<β<π/2, dividing both sides by cos β results in

{circumflex over (φ)}rot(f)=φY(f)−f tan β.

It is straightforward to show that φrot(f) and {circumflex over (φ)}rot(f) peak at the same frequency.

In the modified geometric method, a straight line can be drawn connecting the origin to the maximum value of the integrated spectrum as shown by the segment AB in FIG. 4. The frequency at which the vertical distance from this straight line to the integrated spectrum is maximum (CD) is the estimate of the maximum Doppler shift. Therefore, the vertical distance between the integrated spectrum and the straight line AB at any frequency f is given by

{circumflex over (φ)}rot(f)=φY(f)−f tan β,

where





  
    tan
    ⁢
    
      
    
    ⁢
    β
  
  =
  
    
      
        
          φ
          Y
        
        ⁡
        
          (
          
            f
            N
          
          )
        
      
      
        f
        N
      
    
    =
    
      
        
          
            φ
            X
          
          ⁡
          
            (
            
              f
              max
            
            )
          
        
        
          f
          N
        
      
      +
      
        
          σ
          η
          2
        
        .
      
    
  




It is clear that the modified geometric method is a special case of the rotated integrated spectrum in which the angle of rotation is given by the above equation. The maximum vertical distance occurs at the maximum Doppler shift for angles that satisfy the following smoothness constraint:





  
    
      
        
          φ
          x
        
        ⁡
        
          (
          
            f
            max
          
          )
        
      
      
        f
        N
      
    
    <
    
      
        
          
            φ
            x
          
          ⁡
          
            (
            
              f
              max
            
            )
          
        
        -
        
          
            φ
            x
          
          ⁡
          
            (
            f
            )
          
        
      
      
        
          f
          max
        
        -
        f
      
    
  
  ;
  
    0
    ≤
    f
    ≤
    
      
        f
        max
      
      .
    
  




This condition implies that the average spectrum of the input signal over the whole frequency range is smaller than the average spectrum computed in any interval in which the signal is present. This is a smoothness condition which constrains the signal spectrum to be relatively flat and not to exhibit significant dynamics over its bandwidth. For lowpass signals such that fN>>fmax this assumption is satisfied easily in most practical situations.

The maximum velocity reconstruction algorithm is based on the following lemma.
Lemma: for signal spectra satisfying φY(f)=φX(f)+ση2f the maximum value of {circumflex over (φ)}rot(f) defined above as

{circumflex over (φ)}rot(f)=φY(f)−f tan β

occurs at the maximum Doppler shift for any angle of rotation β that satisfies





  
    σ
    η
    2
  
  <
  
    tan
    ⁢
    
      
    
    ⁢
    β
  
  ≤
  
    
      
        
          φ
          X
        
        ⁡
        
          (
          
            f
            max
          
          )
        
      
      
        f
        N
      
    
    +
    
      
        σ
        η
        2
      
      .
    
  




The lower limit of tan β is chosen such that it is greater than the noise density. This ensures that the rotated integrated spectrum is no longer a monotonically increasing function. Theoretically, fmax can take any value between 0 and fN. Therefore, the maximum limit of the angle of rotation can be chosen such that the reference line passing through the origin and drawn at an angle β to the frequency axis does not intersect the integrated spectrum at any frequency value less than fN. The maximum Doppler frequency can then be estimated as the frequency at which {circumflex over (φ)}rot(f) peaks. The lemma can be proved by contradiction:

As can be seen from




  
    
      tan
      ⁢
      
        
      
      ⁢
      β
    
    =
    
      
        
          
            φ
            Y
          
          ⁡
          
            (
            
              f
              N
            
            )
          
        
        
          f
          N
        
      
      =
      
        
          
            
              φ
              Y
            
            ⁡
            
              (
              
                f
                max
              
              )
            
          
          
            f
            N
          
        
        +
        
          σ
          η
          2
        
      
    
  
  ,




the modified geometric method rotates the integrated spectrum by the maximum allowed angle in





  
    σ
    η
    2
  
  <
  
    tan
    ⁢
    
      
    
    ⁢
    β
  
  ≤
  
    
      
        
          φ
          X
        
        ⁡
        
          (
          
            f
            max
          
          )
        
      
      
        f
        N
      
    
    +
    
      
        σ
        η
        2
      
      .
    
  




Consequently, the frequency estimation using this approach can provide accurate results. This, combined with the ease of implementing the rotation, makes the maximum allowed rotation an excellent choice for maximum velocity estimation.

Reconstruction of Mean Blood Velocity Waveforms

In many applications, the power-weighted mean frequency of a signal x(n) defined as




  
    
      ω
      _
    
    =
    
      
        
          ∫
          
            ω
            min
          
          
            ω
            max
          
        
        ⁢
        
          
            
              ωS
              x
            
            ⁡
            
              (
              ω
              )
            
          
          ⁢
          
            
          
          ⁢
          
            d
            ω
          
        
      
      
        
          ∫
          
            ω
            min
          
          
            ω
            max
          
        
        ⁢
        
          
            
              S
              x
            
            ⁡
            
              (
              ω
              )
            
          
          ⁢
          
            
          
          ⁢
          
            d
            ω
          
        
      
    
  
  ,




where SX(ω) is the power spectrum of x(n) at the frequency ω, and ωmin and ωmax are respectively the minimum and maximum frequencies of the input signal. In current practice, the mean frequency estimation involves three steps: (i) estimate the power spectrum of the input signal, (ii) evaluate the minimum and the maximum frequency values, and (iii) compute the mean frequency using the above equation. Generally, the power spectrum has been estimated using a fast Fourier transform (FFT) algorithm. Once the spectrum is estimated, a threshold level is chosen manually such that it is larger than the noise level. The lowest and the highest frequency values at which the spectral powers equal the threshold level are selected as the minimum and the maximum signal frequencies.

A method which may be superior to using an FFT is the Pisarenko harmonic decomposition algorithm (PHD), wherein the single frequency approximation for a narrowband lowpass signal embedded in white noise using PHD is approximately the power-weighted mean frequency of the signal. In preliminary analysis, only signals whose bandwidth and maximum frequency were much smaller than the sampling frequency of thee digitized signal were used. However, experimental results show that the approximation is valid for signals up to one-fourth of the sampling frequency.
The estimation for complex-valued signals will first be shown and then it will be extended to real-valued signals. For complex-valued signals, let

y(n)=x(n)+η(n),

represent the input signal, where x(n) is the signal of interest and η(n) is an additive circular white noise process with zero-mean value and variance ση2. Here, x(n) is a zero-mean, complex-valued narrowband signal whose bandwidth is much smaller than the sampling frequency. The covariance function rX(k) of x(n) is related to its power spectrum SX(ω) through the expression





  
    
      
        r
        X
      
      ⁡
      
        (
        k
        )
      
    
    =
    
      
        1
        
          2
          ⁢
          π
        
      
      ⁢
      
        
          ∫
          
            ω
            min
          
          
            ω
            max
          
        
        ⁢
        
          
            
              S
              X
            
            ⁡
            
              (
              ω
              )
            
          
          ⁢
          
            e
            
              
                jω
                ⁢
                
                  
                
                ⁢
                k
              
              ⁢
              
                
              
            
          
          ⁢
          
            d
            ω
          
        
      
    
  
  ,




where ωmin and ωmax are respectively, the minimum and the maximum frequencies at which the signal x(n) is present. Consider the 2×2-element covariance matrix of the input signal given by





  
    
      R
      Y
    
    =
    
      [
      
        
          
            
              
                
                  r
                  X
                
                ⁡
                
                  (
                  0
                  )
                
              
              +
              
                σ
                η
                2
              
            
          
          
            
              
                r
                X
              
              ⁡
              
                (
                1
                )
              
            
          
        
        
          
            
              
                r
                X
                *
              
              ⁡
              
                (
                1
                )
              
            
          
          
            
              
                
                  r
                  X
                
                ⁡
                
                  (
                  0
                  )
                
              
              +
              
                σ
                η
                2
              
            
          
        
      
      ]
    
  
  ,




where r*X(1) represents the complex conjugate of rX(1). The eigenvalues of RY are

λ1=rX(0)+ση2+|rX(1)|,

and

λ2=rX(0)+ση2−|rX(1)|.

The eigenvector corresponding to the smallest eigenvalue λ2 is





  
    g
    =
    
      
        m
        ⁡
        
          [
          
            
              -
              
                
                  
                    r
                    X
                  
                  ⁡
                  
                    (
                    1
                    )
                  
                
                
                  
                  
                    
                      r
                      X
                    
                    ⁡
                    
                      (
                      1
                      )
                    
                  
                  
                
              
            
            ⁢
            1
          
          ]
        
      
      T
    
  
  ,




where m is an arbitrary constant.

Now, following the method employed in PHD to estimate the frequency of a single sinusoid that best represents RY, let w=[1 ej{circumflex over (ω)}]H, where {circumflex over (ω)} is the value of the single frequency that satisfies the condition wHg=0 and (·)H and (·)H denotes the Hermitian of the matrix or vector (·). It can be easily shown that wHg can take the value of zero for an appropriate choice of {circumflex over (ω)}. Substituting for w and g we get




  
    
      
        [
        
          1
          ⁢
          
            
          
          ⁢
          
            e
            
              j
              ⁢
              
                ω
                ^
              
            
          
        
        ]
      
      ⁡
      
        [
        
          
            -
            
              
                
                  r
                  x
                
                ⁡
                
                  (
                  1
                  )
                
              
              
                
                
                  
                    r
                    x
                  
                  ⁡
                  
                    (
                    1
                    )
                  
                
                
              
            
          
          ⁢
          1
        
        ]
      
    
    T
  
  =
  
    
      
        -
        
          
            
              r
              x
            
            ⁡
            
              (
              1
              )
            
          
          
            
            
              
                r
                x
              
              ⁡
              
                (
                1
                )
              
            
            
          
        
      
      +
      
        e
        
          j
          ⁢
          
            ω
            ^
          
        
      
    
    =
    0.
  




From this we can write





  
    ω
    ^
  
  =
  
    ∠
    ⁢
    
      
        
          
            r
            x
          
          ⁡
          
            (
            1
            )
          
        
        
          
          
            
              r
              x
            
            ⁡
            
              (
              1
              )
            
          
          
        
      
      .
    
  




The frequency represented by the phase angle of the autocovariance value at lag one of a complex-valued signal is the power weighted mean frequency of the signal. Therefore, the final equation above represents the power weighted mean frequency of the input signal x(n).
Considering the input signal x(n) with both narrowband and lowpass characteristics, for real-valued signals the covariance function rX(k) of x(n) is related to the power spectrum SX(ω) by




  
    
      
        r
        X
      
      ⁡
      
        (
        k
        )
      
    
    =
    
      
        1
        π
      
      ⁢
      
        
          ∫
          
            ω
            min
          
          
            ω
            max
          
        
        ⁢
        
          
            
              S
              X
            
            ⁡
            
              (
              ω
              )
            
          
          ⁢
          
            cos
            ⁡
            
              (
              
                k
                ⁢
                
                  
                
                ⁢
                ω
              
              )
            
          
          ⁢
          
            d
            ω
          
        
      
    
  
  ,




where ωmin and ωmax are respectively, the minimum and the maximum frequencies of the signal. Since each real-valued sinusoid can be represented by two complex-valued sinusoids, 3×3-element covariance matrices are considered. The 3×3-element covariance matrix of y(n) is





  
    
      R
      Y
    
    =
    
      [
      
        
          
            
              
                
                  r
                  x
                
                ⁡
                
                  (
                  0
                  )
                
              
              +
              
                σ
                η
                2
              
            
          
          
            
              
                r
                x
              
              ⁡
              
                (
                1
                )
              
            
          
          
            
              
                r
                x
              
              ⁡
              
                (
                2
                )
              
            
          
        
        
          
            
              
                r
                x
              
              ⁡
              
                (
                1
                )
              
            
          
          
            
              
                
                  r
                  x
                
                ⁡
                
                  (
                  0
                  )
                
              
              +
              
                σ
                η
                2
              
            
          
          
            
              
                r
                x
              
              ⁡
              
                (
                1
                )
              
            
          
        
        
          
            
              
                r
                x
              
              ⁡
              
                (
                2
                )
              
            
          
          
            
              
                r
                x
              
              ⁡
              
                (
                1
                )
              
            
          
          
            
              
                
                  r
                  x
                
                ⁡
                
                  (
                  0
                  )
                
              
              +
              
                σ
                η
                2
              
            
          
        
      
      ]
    
  
  ,




and its eigenvalues are




  
    
      
        λ
        1
      
      =
      
        
          
            r
            x
          
          ⁡
          
            (
            0
            )
          
        
        +
        
          σ
          η
          2
        
        +
        
          
            1
            2
          
          ⁢
          
            (
            
              
                
                  r
                  x
                
                ⁡
                
                  (
                  2
                  )
                
              
              +
              
                
                  
                    
                      r
                      x
                      2
                    
                    ⁡
                    
                      (
                      2
                      )
                    
                  
                  +
                  
                    8
                    ⁢
                    
                      
                        r
                        x
                        2
                      
                      ⁡
                      
                        (
                        1
                        )
                      
                    
                  
                
              
            
            )
          
        
      
    
    ,
    
      

    
    ⁢
    
      
        λ
        2
      
      =
      
        
          
            r
            x
          
          ⁡
          
            (
            0
            )
          
        
        +
        
          σ
          η
          2
        
        -
        
          
            r
            x
          
          ⁡
          
            (
            2
            )
          
        
      
    
    ,
    
      

    
    ⁢
    and
  




  
    
      λ
      3
    
    =
    
      
        
          r
          x
        
        ⁡
        
          (
          0
          )
        
      
      +
      
        σ
        η
        2
      
      +
      
        
          1
          2
        
        ⁢
        
          
            (
            
              
                
                  r
                  x
                
                ⁡
                
                  (
                  2
                  )
                
              
              -
              
                
                  
                    
                      r
                      x
                      2
                    
                    ⁡
                    
                      (
                      2
                      )
                    
                  
                  +
                  
                    8
                    ⁢
                    
                      
                        r
                        x
                        2
                      
                      ⁡
                      
                        (
                        1
                        )
                      
                    
                  
                
              
            
            )
          
          .
        
      
    
  



Considering lowpass signals such that |rX(1)|>rX(2), the eigenvector corresponding to the smallest eigenvalue λ3 is





  
    g
    =
    
      
        m
        ⁡
        
          [
          
            1
            -
            
              
                
                  4
                  ⁢
                  
                    
                      r
                      x
                    
                    ⁡
                    
                      (
                      1
                      )
                    
                  
                
                
                  
                    
                      
                        
                          r
                          x
                          2
                        
                        ⁡
                        
                          (
                          2
                          )
                        
                      
                      +
                      
                        8
                        ⁢
                        
                          
                            r
                            x
                            2
                          
                          ⁡
                          
                            (
                            1
                            )
                          
                        
                      
                    
                  
                  -
                  
                    
                      r
                      x
                    
                    ⁡
                    
                      (
                      2
                      )
                    
                  
                
              
              ⁢
              1
            
          
          ]
        
      
      T
    
  
  ,




where m is an arbitrary constant.

The solution to the relationship wHg=0, where w=[1 ej{circumflex over (ω)} ej2{circumflex over (ω)}]H yields




  
    cos
    ⁢
    
      
    
    ⁢
    
      ω
      ^
    
  
  =
  
    
      
        2
        ⁢
        
          
            r
            x
          
          ⁡
          
            (
            1
            )
          
        
      
      
        
          
            
              
                r
                x
                2
              
              ⁡
              
                (
                2
                )
              
            
            +
            
              8
              ⁢
              
                
                  r
                  x
                  2
                
                ⁡
                
                  (
                  1
                  )
                
              
            
          
        
        -
        
          
            r
            x
          
          ⁡
          
            (
            2
            )
          
        
      
    
    .
  




Substituting for rX2(2)+8rX2(1) gives





  
    
      
        
          
            
              r
              X
              2
            
            ⁡
            
              (
              2
              )
            
          
          +
          
            8
            ⁢
            
              
                r
                X
                2
              
              ⁡
              
                (
                1
                )
              
            
          
        
        =
        
        ⁢
        
          
            
              (
              
                
                  ∫
                  
                    ω
                    min
                  
                  
                    ω
                    max
                  
                
                ⁢
                
                  
                    
                      S
                      X
                    
                    ⁡
                    
                      (
                      ω
                      )
                    
                  
                  ⁢
                  cos
                  ⁢
                  
                    
                  
                  ⁢
                  2
                  ⁢
                  
                    
                  
                  ⁢
                  ω
                  ⁢
                  
                    d
                    ω
                  
                
              
              )
            
            2
          
          +
        
      
    
  
  
    
      
        
        ⁢
        
          
            8
            ⁢
            
              
                (
                
                  
                    ∫
                    
                      ω
                      min
                    
                    
                      ω
                      max
                    
                  
                  ⁢
                  
                    
                      
                        S
                        X
                      
                      ⁡
                      
                        (
                        ω
                        )
                      
                    
                    ⁢
                    cos
                    ⁢
                    
                      
                    
                    ⁢
                    ω
                    ⁢
                    
                      d
                      ω
                    
                  
                
                )
              
              2
            
          
          ,
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            
              ∫
              
                ω
                min
              
              
                ω
                max
              
            
            ⁢
            
              
                ∫
                
                  ω
                  min
                
                
                  ω
                  max
                
              
              ⁢
              
                
                  
                    S
                    X
                  
                  ⁡
                  
                    (
                    ω
                    )
                  
                
                ⁢
                
                  
                    S
                    X
                  
                  ⁡
                  
                    (
                    λ
                    )
                  
                
                ⁢
                cos
                ⁢
                
                  
                
                ⁢
                2
                ⁢
                
                  
                
                ⁢
                ω
                ⁢
                
                  
                
                ⁢
                cos
                ⁢
                
                  
                
                ⁢
                2
                ⁢
                
                  
                
                ⁢
                λ
                ⁢
                
                  d
                  ω
                
                ⁢
                
                  d
                  λ
                
              
            
          
          +
        
      
    
  
  
    
      
        
        ⁢
        
          
            8
            ⁢
            
              
                ∫
                
                  ω
                  min
                
                
                  ω
                  max
                
              
              ⁢
              
                
                  ∫
                  
                    ω
                    min
                  
                  
                    ω
                    max
                  
                
                ⁢
                
                  
                    
                      S
                      X
                    
                    ⁡
                    
                      (
                      ω
                      )
                    
                  
                  ⁢
                  
                    
                      S
                      X
                    
                    ⁡
                    
                      (
                      λ
                      )
                    
                  
                  ⁢
                  cos
                  ⁢
                  
                    
                  
                  ⁢
                  ω
                  ⁢
                  
                    
                  
                  ⁢
                  cos
                  ⁢
                  
                    
                  
                  ⁢
                  λ
                  ⁢
                  
                    d
                    ω
                  
                  ⁢
                  
                    d
                    λ
                  
                
              
            
          
          ,
        
      
    
  
  
    
      
        =
        
        ⁢
        
          
            ∫
            
              ω
              min
            
            
              ω
              max
            
          
          ⁢
          
            
              ∫
              
                ω
                min
              
              
                ω
                max
              
            
            ⁢
            
              
                
                  S
                  X
                
                ⁡
                
                  (
                  ω
                  )
                
              
              ⁢
              
                
                  S
                  X
                
                ⁡
                
                  (
                  λ
                  )
                
              
              ⁢
              
                (
                
                  (
                  
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      
                        cos
                        2
                      
                      ⁢
                      ω
                    
                    -
                    1
                  
                  )
                
              
            
          
        
      
    
  
  
    
      
        
          
            
            ⁢
            
              
                (
                
                  
                    2
                    ⁢
                    
                      
                    
                    ⁢
                    
                      cos
                      2
                    
                    ⁢
                    λ
                  
                  -
                  1
                
                )
              
              +
              
                8
                ⁢
                
                  
                
                ⁢
                cos
                ⁢
                
                  
                
                ⁢
                ω
                ⁢
                
                  
                
                ⁢
                cos
                ⁢
                
                  
                
                ⁢
                λ
              
            
            )
          
          ⁢
          
            d
            ω
          
          ⁢
          
            d
            λ
          
        
        ,
      
    
  
  
    
      
        =
        
        ⁢
        
          
            ∫
            
              ω
              min
            
            
              ω
              max
            
          
          ⁢
          
            
              ∫
              
                ω
                min
              
              
                ω
                max
              
            
            ⁢
            
              
                
                  S
                  X
                
                ⁡
                
                  (
                  ω
                  )
                
              
              ⁢
              
                
                  S
                  X
                
                ⁡
                
                  (
                  λ
                  )
                
              
              ⁢
              
                (
                
                  
                    4
                    ⁢
                    
                      
                    
                    ⁢
                    
                      cos
                      2
                    
                    ⁢
                    ω
                    ⁢
                    
                      
                    
                    ⁢
                    
                      cos
                      2
                    
                    ⁢
                    λ
                  
                  -
                
              
            
          
        
      
    
  
  
    
      
        
          
          ⁢
          
            
              2
              ⁢
              
                (
                
                  
                    
                      cos
                      2
                    
                    ⁢
                    ω
                  
                  +
                  
                    
                      cos
                      2
                    
                    ⁢
                    λ
                  
                
                )
              
            
            +
            
              8
              ⁢
              
                
              
              ⁢
              cos
              ⁢
              
                
              
              ⁢
              ω
              ⁢
              
                
              
              ⁢
              cos
              ⁢
              
                
              
              ⁢
              λ
            
            +
            1
          
          )
        
        ⁢
        
          d
          ω
        
        ⁢
        
          
            d
            λ
          
          .
        
      
    
  




For narrowband signals for which (ωmax−ωmin) is small, (cos ω−cos λ)2≈0. Therefore, for ωmin≦ω, λ≦ωmax 

cos2 ω+cos2 λ≈2 cos ω cos λ,

and

8 cos ω cos λ−2(cos2 ω+cos2 λ)≈2(cos2 ω+cos2 λ).

Substitution yields





  
    cos
    ⁢
    
      
    
    ⁢
    
      ω
      ^
    
  
  ≈
  
    
      
        
          ∫
          
            ω
            min
          
          
            ω
            max
          
        
        ⁢
        
          
            
              S
              X
            
            ⁡
            
              (
              ω
              )
            
          
          ⁢
          cos
          ⁢
          
            
          
          ⁢
          ω
          ⁢
          
            d
            ω
          
        
      
      
        
          ∫
          
            ω
            min
          
          
            ω
            max
          
        
        ⁢
        
          
            
              S
              X
            
            ⁡
            
              (
              ω
              )
            
          
          ⁢
          
            d
            ω
          
        
      
    
    .
  




That is





  
    cos
    ⁢
    
      
    
    ⁢
    
      ω
      ^
    
    ⁢
    
      
        ∫
        
          ω
          min
        
        
          ω
          max
        
      
      ⁢
      
        
          
            S
            X
          
          ⁡
          
            (
            ω
            )
          
        
        ⁢
        
          d
          ω
        
      
    
  
  ≈
  
    
      ∫
      
        ω
        min
      
      
        ω
        max
      
    
    ⁢
    
      
        
          S
          X
        
        ⁡
        
          (
          ω
          )
        
      
      ⁢
      cos
      ⁢
      
        
      
      ⁢
      ω
      ⁢
      
        
          d
          ω
        
        .
      
    
  




Dividing both sides by cos {circumflex over (ω)} results in





  
    
      ∫
      
        ω
        min
      
      
        ω
        max
      
    
    ⁢
    
      
        
          S
          X
        
        ⁡
        
          (
          ω
          )
        
      
      ⁢
      
        d
        ω
      
    
  
  ≈
  
    
      ∫
      
        ω
        min
      
      
        ω
        max
      
    
    ⁢
    
      
        
          S
          X
        
        ⁡
        
          (
          ω
          )
        
      
      ⁢
      
        
          cos
          ⁢
          
            
          
          ⁢
          ω
        
        
          cos
          ⁢
          
            
          
          ⁢
          
            ω
            ^
          
        
      
      ⁢
      
        
          d
          ω
        
        .
      
    
  




Multiplying the numerator and the denominator with sin {circumflex over (ω)} gives




  
    
      ∫
      
        ω
        min
      
      
        ω
        max
      
    
    ⁢
    
      
        
          S
          X
        
        ⁡
        
          (
          ω
          )
        
      
      ⁢
      
        d
        ω
      
    
  
  ≈
  
    
      ∫
      
        ω
        min
      
      
        ω
        max
      
    
    ⁢
    
      
        
          S
          X
        
        ⁡
        
          (
          ω
          )
        
      
      ⁢
      
        
          cos
          ⁢
          
            
          
          ⁢
          ω
          ⁢
          
            
          
          ⁢
          sin
          ⁢
          
            
          
          ⁢
          
            ω
            ^
          
        
        
          cos
          ⁢
          
            
          
          ⁢
          
            ω
            ^
          
          ⁢
          
            
          
          ⁢
          sin
          ⁢
          
            
          
          ⁢
          
            ω
            ^
          
        
      
      ⁢
      
        
          d
          ω
        
        .
      
    
  




Since (ωmax−ωmin) is small and ωmin≦ω≦ωmax, (ω+{circumflex over (ω)})≈2{circumflex over (ω)} and sin(ω−{circumflex over (ω)})≈(ω−{circumflex over (ω)}). Therefore,





  
    
      
        
          
            cos
            ⁢
            
              
            
            ⁢
            ω
            ⁢
            
              
            
            ⁢
            sin
            ⁢
            
              
            
            ⁢
            
              ω
              ^
            
          
          
            cos
            ⁢
            
              
            
            ⁢
            
              ω
              ^
            
            ⁢
            
              
            
            ⁢
            sin
            ⁢
            
              
            
            ⁢
            
              ω
              ^
            
          
        
        =
        
        ⁢
        
          
            
              0.5
              ⁢
              
                
              
              ⁢
              
                sin
                ⁡
                
                  (
                  
                    ω
                    +
                    
                      ω
                      ^
                    
                  
                  )
                
              
            
            -
            
              0.5
              ⁢
              
                
              
              ⁢
              
                sin
                ⁡
                
                  (
                  
                    ω
                    +
                    
                      ω
                      ^
                    
                  
                  )
                
              
            
          
          
            0.5
            ⁢
            
              
            
            ⁢
            sin
            ⁢
            
              
            
            ⁢
            2
            ⁢
            
              
            
            ⁢
            
              ω
              ^
            
          
        
      
    
  
  
    
      
        ≈
        
        ⁢
        
          1
          -
          
            
              
                ω
                +
                
                  ω
                  ^
                
              
              
                sin
                ⁢
                
                  
                
                ⁢
                2
                ⁢
                
                  ω
                  ^
                
              
            
            .
          
        
      
    
  




Substitution gives





  
    
      
        ∫
        
          ω
          min
        
        
          ω
          max
        
      
      ⁢
      
        
          
            S
            X
          
          ⁡
          
            (
            ω
            )
          
        
        ⁢
        
          d
          ω
        
      
    
    ≈
    
      
        ∫
        
          ω
          min
        
        
          ω
          max
        
      
      ⁢
      
        
          
            
              S
              X
            
            ⁡
            
              (
              ω
              )
            
          
          ⁡
          
            [
            
              1
              -
              
                
                  ω
                  -
                  
                    ω
                    ^
                  
                
                
                  sin
                  ⁢
                  
                    
                  
                  ⁢
                  2
                  ⁢
                  
                    
                  
                  ⁢
                  
                    ω
                    ^
                  
                
              
            
            ]
          
        
        ⁢
        
          d
          ω
        
      
    
  
  ,




implying that





  
    
      ∫
      
        ω
        min
      
      
        ω
        max
      
    
    ⁢
    
      
        
          S
          X
        
        ⁡
        
          (
          ω
          )
        
      
      ⁢
      
        
          (
          
            ω
            -
            
              ω
              ^
            
          
          )
        
        
          sin
          ⁢
          
            
          
          ⁢
          2
          ⁢
          
            
          
          ⁢
          
            ω
            ^
          
        
      
      ⁢
      
        d
        ω
      
    
  
  ≈
  0.




Since 2{circumflex over (ω)}≠0, solving for {circumflex over (ω)} yields





  
    ω
    ^
  
  ≈
  
    
      
        
          ∫
          
            ω
            min
          
          
            ω
            max
          
        
        ⁢
        
          
            
              S
              X
            
            ⁡
            
              (
              ω
              )
            
          
          ⁢
          ω
          ⁢
          
            d
            ω
          
        
      
      
        
          ∫
          
            ω
            min
          
          
            ω
            max
          
        
        ⁢
        
          
            
              S
              X
            
            ⁡
            
              (
              ω
              )
            
          
          ⁢
          
            d
            ω
          
        
      
    
    .
  




According to the above derivation, the single frequency estimation for a real-valued narrow-band lowpass signal using the eigendecomposition of the 3×3 element covariance matrix is approximately equal to the power-weighted mean frequency of the signal. Experimental results of this method indicate that it is relatively unbiased and provides accurate results for SNRs above 0 dB.
Magnitude Squared Coherence Function and Coherence Value
The maternal blood flow is a driving force to the placental-fetal system. Therefore, functional or structural changes that affect the maternal circulation should affect the flow through the uterine artery and the flow through spiral arteries. In addition to the inadequate secondary trophoblast invasion and impaired spiral arterial blood flow, in preeclamptic pregnancies there are significant lesions in the maternal uteroplacental blood flow resulting in increased placental insufficiency. It has been implied that the relationship of maternal-fetal blood flow may be impaired in preeclamptic pregnancies due to the changes in the maternal blood flow. Such changes will be most prominent in the frequency domain at the maternal heart rate, since the maternal blood flow drives the placental-fetal circulation system. Therefore, the magnitude-square coherence function (MSC) can be used to analyze differences between the uterine and umbilical blood velocity waveforms at the maternal heart rate. Here, the MSC function is always less than one because the fetal cardiovascular system is independently controlled by its autonomic nervous system.
Mathematically, the MSC function |γ(ω)|2 of the two blood velocity waveforms x(n) and y(n) at frequency ω is given by




  
    
      
        
        
          γ
          ⁡
          
            (
            ω
            )
          
        
        
      
      2
    
    =
    
      
        
          
          
            
              S
              xy
            
            ⁡
            
              (
              ω
              )
            
          
          
        
        2
      
      
        
          
            S
            xx
          
          ⁡
          
            (
            ω
            )
          
        
        ⁢
        
          
            S
            yy
          
          ⁡
          
            (
            ω
            )
          
        
      
    
  
  ,




where Sxy(ω) is the cross-spectrum, and Sxx(ω) and Syy(ω) are the power spectra of x(n) and y(n) respectively. The MSC is a function with values between 0 and 1 that indicates how well two waveforms correspond to each other in the frequency domain. If the two signals are linearly related, the coherence value can be equal to one. If the two signals are uncorrelated, the coherence between the two signals can be equal to zero and it can be less than one for non-linearly related signals. Similarly, the coherence value can be less than one whenever uncorrelated noise is present in the two signals.

The cross spectrum between x(n) and y(n) at frequency ω is estimated by taking the Fourier transform of the crosscorrelation function, and is given by




  
    
      S
      xy
    
    ⁡
    
      (
      ω
      )
    
  
  =
  
    
      ∑
      
        k
        =
        
          -
          
            (
            
              N
              -
              1
            
            )
          
        
      
      
        N
        -
        1
      
    
    ⁢
    
      
        
          
            r
            ^
          
          xy
        
        ⁡
        
          (
          k
          )
        
      
      ⁢
      
        e
        
          
            -
            j
          
          ⁢
          
            
          
          ⁢
          ω
          ⁢
          
            
          
          ⁢
          k
        
      
    
  




where {circumflex over (r)}xy(k) is the estimated value of the crosscorrelation function between x(n) and y(n) at lag k, N is the number of lag values used in estimating the cross spectrum. The crosscorrelation function between the two signals is estimated as





  
    
      
        r
        ^
      
      xy
    
    ⁡
    
      (
      k
      )
    
  
  =
  
    
      1
      M
    
    ⁢
    
      
        ∑
        
          n
          =
          1
        
        M
      
      ⁢
      
        
          x
          ⁡
          
            (
            n
            )
          
        
        ⁢
        
          y
          ⁡
          
            (
            
              n
              -
              k
            
            )
          
        
      
    
  




where M is the number of samples of the data sequences employed to make the estimate.

The power spectrum of x(n) is estimated as




  
    
      S
      xx
    
    ⁡
    
      (
      ω
      )
    
  
  =
  
    
      ∑
      
        k
        =
        
          -
          
            (
            
              N
              -
              1
            
            )
          
        
      
      
        N
        -
        1
      
    
    ⁢
    
      
        
          
            r
            ^
          
          xx
        
        ⁡
        
          (
          k
          )
        
      
      ⁢
      
        e
        
          
            -
            j
          
          ⁢
          
            
          
          ⁢
          ω
          ⁢
          
            
          
          ⁢
          k
        
      
    
  




where {circumflex over (r)}xx(k) is the autocorrelation function of x(n) given by





  
    
      
        r
        ^
      
      xx
    
    ⁡
    
      (
      k
      )
    
  
  =
  
    
      1
      M
    
    ⁢
    
      
        ∑
        
          n
          =
          1
        
        M
      
      ⁢
      
        
          x
          ⁡
          
            (
            n
            )
          
        
        ⁢
        
          
            x
            ⁡
            
              (
              
                n
                -
                k
              
              )
            
          
          .
        
      
    
  




Similarly, the power spectrum of y(n) is estimated as





  
    
      S
      yy
    
    ⁡
    
      (
      ω
      )
    
  
  =
  
    
      ∑
      
        k
        =
        
          -
          
            (
            
              N
              -
              1
            
            )
          
        
      
      
        N
        -
        1
      
    
    ⁢
    
      
        
          
            r
            ^
          
          yy
        
        ⁡
        
          (
          k
          )
        
      
      ⁢
      
        e
        
          
            -
            j
          
          ⁢
          
            
          
          ⁢
          ω
          ⁢
          
            
          
          ⁢
          k
        
      
    
  




where {circumflex over (r)}yy(k) is the autocorrelation function of y(n) given by





  
    
      
        r
        ^
      
      yy
    
    ⁡
    
      (
      k
      )
    
  
  =
  
    
      1
      M
    
    ⁢
    
      
        ∑
        
          n
          =
          1
        
        M
      
      ⁢
      
        
          y
          ⁡
          
            (
            n
            )
          
        
        ⁢
        
          
            y
            ⁡
            
              (
              
                n
                -
                k
              
              )
            
          
          .
        
      
    
  




Predetermined Threshold

An optimum threshold level that can be employed to identify pregnancies at risk for subsequent hypertension from normotensive pregnancies can be determined using an algorithm and is selected such that it maximizes the weighted sum of the sensitivity and the specificity of the algorithm. Sensitivity measures the proportion of the diseased subjects correctly identified by the test. Specificity of the test measures the proportion of the normal subjects correctly identified by the test. In order to select the threshold value that is applicable to all gestational ages, the gestational age dependent mean distribution is first removed from the estimated MSC γ(ωo)) values. We denote the mean removed values of γ(ωo) as {tilde over (γ)}(ωo).
Let Ho denote the null hypothesis that a given value of {tilde over (γ)}(ωo) is within the hypertensive region and the alternative hypothesis H1 be that {tilde over (γ)}(ωo) is within the normotensive region. Let pnorm({tilde over (γ)}(ωo)) and phyp({tilde over (γ)}(ωo)) be the probability density functions of {tilde over (γ)}(ωo) estimated for the normotensive and the hypertensive groups, respectively.
The conditional probability that {tilde over (γ)}(ωo) belongs to the null hypothesis is denoted by

p({tilde over (γ)}(ωo)|Ho)=phyp({tilde over (γ)}(ωo)).

Similarly, we define

p({tilde over (γ)}(ωo)|H1)=pnorm({tilde over (γ)}(ωo))

as the conditional probability the {tilde over (γ)}(ωo) belongs to the alternate hypothesis. Then, the threshold of detection {tilde over (γ)}(ωo) is chosen such that the cost function





  
    
      
        J
        =
        
        ⁢
        
          
            
              w
              l
            
            ⁢
            
              
                ∫
                
                  -
                  ∞
                
                
                  
                    
                      γ
                      ~
                    
                    th
                  
                  ⁡
                  
                    (
                    
                      ω
                      0
                    
                    )
                  
                
              
              ⁢
              
                
                  Phyp
                  ⁡
                  
                    (
                    
                      
                        γ
                        ~
                      
                      ⁡
                      
                        (
                        
                          ω
                          0
                        
                        )
                      
                    
                    )
                  
                
                ⁢
                
                  d
                  
                    γ
                    ~
                  
                
                ⁢
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
            
          
          +
        
      
    
  
  
    
      
        
        ⁢
        
          
            
              ω
              2
            
            ⁢
            
              
                ∫
                
                  
                    
                      γ
                      ~
                    
                    th
                  
                  ⁡
                  
                    (
                    
                      ω
                      0
                    
                    )
                  
                
                ∞
              
              ⁢
              
                
                  Pnorm
                  ⁡
                  
                    (
                    
                      
                        γ
                        ~
                      
                      ⁡
                      
                        (
                        
                          ω
                          0
                        
                        )
                      
                    
                    )
                  
                
                ⁢
                
                  d
                  
                    
                      γ
                      ~
                    
                    ⁡
                    
                      (
                      
                        ω
                        0
                      
                      )
                    
                  
                
              
            
          
          ,
        
      
    
  




is maximized, where ω1 and ω2 are positive weighting constants. Note that the first integral in the above expression represents the sensitivity and the second integral defines the specificity of the algorithm.

At the optimum threshold level {tilde over (γ)}th(ωo),




  
    
      
        ∂
        J
      
      
        ∂
        
          
            γ
            ~
          
          ⁡
          
            (
            
              ω
              0
            
            )
          
        
      
    
    =
    
      
        
          
            ω
            1
          
          ⁢
          
            Phyp
            ⁡
            
              (
              
                
                  γ
                  ~
                
                ⁢
                
                  
                
                ⁢
                
                  th
                  ⁡
                  
                    (
                    
                      ω
                      0
                    
                    )
                  
                
              
              )
            
          
        
        -
        
          
            ω
            2
          
          ⁢
          
            Pnorm
            ⁡
            
              (
              
                
                  γ
                  ~
                
                ⁢
                
                  
                
                ⁢
                
                  th
                  ⁡
                  
                    (
                    
                      ω
                      0
                    
                    )
                  
                
              
              )
            
          
        
      
      =
      0
    
  
  ,




implying that

ω1Phyp({tilde over (γ)}th(ω0))=ω2Pnorm({tilde over (γ)}th(ω0)).

The solution to the above equation provides the optimum threshold level {tilde over (γ)}th(ωo).

Preliminary Studies

In preliminary studies it was assumed that there was no significant variation in the environmental, fetal-umbilical, and maternal-uterine waveform characteristics within short intervals of time. The blood velocity waveforms were recorded sequentially within a 5-minute interval. A delay in one of the signals changes only the phase of the cross-spectrum of the two signals. Since the MSC function contains no information about the relative phase between the two signals, the estimated MSC function from the sequentially acquired data will not show much difference from the estimates obtained from simultaneously acquired data. This assumes that the joint statistics of the two signals do not change during the time frame in which the data were collected.
A total of 44 maternal-fetal data sets were collected during the 10th through the 22nd week of gestation, of which 34 resulted in normal pregnancy outcome. The rest of the patients (10 data sets) subsequently developed PIH. Transvaginal (5 MHz) and transabdominial (3.75 MHz) Doppler recording methods were used at 10-13 and 14-22 weeks of gestation, respectively. All Doppler studies were performed with the women in the semi-recumbent position and during fetal apnea. The umbilical arterial velocity waveforms were obtained from the free floating loop of the umbilical cord. The angles of incidence of the ultrasound beams in the measurements were less than 20°. The sampling frequency of the blood velocity waveforms was 93.75 samples/sec and each waveform was recorded for approximately 20 seconds.
In estimating the coherence values at the maternal heart rates for each pair of waveforms, first, the mean values of the maternal heart rate was estimated for each waveform as the frequency corresponding to the peak of the amplitude spectrum. The amplitude spectrum was estimated as the absolute value of the Fourier transform estimated over the entire length of the blood velocity waveform. The power spectra estimated over 16 or more 256-sample data segments. Each data segment was overlapped by 75% by the adjacent segment. Then, the MSC values were estimated using the equation: |γ(ω)|2=(|Sxy(ω)|2)/(Sxx(ω)Syy(ω)).
The distribution of γ(ωo) is shown in FIG. 5. The open circles denote results for normotenisive pregnancies and asterisks represent the corresponding results for hypertensive pregnancies. The average distribution of γ(ωo) estimated as a function of the gestational age for normotenisive pregnancies is shown by the solid line in FIG. 5. This curve was obtained using a first order polynomial fit, and satisfies the relationship

{tilde over (γ)}(ω0, θ)=−0.0092θ+0.2885,

where γ(ωo, θ) is the mean value of γ(ωo) at the gestational age θ. FIG. 5 demonstrates that γ(ωo) values estimated for pregnant women who subsequently developed PIH are in general well below the average for gestational age. Therefore, the average distribution of γ(ωo) at maternal heart rate can be used as a marker for identifying complicated pregnancies.

Selection of a threshold function to detect pregnancies that may later develop preeclampsia is done using

ω1Phyp({tilde over (γ)}th(ω0))=ω2Pnorm({tilde over (γ)}th(ω0)).

First, {tilde over (γ)}(ωo) values were estimated by subtracting the gestational age dependent mean distribution from γ(ωo) values. It was assumed that {tilde over (γ)}(ωo) in both the normotensive and hypertensive pregnancies belong to truncated Gaussian distributions in the interval [{tilde over (γ)}min(ωo), {tilde over (γ)}max(ωo)]. That is,





  
    Phyp
    ⁡
    
      (
      
        
          γ
          ~
        
        ⁡
        
          (
          
            ω
            0
          
          )
        
      
      )
    
  
  =
  
    {
    
      
        
          
            
              
                1
                
                  C
                  hyp
                
              
              ⁢
              
                N
                ⁡
                
                  (
                  
                    
                      μ
                      hyp
                    
                    ,
                    
                      σ
                      hyp
                      2
                    
                  
                  )
                
              
            
            ;
          
        
        
          
            
              
                
                  γ
                  ~
                
                min
              
              ⁡
              
                (
                
                  ω
                  0
                
                )
              
            
            ≤
            
              
                γ
                ~
              
              ⁡
              
                (
                
                  ω
                  0
                
                )
              
            
            ≤
            
              
                
                  γ
                  ~
                
                max
              
              ⁡
              
                (
                
                  ω
                  0
                
                )
              
            
          
        
      
      
        
          
            0
            ;
          
        
        
          
            
              
                
                  
                    γ
                    ~
                  
                  min
                
                ⁢
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
              >
              
                
                  γ
                  ~
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
              >
              
                
                  
                    γ
                    ~
                  
                  max
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
            
            ,
          
        
      
    
  




where (μhyp, σhyp2) represents a Gaussian density function with mean μhyp and variance σhyp2, and





  
    C
    hyp
  
  =
  
    
      1
      
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          
            σ
            hyp
            2
          
        
      
    
    ⁢
    
      
        ∫
        
          
            
              t
              ~
            
            mni
          
          ⁡
          
            (
            
              ω
              0
            
            )
          
        
        
          
            t
            ~
          
          
            max
            ⁡
            
              (
              
                ω
                0
              
              )
            
          
        
      
      ⁢
      
        exp
        ⁢
        
          {
          
            -
            
              
                
                  (
                  
                    
                      
                        γ
                        ~
                      
                      ⁡
                      
                        (
                        
                          ω
                          0
                        
                        )
                      
                    
                    -
                    
                      μ
                      hyp
                    
                  
                  )
                
                2
              
              
                2
                ⁢
                
                  σ
                  hyp
                  2
                
              
            
          
          }
        
        ⁢
        
          
        
        ⁢
        
          
            d
            
              
                γ
                ~
              
              ⁡
              
                (
                
                  ω
                  0
                
                )
              
            
          
          .
        
      
    
  




Similarly,





  
    Pnorm
    ⁡
    
      (
      
        
          γ
          ^
        
        ⁡
        
          (
          
            ω
            0
          
          )
        
      
      )
    
  
  =
  
    {
    
      
        
          
            
              
                
                  1
                  
                    C
                    norm
                  
                
                ⁢
                
                  N
                  ⁡
                  
                    (
                    
                      
                        μ
                        norm
                      
                      ,
                      
                        σ
                        norm
                        2
                      
                    
                    )
                  
                
              
              ;
            
          
          
            
              
                
                  
                    γ
                    ~
                  
                  min
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
              ≤
              
                
                  γ
                  ~
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
              ≤
              
                
                  
                    γ
                    ~
                  
                  max
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
            
          
        
        
          
            
              
                0
                ;
              
              ⁢
              
                
              
            
          
          
            
              
                
                  
                    γ
                    ~
                  
                  min
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
              >
              
                
                  γ
                  ~
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
              >
              
                
                  
                    γ
                    ~
                  
                  max
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
            
          
        
      
      ,
    
  




where Cnorm, μnorm and σnorm2 are defined as before for normotensive pregnancies. Substituting the above equations and solving for {tilde over (γ)}th(ωo) results in




  
    
      
        
          
            γ
            th
          
          ~
        
        ⁡
        
          (
          
            ω
            0
          
          )
        
      
      =
      
        
          
            
              σ
              hyp
              2
            
            ⁢
            
              
                
                  γ
                  ~
                
                norm
              
              ⁡
              
                (
                
                  ω
                  0
                
                )
              
            
          
          -
          
            
              σ
              norm
              2
            
            ⁢
            
              
                
                  γ
                  ~
                
                hyp
              
              ⁡
              
                (
                
                  ω
                  0
                
                )
              
            
          
          -
          
            K
          
        
        
          
            σ
            hyp
            2
          
          -
          
            σ
            norm
            2
          
        
      
    
    ,
    
      

    
    ⁢
    where
  




  
    K
    =
    
      
        
          (
          
            
              
                σ
                hyp
                2
              
              ⁢
              
                
                  
                    γ
                    ~
                  
                  norm
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
            
            -
            
              
                σ
                norm
                2
              
              ⁢
              
                
                  
                    γ
                    ~
                  
                  hyp
                
                ⁡
                
                  (
                  
                    ω
                    0
                  
                  )
                
              
            
          
          )
        
        2
      
      -
      
        (
        
          
            σ
            hyp
            2
          
          -
          
            σ
            norm
            2
          
        
        )
      
    
  




  
    
      
    
    ⁢
    
      
        (
        
          
            
              σ
              hyp
              2
            
            ⁢
            
              
                
                  γ
                  ~
                
                norm
                2
              
              ⁡
              
                (
                
                  ω
                  0
                
                )
              
            
          
          -
          
            
              σ
              norm
              2
            
            ⁢
            
              
                
                  γ
                  ~
                
                hyp
                2
              
              ⁡
              
                (
                
                  ω
                  0
                
                )
              
            
          
          -
          
            2
            ⁢
            
              σ
              hyp
              2
            
            ⁢
            
              σ
              norm
              2
            
            ⁢
            
              ln
              ⁡
              
                [
                
                  
                    
                      σ
                      hyp
                    
                    
                      σ
                      norm
                    
                  
                  ⁢
                  
                    
                      ω
                      2
                    
                    
                      ω
                      1
                    
                  
                  ⁢
                  
                    
                      C
                      hyp
                    
                    
                      C
                      norm
                    
                  
                
                ]
              
            
          
        
        )
      
      .
    
  



The mean-removed MSC function at the maternal heart rate ranged from {tilde over (γ)}min(ωo)=−0.15 to {tilde over (γ)}max(ωo)=0.2 for the data used experimentally. The variance and mean values of {tilde over (γ)}(ωo) were σnorm2=0.0059 and μnorm=0, respectively, for normotensive pregnancies and σhyp2=0.0046 and μhyp=−0.00764, respectively, for hypertensive pregnancies. The results may change for larger test groups.
The variation of the threshold level {tilde over (γ)}th(ωo) estimated using the equations above with various choices of the two weighting factors ω1 and ω2 is shown in FIG. 6. With the increase in the ω1, the weight for sensitivity, the threshold level increases whereas with the increase in ω2, the weight factor for specificity, the threshold level decreases. The distribution of the sensitivity Psen and specificity Pspe are shown in FIG. 7 for the above ranges of ω1 and ω2. The distribution of mean-removed coherence values and the threshold level for ω1=5 and ω2=4 is shown in FIG. 8. In the preliminary testing, the optimum threshold level was −0.0173. The theoretically calculated sensitivity was 80.03% and the corresponding specificity was 58.11%. The experimentally evaluated sensitivity and specificity of the test were 90% and 52.94%, respectively. With 22.73% prevalence of the condition in the acquired data, the positive predictive value of the test was 36% and the negative predictive value was 94.74%
Although the results stated above were obtained for the data on which the algorithm was trained to select the threshold, in practice, the algorithm should be tested on data that were not employed to train the system. Another round of testing was performed on the data sets to see how well the method works when the training data was not the testing data. Of the original 44 data sets, 34 (27 normal and 7 PIH) were selected to serve as the training data sets and the remaining 10 (7 normal and 3 PIH) were used for testing the method. The mean value and variance of the threshold value, sensitivity and specificity associated with all possible combinations of training and test data sets for ω1=5 and ω2=4, and 80% theoretical sensitivity are displayed in Table I in FIG. 9. The results were close to 80% sensitivity for ω1=5 and ω2=4. With 80% sensitivity, the results were approximately similar. From this, it is expected that the method will perform well in practical situations.
Comparing these results with other available methods that use Doppler ultrasound data reveals the other methods to be unreliable. With respect to using the uterine-arterial pulsatility index (PI), there is no significant variation in the resistance to the uterine arterial blood flow at 10-22 weeks in hypertensive pregnancies from normotensive pregnancies. Also, with respect to the umbilical-arterial blood velocity variability analysis, there was no significant difference in the variance of the blood velocity variability in normotensive and hypertensive pregnancies. Therefore, it appears that neither PI analysis nor blood velocity variability are reliable indicators of subsequent development of preeclampsia.
The main advantages of this coherence analysis for the prediction of preeclampsia are: 1) it can be used to identify women at risk for hypertensive disorders at an earlier stage of pregnancy than currently possible, 2) it is non-invasive and simple to perform, and 3) it is highly sensitive.
The embodiments of the method 100 and apparatus 200 of the present invention may be used to determine other bodily parameters such as, for example and not by way of limitation, maternal and fetal heart rate, HRV and electrocardiogram (EKG). Additionally, the method 100 and apparatus of the present invention may be applied in other circumstances where determining a coherence value relating to two waveforms using a magnitude-squared coherence function and is therefore not limited to predicting preeclampsia. For example, milder forms of the disease such as pregnancy induced hypertension or generally maternal hypertension during pregnancy may also be predicted according to embodiments of the invention described herein.
It is to be understood that the above-referenced arrangements are illustrative of the application for the principles of the present invention. Numerous modifications and alternative arrangements can be devised without departing from the spirit and scope of the present invention while the present invention has been shown in the drawings and described above in connection with the exemplary embodiments of the invention. It will be apparent to those of ordinary skill in the art that numerous modifications can be made without departing from the principles and concepts of the invention as set forth in the claims.

## Claims

1. A system for predicting maternal hypertension during pregnancy, comprising:
a processor for calculating a coherence value using a magnitude-squared coherence function between a maternal blood velocity waveform and a fetal blood velocity waveform, the magnitude-squared coherence function defined by:


|γ(ω)|2=(|Sxy(ω)|2)/(Sxx((ω)Syy((ω)),



wherein |Sxy(ω)|2 is a cross spectrum of the maternal and fetal blood velocity waveforms and Sxx(ω) and Syy(ω) represent a power spectra of the maternal and fetal blood velocity waveforms;
a maternal sensor in communication with the processor for measuring the maternal blood velocity waveform;
a fetal sensor in communication with the processor for measuring the fetal blood velocity waveform, wherein the maternal sensor and the fetal sensor are the same sensor; and
a comparator configured to compare the coherence value to a predetermined threshold value to identify pregnancies at risk for subsequent hypertension.

2. The system of claim 1, wherein the maternal and fetal sensors are in communication with at least one Doppler ultrasound apparatus.

3. The system of claim 1, wherein the maternal and fetal sensors are Doppler ultrasound sensors.

4. The system of claim 1, wherein the predetermined threshold value can be calculated according to the equation:


ω1Phyp({tilde over (γ)}th(ω0))=ω2Pnorm({tilde over (γ)}th(ω0)).

5. The system of claim 1, wherein the magnitude-squared coherence function can further be used to determine a coherence value between a maternal heart rate and a fetal heart rate.

6. The system of claim 1, wherein the magnitude-squared coherence function can further be used to determine a coherence value between a maternal electrocardiogram and a fetal electrocardiogram.

7. A system for predicting maternal hypertension during pregnancy, comprising:
a processor for calculating a coherence value using a magnitude-squared coherence function between a maternal blood velocity waveform and a fetal blood velocity waveform, the magnitude-squared coherence function defined by:


|γ(ω)|2=(|Sxy(ω)|2)/(Sxx(ω)Syy(ω)),



wherein |Sxy(ω)|2 is a cross spectrum of the maternal and fetal blood velocity waveforms and Sxx(ω) and Syy(ω) represent a power spectra of the maternal and fetal blood velocity waveforms;
a maternal sensor in communication with the processor for measuring the maternal blood velocity waveform; and
a fetal sensor in communication with the processor for measuring the fetal blood velocity waveform;
a comparator configured to compare the coherence value to a predetermined threshold value to identify pregnancies at risk for subsequent hypertension, wherein the predetermined threshold value can be calculated according to the equation:


ω1Phyp({tilde over (γ)}th(ω0))=ω2Pnorm({tilde over (γ)}th(ω0)).

8. The system of claim 7, wherein the maternal and fetal sensors are in communication with at least one Doppler ultrasound apparatus.

9. The system of claim 7, wherein the magnitude-squared coherence function can further be used to determine a coherence value between a maternal heart rate and a fetal heart rate.

10. The system of claim 7, wherein the magnitude-squared coherence function can further be used to determine a coherence value between a maternal electrocardiogram and a fetal electrocardiogram.