The fact that fPCG acquisition is highly susceptible to noises paves the way to innovative hardware oriented methods for a less noisy acquisition. In this study, fPCG are recorded using vibration sensors (cost $1 each) embedded in high definition 3D printed plastic harnesses. Each harness holds a ceramic piezo vibration sensor (35 mm diameter) on the maternal abdomen with rubber made cushion to minimize the shear noise. The 3D printed harness is designed with precise parameters that rigidly mount the piezo sensor. The sketch in Figure 1 shows the setup of the sensors, where each sensor picks fPCG signals through a coaxial cable having very high insulating resistance. Power lab data acquisition system by ADinstrument (www.adinstruments.com) was used to record the abdominal phonograms at a sampling frequency of f_s = 1, 000Hz. In order to validate the extracted fHR, 14 simultaneous CTG recordings were collected by Monica wireless CTG AN24 (www.monicahealthcare.com) and, 6 by Phillips Avalon FM300 CTG (www.usa.philips.com) devices for 20 subjects in total.

The prototype was tested on 20 pregnant mothers in Al-Ain Hospital, United Arab Emirates. Al Ain District Ethics Committee approved this study (ref: AAHEC-09-14-013) and informed signed consent forms were obtained from all pregnant volunteers. Figure 2 shows an example of recorded raw fPCG signal over time (seconds), which contains clear maternal breathing patterns (exhaling and inhaling) and fHS, as well, modulated, however, by the maternal breathing patterns, setting the challenge for efficient source separation. The total weight of the prototype (four channel sensors with harness) is 200 g and it was comfortable to wear which was confirmed by end of measurement survey questionnaire.

The special problem of extracting sources from a set of linear mixtures without any knowledge of the mixing channel has been widely studied by assuming different sources statistical conditions (Parra and Sajda, 2003; Naik and Wang, 2014). In its simplest form, it can be expressed as the problem of identifying the factorization of P-dimensional observations x into a mixing channel A and Q-dimensional sources s.

Let x^o[n] be a multichannel discrete-time signal resulting from the synchronous sampling of P continuous-time signals (channels), where n denotes discrete time. The multichannel signal x^o[n] corresponds to the linear combination of Q unknown source signals s[n] according to

where A is the P × Q mixing matrix. The observation signals in x^o[n] are given in row format (as well as the source signals in s[n]), such that xio[n] denotes the ith channel in x^o[n]. The problem we deal with is the recovery of the sources s[n] from the observations x^o[n], without any knowledge of the mixing matrix A.

In our particular problem, the observation signals x^o[n] correspond to the synchronous reading of P = 4 vibration sensors. Given the sensor arrangement (Figure 1), placed 16 cm apart from each other, a distance much smaller than the minimum wavelength (λ = 3m in liquid for the specific f_s), each source signal s_i[n] arrives to every sensor with a different intensity, usually based on how close the sensor is to its respective physical source. We assume s[n] to be related to the mother's breathing, mother's heart, fetal heart, and noise, hence Q ≲ P. In consequence, the mixing matrix A is assumed to be of rank P, but not necessarily well-conditioned.

Let e be a diagonal control matrix to disable bad channels of M × P multiplied by x^o[n], where M ≤ 4 is the number of strongly correlated channels, i.e.,

We define the observation cross-covariance matrix as

where for such operation a N-point window is involved, that is, the ith row and jth column of matrix X[n] is obtained as follows

Based on the mixing scenario (1), we can express matrix (3) in terms of the mixing matrix A as

where S[n] is the source cross-covariance (unknown) matrix.

We are interested in a demixing matrix W^T such that W^TA = I or A^TW = I. This would separate the mixed signals as,

The expression in Equation (5) aims to provide general cross statistics between shifted versions of the collected sensors readings. Using Equation (3), we can refer to X[0] as the zeroth lag covariance matrix of the sensors recordings, while X[1] as the 1st lag cross-covariance matrix and X[n] as the nth lag cross-covariance matrix. These matrices will be used to devise an algorithm to recover the original components of a specific set of mixtures of data. In particular, starting from the zeroth lag it can be written

In Equation (7), the sensors recordings zeroth lag covariance matrix X[0] is written in terms of the sources zeroth lag covariance matrix S[0] and the mixing matrix A. By multiplying each right side of Equation (7) by W it would result in

This process is repeated for n = 1, 2…k as shown below

Blind source separation utilizing only the first lag covariance matrix X[1] has been shown and discussed in Parra and Sajda (2003), Weinstein et al. (1993), Parra and Spence (2000). Belouchrani et al. (1997) introduced a source separation technique utilizing the time coherence of the source signals. This relies only on stationary second-order statistics that are based on a joint diagonalization of a set of covariance matrices. In our mixing scenario, we assume that the sources are non-white, stationary and decorrelated. The aim here is to provide a simplified understanding of the stationary assumption by studying sums of multi-lag cross covariance up to kth-lag sample in a form of eigenvalue decomposition problem. As Equations (9–11) share the same mixing matrix A, let the kth sample be the discrete time lag over a very short period of time; the sum of the nth lag covariance matrices up to k is given by

We multiply each right side of Equation (12) by a desired demixing matrix W to result in the following

and the original mixing matrix A can be written as,

The goal here is to find the minimum sum to kth lag covariance matrix that constructs an eigenvalue decomposition problem to best approximate the demixing matrix W^T. With the substitution of the original mixing matrix A derived in Equation (14) in the zeroth covariance form derived in Equation (8) we get

The above expression can be further written as

where Λ=(∑n=1kS[n])-1S[0] is a diagonal matrix and S[0], S[1] …S[k] have the same diagonalization property over a short period of time.

For simplicity, we want to produces a diagonal matrix D=Λ of generalized eigenvalues and a full matrix V=W whose columns are the corresponding eigenvectors where A=X[0] and B in this case is nothing but ∑n=1kX[n] so that the problem is in the form of an eigenvalue decomposition problem:

The idea in Equation (16) is that the full matrix W is nothing but the eigenvectors matrix constructed not only from the covariance matrix X[0], but also from the sum of the minimum kth lag cross covariance matrix ∑n=1kX[n]. Because the sources are non-white, this method is described as using second order statistics in the form of cross-covariance for different time lags. In this way, the issue is transfered to the correlation maximization problem.

As previously discussed and shown in Equation (2), e matrix disables any bad channels or dismiss any operating window if they are not sufficiently correlated. A linear set of channels tend to have high degree of correlation because the independent components are available in every channel but with different weights. Let C(i, j) be the zeroth lag of a normalized covariance matrix between the ith and the jth channels; the correlation coefficient R(i, j) is, then, defined as

where a matrix R written in a form of table. Table 1 gives the degree of correlation between all the channels. Performing such practice would give us valuable information on the channels similarities. The goal here is to extract a binary control matrix e to be multiplied by the raw collected signals x^o[n] and disable any noisy (less correlated channels).

The upper triangle in Table 1 is converted to binary, based on a threshold value (initially set as 0.3); hence, any |R(i, j)| < 0.3 is set to 0 while those ≥ 0.3 are set to 1. Such practice would disable channels that are less correlated and would improve the blind source separation results. The control matrix e can be extracted from the logic OR gate of the R(i, j) values of the current channel with the remaining channels. The control matrix e is constructed as

where if any e(i) value is set to zero the raw containing it is removed. This method can disable noisy channels and dismiss all of the sensors if all of them are noisy. Table 1 summarizes the proposed process.

The algorithm described in this section gives a close approximation of the original mixing matrix. In other words, the estimated demixing matrix W^T is able to recover the independent components, which are the maternal breathing, mHS and the desired fHS, respectively.

Steps through application:

An actual maternal breathing pattern (recorded through vibration sensor), fetal ECG (fECG), and maternal ECG (mECG) signals were recorded to create synthetic channels. Because of the nature of the abdominal fading coefficients, the abdominal channel between the sources and the sensors are in continuous change due to mainly lung inhalation and exhalation and accordingly, its effect on the medium. In an attempt to mimic this scenario, 50 iterations of different positive normally distributed pseudorandom mixing matrices are used to mix the original sources as seen in Figure 3. The constructed mixing signals were then used to test the derived algorithm and to estimate the demixing matrix at each iteration with different k-values.

The source separation process for abdominal phonogram signals is described in Figure 4. The proposed method was tested on abdominal phonogram signals from 20 pregnant women; their anthropometric characteristics are listed in Table 2. Initially, the electric hum of 50 Hz was removed, using a notch filter with quality factor of 60. Then, the proposed source separation technique (Section 2.2) was applied to the filtered data and the mHS and fHS were extracted. These were further denoised using a wavelet transform-based stationary-non-stationary filter (WTST-NST) (Hadjileontiadis and Panas, 1997). The WTST-NST filter basically utilizes the multi-resolution wavelet decomposition with dynamic threshold values at each level to separate non-stationary signal components (noise) from stationary signal components (heart sounds). The peaks of the heart sounds were identified by finding the local maxima in a pre-specified peak-to-peak time distance of 0.3 s (assuming maximum fHR of 180 bpm). In this window, the data samples were compared to their neighboring values, where the peak of maximum energy was highlighted. If any element of data was larger than its neighbors, the element was considered as a local heart sound peak. The separated denoised fHS were then used to estimate the fHR for comparison with the same from the CTG machine.

Figures 5A–D show a comparison between the approximated mixing vectors (plotted in red) using the expression derived in Equation (16) and the original mixing vectors (plotted in blue), for a full rank mixing matrix of 3 and for k= 1, 6, 50, and 100, respectively. The main purpose here is to overlap the vectors as much as possible or, in other words, minimize the 2-norm between the approximated mixing matrix and the original one. In addition, Figure 6 shows the abs 2-norm difference between the approximated mixing matrix W and the original one to indicate the suitable sum of the kth lag. The least 2-norm range occurs in the range from 6th lag (equivalent to 0.006 s) to a max of 60th lag (equivalent to 0.06 s) for the given f_s (1,000 Hz). The 6th lag sum gives the best approximation results, while minimizing the computational process. The 0.006–0.06 s gives the indication that in practical scenarios the time lagged cross-covariance matrices can not share across all the time lags the same A, because practically we can partially assume that the signal is stationary over a short period of time just as the 0.06 s. The sampling frequency here will not affect the ideology as generally we define its sufficient value to make sure that we detect the maximum needed frequency in our collected data. To be more specific, the advantage of using a 1,000 Hz sampling frequency appears when it comes to record the peaks of the detected fHS, a less sampling frequency would still collect the desired frequency component but would miss a peak point in fHS. It is very critical to obtain clear set of signals. The ECG data set that was collected to construct our synthetic model was recorded using AD instrument power lab four channel data acquisition system at a 1,000 Hz sampling rate, which produced smooth heart sounds. In industry, an ECG sampling interval of 1 ms is recommended to get accurate time domain measures. Accordingly, our PCG data were recorded at a 1,000 Hz, to acquire enough time domain measure. It is worth mentioning that a less sampling rate would result in a less number of time lags, but the time sample which happens to be at sample number 60 in Figure 6 would be missed in a less number of time lags. So, a less number of time lags due to a reduced sampling frequency would affect the resolution of Figure 6 and would not give us a clear picture of our kth-lag range.

Figure 7 shows the results of the proposed source separation approach when applied to the simulated data, using the 6th lag. When comparing the initial sources used for mixing (Figure 3) with the ones estimated in Figure 7, a clear source separation can be noticed, since the maternal breathing pattern, mHS and the fHS are correctly separated. In approximately all of the 50 iterations (program runs), the suitable k-value range was concluded the same.

The raw data shown in Figure 2 of subject ID 145 were separated intro their source components, as seen in Figure 8, i.e., the maternal breathing, mHS, and fHS. The extracted fetal heart rates were then compared with the same obtained from Monica wireless CTG AN24 over the same time span, as shown in Figure 9. In general, the PCG-based fHR shows good similarity with the one obtained by Monica wireless CTG AN24, except from some sudden changes noticed in baseline (e.g., around 12th second in Figure 9). This is however expected, as Monica CTG AN24 uses a 3.75 s averaging window; hence, the original variation in fHR might have been lost due to this averaging process.

Furthermore, fHR extracted from fPCG were also compared with the same from the Phillips Avalon FM300 CTG device. In Figure 10, the fHR has been obtained from the same fHS panel but extended to 60 s (inset panel) to produce the same CTG portion highlighted within a square. The 60 s window is used, since the Phillips Avalon FM300 CTG device has a resolution of 60 s data/cm as the drawn 1 cm red square shows in Figure 10. The corresponding fPCG-based fHR (in the embedded subfigure of Figure 10 connected with the red square) shows good similarity with the same from the Philips device, as seen in the drawing pattern of the fHR, where the variability pattern is approximately similar. Acceleration and deceleration of the fHR are important signs for fetal well being. Figure 11 shows the denoised fHS detected of subject ID 195. Using fPCG-based monitoring, clear signs of acceleration and deceleration were noticed. The fPCG-based fHR were plotted and compared with the same from its Monica wireless CTG AN24 version, exhibiting a good matching in the acceleration/deceleration periods.

From an overall comparative perspective, Figure 12 shows the generated Bland-Altman plot of the mean fHR values (fPCG vs. CTG) in all 20 subjects, where Bland–Altman mean = −0.21 BPM and ±2 SD = ±3 BPM are identified. Moreover, Figure 13 depicts the Spearman correlation between the fHR from fPCG and fECG, respectively, exhibiting a ρ value of 0.95 (significance level of p < 0.001). Apparently, these results confirm that the fPCG-based fHR is comparable to those obtained by Monica wireless CTG AN24 and by Phillips Avalon FM300 CTG device. Such sensors setup and signal processing algorithm is of low cost and can be run on any PC and it is safe for long term monitoring. There is be no need for skilled operator to be able to operate the system, as the sensors harness can easily be placed on the abdomen, which could open up the possibility for home monitoring of fetal well-being.

The acceptable matches of mean fHR values and fHR baseline between the PCG and CTG data, verified by Bland–Altman and the Spearman correlation plots, support the application of PCG-based obstetric medical assessment rules on all the tested subjects. However, both PCG and CTG scan-based obstetric medical assessment suggested that the tested subjects had no urgent clinical issues in this study.

The novelty of this study includes the proof of concept that multi-channel low-cost vibration sensors harnessed in high definition 3D printed casings were able to reliably pick the fetal heart sounds. The control matrix e extracted from the OR logic of the zeroth lag covariance has been used to disable bad channels before applying any processing. Our proposed signal processing technique introduces the concept of summing multi-lag cross covariance matrices in an eigen value decomposition form. The multi-lag cross covariance sum is up to kth-lag sample. Depending on the final kth-lag cross covariance sum, the approximated demixing matrix performance changes, where we showed that multi-lag cross covariance sum up to kth-lag = 6–60 is reliable to separate fetal heart sounds from abdominal phonograms. In this acceptable range, choosing k = 6 gives the best optimization in terms of extraction performance and computational power as less sums will be computed. Besides, the multi-resolution wavelet decomposition with dynamic threshold values at each level to separate non-stationary signal components (noise) from stationary signal components (heart sounds) proved to be effective in denoising extracted noisy fetal heart sounds. The proposed method does not assume independent signals as this is a hardly met condition on abdominal surface due to different sources of noises. However, it looks at the best way to maximally decorrelate the sources. Due to the fact that it involves decomposition of the sources using extracted eigenvectors, puts it as a good candidate for fast processing specially in real-time processes. The length of the extracting vector equals the number of the maximally decorrelated channels, which imposes no storage limitation. The algorithm does not need any complex operations, while it only involves multiplications and additions which would not impose any heavy computation on the CPU. These advantages would qualify the algorithm to be run on any smart phone for domestic use.

The setup of our data acquisition system helped our recording environment to be as close to non-noisy as possible. The high definition 3D printed sensors holders harness robustly the vibration sensors on the abdomen that isolate them from possible unneeded external BSS context; such as, any speech signal. This along with our e control matrix narrow down the channels into maximally correlated ones. This helps our four sensors or less (depending on the e matrix) to receive the desired BSS context but with different linear combination due to the medium.

Possible unneeded abdominal BSS context; such as, muscle noise and baby kicks are not a big concern. These events don't happen all the time and can be tracked if a loss of fHS is witnessed. This is similar to the methodology used in the CTG scans.

However, further validations of PCG-based fHR on a variety of non-reassuring patterns of fHR, as seen in clinical scenarios, such as early decelerations, late decelerations, prolonged deceleration, recurrent, sinusoidal fHR, and variable decelerations are needed before the proposed technique could be used for a potential low cost antenatal care system. It is important to mention that we did not investigate how maternal weight such abdominal adiposity could influence the phonocardiogram signal quality in this study.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer LE and handling Editor declared their shared affiliation.