Fusion signal processing for maternal uterine activity detection

Abstract

The computer-implemented method includes providing, by at least one computer processor, a plurality of signal channels, the plurality of signal channels including a plurality of electrical uterine monitoring signal channels and a plurality of acoustic uterine monitoring signal channels; determining, by the at least one computer processor, a plurality of channel weights, each of the channel weights corresponding to a particular one of the signal channels; and generating, by the at least one computer processor, a combined uterine monitoring signal channel by calculating a weighted average of the signal channels based on the channel weight for each of the signal channels.

Description

[Technical Field]

[0001] Cross-Reference to Related Applications This is an international (PCT) patent application related to and claims the benefit of co-filed and co-pending U.S. Provisional Patent Application No. 62 / 970,585, filed February 5, 2020, and entitled "Fused Signal Processing for Maternal Uterine Activity Detection," the contents of which are incorporated herein by reference in their entirety.

[0002] Technical field to which the invention belongs The present invention relates generally to maternal monitoring, and more particularly to the analysis of sensed biopotential and / or acoustic data to generate calculated representations of uterine activity, such as uterine contractions. [Background technology]

[0003] background Uterine contractions are a temporary process during which the uterine muscles shorten and the space between the muscle cells decreases. These structural changes in the muscles cause an increase in uterine cavity pressure to allow the fetus to be pushed down toward labor. During uterine contractions, the structure of uterine muscle cells (i.e., uterine cells) changes, causing the uterine walls to thicken. Figure 1A is an illustration of an atonic uterus with relaxed uterine muscle walls, and Figure 1B is an illustration of an atonic uterus with relaxed uterine muscle walls. Figure 1B is an illustration of a uterine contraction in which the uterine muscle walls contract to press the fetus against the cervix.

[0004] Uterine contractions are monitored to assess the progress of labor. Typically, labor progress is monitored using two sensors: a tocodynamometer, a strain-gauge-based sensor placed on the pregnant woman's abdomen, and an ultrasound transducer, also placed on the abdomen. Uterine contractions can be identified by analyzing the tocodynamometer signal to obtain a tocograph (TOCO), and fetal and maternal heart rates and fetal movements can be detected by analyzing the ultrasound transducer signal. However, these sensors can be uncomfortable to wear and can produce unreliable data when worn by obese pregnant women. [Prior art documents] [Patent documents]

[0005] [Patent Document 1] U.S. Patent No. 9,713,430 [Patent Document 2] U.S. Patent No. 9,392,954 [Non-patent literature]

[0006] [Non-Patent Document 1] Hyvarinen et al., “Independent Component Analysis: Algorithms and applications,” Neural Networks 13(4-5):411-430 (2000) [Non-patent document 2] Zong et al., “A QT Interval Detection Algorithm Based On ECG Curve Length Transform,” Computers In Cardiology 33:377-380 (October 2006) Summary of the Invention [Means for solving the problem]

[0007] overview In some embodiments, the present invention provides a tangibly programmed computer system including at least the following components: a non-transitory memory that electronically stores computer-executable program code; and at least one computer processor that, when executing the program code, becomes a tangibly programmed computing processor configured to perform at least the following operations: receive a plurality of biopotential signals collected at a plurality of locations on a pregnant mother's abdomen, detect R-wave peaks in the biopotential signals, extract maternal electrocardiogram ("ECG") signals from the biopotential signals, determine R-wave amplitudes of the maternal ECG signals, create an R-wave amplitude signal for each of the maternal ECG signals, calculate an average of all the R-wave amplitude signals, and normalize the average to generate an electrical uterine monitoring ("EUM") signal. In some embodiments, the operations also include identifying at least one uterine contraction based on a corresponding at least one peak in the EUM signal.

[0008] In some embodiments, the present invention provides a method that includes receiving multiple biopotential signals collected at multiple locations on a pregnant mother's abdomen, detecting R-wave peaks in the biopotential signals, extracting maternal ECG signals from the biopotential signals, determining R-wave amplitudes in the maternal ECG signals, creating an R-wave amplitude signal for each of the maternal ECG signals, calculating an average of all the R-wave amplitude signals, and normalizing the average to generate an EUM signal. In some embodiments, the method also includes identifying at least one uterine contraction based on a corresponding at least one peak in the EUM signal.

[0009] In one embodiment, a computer-implemented method includes receiving, by at least one computer processor, a plurality of raw biopotential inputs, each raw biopotential input received from a corresponding one of a plurality of electrodes, each of the plurality of electrodes positioned to measure a respective one of the raw biopotential inputs of the pregnant human subject; generating, by the at least one computer processor, a plurality of signal channels from the plurality of raw biopotential inputs, the plurality of signal channels including at least three signal channels; pre-processing, by the at least one computer processor, signal channel data of each of the signal channels to create a plurality of pre-processed signal channels, each of the pre-processed signal channels including a respective pre-processed signal channel data; and extracting, by the at least one computer processor, a plurality of respective R-wave peaks from the pre-processed signal channel data of each of the pre-processed signal channels to create a plurality of R-wave peak data sets, removing, by at least one computer processor, at least one of (a) at least one signal artifact or (b) at least one outlier data point from the plurality of R-wave peak data sets, where the at least one signal artifact is either an electromyogram artifact or a baseline artifact; replacing, by the at least one computer processor, the at least one signal artifact, the at least one outlier data point, or both, with at least one statistical value determined based on a corresponding one of the R-wave peak data sets from which the at least one signal artifact, the at least one outlier data point, or both have been removed; generating, by the at least one computer processor, a respective R-wave signal data set for each R-wave signal channel at a predetermined sampling rate based on the respective R-wave peak data sets, thereby creating a plurality of R-wave signal channels;selecting at least one first specific R-wave signal channel and at least one second specific R-wave signal channel from the plurality of R-wave channels based on at least one correlation between (a) a respective R-wave signal data set of at least one first specific R-wave signal channel and (b) a respective R-wave signal data set of at least one second specific R-wave signal channel; and generating, by at least one computer processor, electrical uterine monitoring data representative of the electrical uterine monitoring signal based at least on the respective R-wave signal data set of the first selected R-wave signal channel and the respective R-wave signal data set of the second selected R-wave signal channel.

[0010] In embodiments, the computer-implemented method also includes sharpening, by at least one computer processor, the electrical uterine monitoring data to create a sharpened electrical uterine monitoring signal. In embodiments, if the electrical uterine monitoring data is calculated based on a selected one of the electrical uterine monitoring signal channels that is the corrupted electrical uterine monitoring channel, the sharpening is omitted. In embodiments, the computer-implemented method also includes post-processing the sharpened electrical uterine monitoring signal data to create a post-processed electrical uterine monitoring signal. In embodiments, the sharpening includes identifying a set of peaks in the electrical uterine monitoring signal data, determining a prominence of each of the peaks, removing from the set of peaks peaks that have a prominence less than at least one threshold prominence value, calculating a mask based on the remaining peaks of the set of peaks, smoothing the mask based on a moving average window to create a smoothed mask, and adding the smoothed mask to the electrical uterine monitoring signal data to create the sharpened electrical uterine monitoring signal data. In an embodiment, the at least one threshold prominence value comprises at least one threshold prominence value selected from the group consisting of: an absolute prominence value and a relative prominence value calculated based on a maximum prominence of the peaks in the set of peaks, In an embodiment, the mask comprises zero values ​​outside the region of the remaining peaks and non-zero values ​​inside the region of the remaining peaks, the non-zero values ​​being calculated based on a Gaussian function.

[0011] In an embodiment, the filtering step of at least one of the pre-processing steps comprises applying at least one filter selected from the group consisting of a DC removal filter, a power line filter, and a high pass filter.

[0012] In an embodiment, the extracting step includes receiving a set of maternal ECG peaks of the pregnant human subject and identifying an R-wave peak of each of the pre-processed signal channels within a predetermined time window before and after each of the maternal ECG peaks in the set of maternal ECG peaks as the maximum absolute value of each of the pre-processed signal channels within the predetermined time window.

[0013] In an embodiment, the step of removing at least one of the signal artifacts or outlier data points includes removing at least one electromyographic artifact by a process that includes identifying at least one corrupted peak in one of the plurality of R-wave peak data sets based on the at least one corrupted peak having a peak-to-peak root mean square value greater than a threshold, and replacing the corrupted peak with a median, where the median is one of a local median or a global median.

[0014] In an embodiment, the step of removing at least one of the signal artifacts or outlier data points includes removing at least one baseline artifact by a process including: identifying a change point of an R-wave peak in one of the plurality of R-wave peak data sets, subdividing the one of the plurality of R-wave peak data sets into a first portion located before the change point and a second portion located after the change point, determining a first root-mean-square value for the first portion, determining a second root-mean-square value for the second portion, determining an equalization factor based on the first root-mean-square value and the second root-mean-square value, and modifying the first portion by multiplying the R-wave peak in the first portion by the equalization factor.

[0015] In an embodiment, the step of removing at least one of signal artifacts or outlier points comprises removing at least one outlier according to a Grubbs test for outliers.

[0016] In an embodiment, generating each R-wave data set based on each R-wave peak data set includes interpolating between R-wave peaks of each R-wave peak data set, and interpolating between the R-wave peaks includes interpolating using an interpolation algorithm selected from the group consisting of a cubic spline interpolation algorithm and a shape-preserving piecewise cubic interpolation algorithm.

[0017] In an embodiment, the step of selecting at least one first R-wave signal channel and at least one second R-wave signal channel includes: selecting candidate R-wave signal channels from the R-wave signal channels based on a percentage of previous intervals in which each of the R-wave signal channels experienced a contact problem; grouping the selected candidate R-wave signal channels into a plurality of couples, each couple including two of the selected candidate R-wave channels that are independent of each other; calculating a correlation value for each of the couples; and selecting at least one candidate R-wave signal channel of the couple as the selected at least one first R-wave signal channel and the selected at least one second R-wave signal channel based on at least one of the couples having a correlation value that exceeds a threshold correlation value.

[0018] In an embodiment, the step of calculating the electrical uterine monitoring signal includes calculating a signal that is a predetermined percentile of at least one first selected one of the R-wave signal channels and at least one second selected one of the R-wave signal channels, in an embodiment, the predetermined percentile is the 80th percentile.

[0019] In an embodiment, the statistic is one of a local median, a global median, or a mean.

[0020] In some embodiments, a computer-implemented method includes providing, by at least one computer processor, a plurality of signal channels, the plurality of signal channels including a plurality of electrical uterine monitoring signal channels and a plurality of acoustic uterine monitoring signal channels; determining, by the at least one computer processor, a plurality of channel weights, each of the channel weights corresponding to a particular one of the signal channels; and generating a combined uterine monitoring signal channel by calculating, by the at least one computer processor, a weighted average of the signal channels based on the channel weight for each of the signal channels.

[0021] In some embodiments, the plurality of channel weights are determined based on a machine learning algorithm, hi some embodiments, the machine learning algorithm comprises a gradient descent optimization process.

[0022] In some embodiments, the plurality of channel weights are determined by a process that includes: defining, by at least one computer processor, a plurality of channel sets, each of the plurality of channel sets including at least some of the plurality of signal channels; defining, by the at least one computer processor, a plurality of initial weight sets, each of the plurality of initial weight sets corresponding to a particular one of the plurality of channel sets; optimizing, by the at least one computer processor, the plurality of initial weight sets to generate a plurality of optimized weight sets, each of the plurality of optimized weight sets corresponding to a particular one of the plurality of channel sets; and selecting, by the at least one computer processor, a best one of the plurality of optimized weight sets as the plurality of channel weights.

[0023] In some embodiments, optimizing the multiple initial weight sets comprises a gradient descent process.

[0024] In some embodiments, the step of selecting the best one of the multiple optimized weight sets is performed by a process that includes: generating, by at least one computer processor, a plurality of intermediate uterine activity traces, each of the plurality of intermediate uterine activity traces corresponding to a particular one of a plurality of optimized weight sets; calculating, by the at least one computer processor, for each of the plurality of optimized weight sets: (a) a signal-to-noise ratio of a particular one of the intermediate uterine activity traces corresponding to each of the plurality of optimized weight sets, (b) a cost function, (c) a contraction reliability index, and (d) a dissimilarity index; calculating, by the at least one computer processor, an optimized weight set mean for each of the plurality of optimized weight sets, the optimized weight set mean being an average of (a) the signal-to-noise ratio of the particular one of the optimized weight sets, (b) the cost function of the particular one of the optimized weight sets, (c) the contraction reliability index of the particular one of the optimized weight sets, and (d) the dissimilarity index of the particular one of the optimized weight sets; and selecting, by the at least one computer processor, the optimized weight set having the best optimized weight set mean from the plurality of optimized weight sets as the best one of the plurality of optimized weight sets.

[0025] In some embodiments, the computer-implemented method also includes generating, by at least one computer processor, a first intermediate uterine activity trace and a second intermediate uterine activity trace corresponding to a particular one of a plurality of channel sets, wherein the first intermediate uterine activity trace corresponds to a first one of a plurality of optimized weight sets for the particular one of the plurality of channel sets and the second intermediate uterine activity trace corresponds to a second one of a plurality of optimized weight sets for the particular one of the plurality of channel sets; calculating, by the at least one computer processor, for the first one of the plurality of optimized weight sets: (a) a signal-to-noise ratio of the first intermediate uterine activity trace, (b) a cost function, (c) a contraction reliability index, and (d) a difference index; and calculating, by the at least one computer processor, for the second one of the plurality of optimized weight sets: (a) a signal-to-noise ratio of the second intermediate uterine activity trace, (b) a cost function, (c) a contraction reliability index, and (d) a difference index. calculating, for a first of the plurality of optimized weight sets, a first average that is the average of (a) the signal-to-noise ratio of the first intermediate uterine activity traces, (b) the cost function of the first of the plurality of optimized weight sets, (c) the contraction reliability index of the first of the plurality of optimized weight sets, and (d) the difference index of the first of the plurality of optimized weight sets; and calculating, by at least one computer processor, a second average that is the average of (a) the signal-to-noise ratio of the second intermediate uterine activity traces, (b) the cost function of the second of the plurality of optimized weight sets, (c) the contraction reliability index of the second of the plurality of optimized weight sets, and (d) the difference index of the second of the plurality of optimized weight sets; and selecting, by the at least one computer processor, the first of the plurality of weight sets as the best weight set for the particular one of the plurality of channel sets based on a determination that the first average is better than the second average;and selecting, by the at least one computer processor, a second one of the plurality of weight sets as the best weight set for the particular one of the plurality of channel sets based on determining that the second average is better than the first average. In some embodiments, the computer-implemented method also includes, prior to the step of selecting, by the at least one computer processor, the best one of the plurality of optimized weight sets as the plurality of channel weights, by at least the computer processor.

[0026] In some embodiments, defining the plurality of channel sets includes defining a contraction-based channel set, where the contraction-based channel set is determined by a process including: identifying, by at least one computer processor, a set of contractions in each of a plurality of signal channels; extracting, by at least one computer processor, contraction features for each of the plurality of signal channels based on the set of contractions identified for each of the plurality of signal channels; clustering, by at least one computer processor, the plurality of signal channels into a plurality of clusters; and selecting, by at least one computer processor, a best one of the plurality of clusters as the contraction-based channel set. In some embodiments, defining the plurality of channel sets further includes refining, by at least one computer processor, the best one of the plurality of clusters. In some embodiments, defining the plurality of channel sets also includes adding, by at least one computer processor, a portion of one of the signal channels that is not included in the best one of the plurality of clusters to the best one of the plurality of clusters. In some embodiments, identifying a set of contractions in each of the plurality of signal channels is performed by a process including, for each one of the plurality of signal channels, generating, by at least one computer processor, an enhanced version of one of the plurality of signal channels; detecting, by at least one computer processor, a candidate set of contractions in the enhanced one of the plurality of signal channels, the candidate set of contractions including a plurality of contraction candidates; calculating, by at least one computer processor, a plurality of confidence indicators for each candidate contraction; and removing, by at least one computer processor, at least one of the candidate contractions from the set of candidate contractions based on the confidence indicator corresponding to the removed one of the at least one candidate contractions, the removal generating a set of contractions.

[0027] In some embodiments, the step of defining, by the at least one computer processor, the plurality of initial weight sets includes generating, by the at least one computer processor, a channel voting weight set and a natural equivalent weight set for each of the channel sets.

[0028] In some embodiments, the step of providing a plurality of signal channels includes generating, by at least one computer processor, at least one of the plurality of electrical uterine monitoring signal channels, wherein at least one of the plurality of electrical uterine monitoring signal channels is generated by a process including: receiving, by the at least one computer processor, a plurality of raw biopotential inputs, each of the raw biopotential inputs received from a corresponding one of a plurality of electrodes, each of the plurality of electrodes positioned to measure a respective one of the raw biopotential inputs of the pregnant human subject; generating, by the at least one computer processor, a plurality of signal channels from the plurality of raw biopotential inputs, the plurality of signal channels including at least three signal channels; pre-processing, by the at least one computer processor, respective signal channel data for each of the signal channels to create a plurality of pre-processed signal channels, each of the pre-processed signal channels including a respective pre-processed signal channel data; extracting, by the at least one computer processor, a respective plurality of R-wave peaks from the pre-processed signal channel data for each of the pre-processed signal channels to create a plurality of R-wave peak data sets, each of the R-wave peak data sets including a respective plurality of R-wave peaks; removing, by the at least one computer processor, at least one of: (a) at least one signal artifact; or (b) at least one outlier data point from the plurality of R-wave peak data sets, wherein the at least one signal artifact is one of an electromyogram artifact or a baseline artifact; replacing, by the at least one computer processor, the at least one signal artifact, the at least one outlier data point, or both, with at least one statistical value determined based on a corresponding one of the R-wave peak data sets from which the at least one signal artifact, the at least one outlier data point, or both have been removed, to create a plurality of interpolated R-wave peak data sets; generating, by the at least one computer processor, a respective R-wave signal data set for each R-wave signal channel at a predetermined sampling rate based on each interpolated R-wave peak data set to create a plurality of R-wave signal channels; selecting, by the at least one computer processor, at least one first selected R-wave signal channel and at least one second selected R-wave signal channel from the plurality of R-wave signal channels based on at least one correlation between (a) a respective R-wave signal data set of at least one first specific R-wave signal channel and (b) a respective R-wave signal data set of at least one second specific R-wave signal channel; generating, by the at least one computer processor, electrical uterine monitoring data representative of an electrical uterine monitoring signal based at least on the respective R-wave signal data sets of the first selected R-wave signal channel and the respective R-wave signal data sets of the second selected R-wave signal channel, thereby creating the at least one electrical uterine monitoring signal channel.

[0029] In some embodiments, the step of providing a plurality of signal channels includes generating, by at least one computer processor, at least one of the plurality of acoustic uterine monitoring signal channels, wherein at least one of the plurality of acoustic uterine monitoring signal channels is generated by a process including: receiving, by the at least one computer processor, a plurality of raw acoustic inputs, each of the raw acoustic inputs received from a corresponding one of the plurality of acoustic sensors, each of the plurality of acoustic sensors positioned to measure a respective one of the raw acoustic inputs of the pregnant human subject; generating, by the at least one computer processor, a plurality of signal channels from the plurality of raw acoustic inputs, the plurality of signal channels including at least three signal channels; pre-processing, by the at least one computer processor, respective signal channel data for each of the signal channels to create a plurality of pre-processed signal channels, each of the pre-processed signal channels including a respective pre-processed signal channel data; extracting, by the at least one computer processor, a respective plurality of S1-S2 peaks from the pre-processed signal channel data for each of the pre-processed signal channels to create a plurality of S1-S2 peak data sets, each of the S1-S2 peak data sets including a respective plurality of S1-S2 peaks; removing, by the at least one computer processor, at least one of: (a) at least one signal artifact; or (b) at least one outlier data point from the plurality of S1-S2 peak data sets, wherein the at least one signal artifact is one of a motion-related artifact or a baseline artifact; replacing, with the at least one computer processor, the at least one signal artifact, the at least one outlier data point, or both, with at least one statistical value determined based on a corresponding one of the S1-S2 peak data sets from which the at least one signal artifact, the at least one outlier data point, or both have been removed, to create a plurality of interpolated S1-S2 peak data sets; generating, by the at least one computer processor, a respective S1-S2 signal data set for each S1-S2 signal channel at a predetermined sampling rate based on the respective interpolated S1-S2 peak data set to create a plurality of S1-S2 signal channels; selecting, by the at least one computer processor, at least one first selected S1-S2 signal channel and at least one second selected S1-S2 signal channel from the plurality of S1-S2 signal channels based on at least one correlation between (a) a respective S1-S2 signal data set of at least one first specific S1-S2 signal channel and (b) a respective S1-S2 signal data set of at least one second specific S1-S2 signal channel; and generating, by the at least one computer processor, acoustic uterine monitoring data representative of an acoustic uterine monitoring signal based at least on the respective S1-S2 signal data sets of the first selected S1-S2 signal channel and the respective S1-S2 signal data sets of the second selected S1-S2 signal channel, thereby creating the at least one acoustic uterine monitoring signal channel. [Brief explanation of the drawings]

[0030] [Figure 1A] FIG. 1A shows a representative uterus in a non-contractile state. [Figure 1B] FIG. 1B shows a representative uterus in a contracted state. [Figure 2]FIG. 2 is a flow chart of an exemplary method. [Figure 3] FIG. 3 illustrates an exemplary garment including multiple biopotential sensors that may be used to sense data to be analyzed according to the exemplary method of FIG. [Figure 4A] FIG. 4A is a front view illustrating the positioning of ECG sensor pairs on the abdomen of a pregnant woman according to some embodiments of the present invention. [Figure 4B] FIG. 4B is a side view of the location of an ECG sensor pair on the abdomen of a pregnant woman, according to some embodiments of the present invention. [Figure 5] FIG. 5 shows an exemplary biopotential signal before and after preprocessing. [Figure 6A] FIG. 6A is an exemplary biopotential signal after preprocessing, showing the detected R-wave peak. [Figure 6B] FIG. 6B shows the example biopotential signal of FIG. 6A after peak redetection. [Figure 6C] FIG. 6C shows an expanded view of a portion of the signal in FIG. 6B. [Figure 6D] FIG. 6D illustrates the example biopotential signal of FIG. 6B after inspection of the detected peaks. [Figure 7A] FIG. 7A shows a portion of an exemplary biopotential signal with an identified R-wave peak. [Figure 7B] FIG. 7B shows a portion of an exemplary biopotential signal with the P wave, QRS complex, and T wave identified therein. [Figure 7C] FIG. 7C shows an exemplary biopotential signal containing mixed maternal and fetal data. [Figure 7D] FIG. 7D shows a portion of the signal of FIG. 7C together with the initial template. [Figure 7E] FIG. 7E shows a portion of the signal of FIG. 7C with an adapted template. [Figure 7F] FIG. 7F shows a portion of the signal of FIG. 7C with the current template and the zeroth iteration of adaptation. [Figure 7G]FIG. 7G shows a portion of the signal of FIG. 7C along with the current template, the zeroth iteration of adaptation, and the first iteration of adaptation. [Figure 7H] FIG. 7H shows a portion of the signal of FIG. 7C along with the current template and a maternal ECG signal reconstructed based on the current template. [Figure 7I] FIG. 7I shows the progress of the adaptation, represented by the logarithm of the error signal plotted against the number of iterations. [Figure 7J] FIG. 7J shows the extracted maternal ECG signal. [Figure 8] FIG. 8 shows an exemplary filtered maternal ECG signal. [Figure 9] FIG. 9 shows an exemplary maternal ECG signal with the R-wave peaks annotated. [Figure 10A] FIG. 10A shows an exemplary R-wave amplitude signal. [Figure 10B] FIG. 10B shows an exemplary modulated R-wave amplitude signal. [Figure 11A] FIG. 11A illustrates an exemplary modulated R-wave amplitude signal along with the results of applying a moving average filter. [Figure 11B] FIG. 11B shows exemplary filtered R-wave amplitude signals for multiple channels over the same time window. [Figure 12A] FIG. 12A is an illustration of a first exemplary normalized electrical uterine signal generated based on the exemplary filtered R-wave amplitude signal shown in FIG. 11B. [Figure 12B] FIG. 12B shows the first tocograph signal recorded over the same time period as an exemplary normalized electrical uterine signal, with self-reported contractions annotated. [Figure 13] FIG. 13 is a flowchart of a second exemplary method. [Figure 14A] FIG. 14A shows a second tocograph signal with self-reported contractions annotated. [Figure 14B]FIG. 14B shows a second exemplary uterine electrical signal obtained from biopotential data recorded during the same time period as shown in FIG. 14A. [Figure 15A] FIG. 15A shows a third tocograph signal with self-reported contractions annotated. [Figure 15B] FIG. 15B shows a third exemplary uterine electrical signal derived from biopotential data recorded during the same time period as shown in FIG. 15A. [Figure 16A] FIG. 16A shows a fourth tocograph signal with self-reported contractions annotated. [Figure 16B] FIG. 16B shows a fourth exemplary uterine electrical signal derived from biopotential data recorded during the same time period as shown in FIG. 16A. [Figure 17A] FIG. 17A shows a fifth tocograph signal with self-reported contractions annotated. [Figure 17B] FIG. 17B shows a fifth exemplary uterine electrical signal derived from biopotential data recorded during the same time period as shown in FIG. 17A. [Figure 18A] FIG. 18A shows an exemplary raw biopotential data set. [Figure 18B] FIG. 18B shows an exemplary filtered data set based on the exemplary raw data set of FIG. 18A. [Figure 18C] FIG. 18C shows an exemplary raw biopotential data set. [Figure 18D] FIG. 18D shows an exemplary filtered dataset based on the exemplary raw dataset of FIG. 18C. [Figure 18E] FIG. 18E shows an exemplary raw biopotential data set. [Figure 18F] FIG. 18F shows an exemplary filtered dataset based on the exemplary raw dataset of FIG. 18E. [Figure 18G] FIG. 18G shows an exemplary raw biopotential data set. [Figure 18H]FIG. 18H shows an exemplary filtered dataset based on the exemplary raw dataset of FIG. 18G. [Figure 19A] FIG. 19A shows an exemplary filtered data set with input peak locations. [Figure 19B] FIG. 19B shows the exemplary filtered data set of FIG. 19A with the peak locations extracted. [Figure 20A] FIG. 20A shows an exemplary filtered data set. [Figure 20B] FIG. 20B shows the exemplary filtered data set of FIG. 20A with the corresponding maternal motor envelope and peak-to-peak absolute sum representations. [Figure 20C] FIG. 20C shows an exemplary corrected data set produced by removing the electromyographic artifact from the filtered data set of FIG. 20A. [Figure 21A] FIG. 21A shows an exemplary corrected data set including baseline artifacts. [Figure 21B] FIG. 21B shows the exemplary corrected data set of FIG. 21A following removal of the baseline artifact. [Figure 22A] FIG. 22A shows an exemplary corrected data set including an outlier data point. [Figure 22B] FIG. 22B shows the exemplary corrected data set of FIG. 22A following removal of the outlier data points. [Figure 23A] FIG. 23A shows an exemplary R-wave peak signal. [Figure 23B] FIG. 23B is an exemplary R-wave signal generated based on the exemplary R-wave peak signal of FIG. 23A. [Figure 24A] FIG. 24A shows an exemplary set of candidate R-wave signal channels. [Figure 24B] FIG. 24B shows an exemplary set of selected signal channels based on the exemplary set of candidate R-wave signal channels of FIG. 24A. [Figure 25A]FIG. 25A illustrates an exemplary electrical uterine monitoring signal generated based on the set of selected signal channels shown in FIG. 24B. [Figure 25B] FIG. 25B shows an exemplary corrected electrical uterine monitoring signal produced by applying wandering baseline removal to the exemplary electrical uterine monitoring signal of FIG. 25A. [Figure 26] FIG. 26 illustrates an exemplary normalized electrical uterine monitoring signal generated based on the exemplary corrected electrical uterine monitoring signal of FIG. 25B. [Figure 27A] FIG. 27A shows an exemplary normalized electrical uterine monitoring signal. [Figure 27B] FIG. 27B illustrates an exemplary sharpening mask generated based on the exemplary normalized electrical uterine monitoring signal of FIG. 27A. [Figure 27C] FIG. 27C shows an example sharpened electrical uterine monitoring signal generated based on the example normalized electrical uterine monitoring signal of FIG. 27A and the example sharpening mask of FIG. 27B. [Figure 28] FIG. 28 illustrates an exemplary post-processed electrical uterine monitoring signal. [Figure 29] FIG. 29 shows a tocograph signal corresponding to the exemplary processed electrical uterine monitoring signal of FIG. [Figure 30] FIG. 30 is a flowchart of a third exemplary method. [Figure 31A] FIG. 31A shows an exemplary pre-processed data set. [Figure 31B] FIG. 31B shows an expanded view of a portion of the exemplary pre-processed data set of FIG. 31A. [Figure 32A] FIG. 32A illustrates extracted R-wave peaks in an exemplary pre-processed data set. [Figure 32B] FIG. 32B is a zoomed-in view of an extracted R-wave peak in an exemplary pre-processed data set. [Figure 32C]FIG. 32C illustrates an exemplary R-wave amplitude signal. [Figure 32D] FIG. 32D shows an exemplary R-wave amplitude signal over a larger time window. [Figure 33] FIG. 33 shows an exemplary R-wave amplitude signal after filtering. [Figure 34] FIG. 34 shows four data channels of exemplary R-wave data. [Figure 35A] FIG. 35A shows a sixth tochograph signal. [Figure 35B] FIG. 35B shows a first exemplary acoustic uterine signal obtained from acoustic data recorded during the same time period as shown in FIG. 35A. [Figure 36A] FIG. 36A shows the seventh tochograph signal. [Figure 36B] FIG. 36B shows a second exemplary acoustic uterine signal obtained from acoustic data recorded during the same time period as shown in FIG. 36A. [Figure 37A] FIG. 37A shows an eighth tochograph signal. [Figure 37B] FIG. 37B shows a third exemplary acoustic uterine signal obtained from acoustic data recorded during the same time period as shown in FIG. 37A. [Figure 38] FIG. 38 shows a set of ECG-based EUM processed signals and PCG-based processed signals. [Figure 39] FIG. 39 shows the fusion process for generating a uterine activity signal from an ECG-based EUM processed signal and a PCG-based processed signal. [Figure 40] FIG. 40 illustrates a process for generating a weight set for use in generating a uterine activity signal based on electrical and acoustic uterine monitoring data. [Figure 41] FIG. 41 illustrates a process for defining an initial channel set for use in the method of FIG. [Figure 42] FIG. 42 illustrates a process for identifying contractions in a data set. [Figure 43]FIG. 43 illustrates an exemplary uterine monitoring signal and an exemplary smoothed and enhanced signal produced by FIG. [Figure 44] FIG. 44 illustrates an exemplary electrical uterine monitoring data signal, an exemplary acoustic uterine monitoring data signal, and an exemplary output uterine monitoring signal generated by the method of FIG. DETAILED DESCRIPTION OF THE INVENTION

[0031] Among those benefits and improvements that have been disclosed, other objects and advantages of the present invention will become apparent from the following description taken in conjunction with the accompanying drawings. Although detailed embodiments of the present invention are disclosed herein, it should be understood that the disclosed embodiments are merely exemplary of the invention, which may be embodied in various forms. Additionally, the examples given in connection with various embodiments of the present invention are intended to be illustrative, not limiting.

[0032] Throughout this specification and claims, the following terms take the meanings expressly associated therewith herein, unless the context clearly dictates otherwise. As used herein, the phrases "in one embodiment," "in one embodiment," and "in some embodiments" do not necessarily refer to the same embodiment(s), although they may. Furthermore, as used herein, the phrases "in another embodiment" and "in some other embodiments" do not necessarily refer to different embodiments, although they may. Thus, as described below, various embodiments of the invention can be readily combined without departing from the scope or spirit of the invention.

[0033] As used herein, the term "based on" is not exclusive and allows for based on additional unrecited factors unless the context clearly dictates otherwise. Additionally, throughout this specification, the meanings of "a," "an," and "the" include plural references. Also, the meaning of "in" includes "in" and "on." Ranges discussed herein are inclusive (e.g., a range "between 0 and 2" includes not only the values ​​0 and 2, but all values ​​therebetween).

[0034] As used herein, the term "contact area" includes the contact area between the skin of the pregnant human subject and the skin contact, i.e., the surface area through which the flow of electrical current can pass between the skin of the pregnant human subject and the skin contact.

[0035] In some embodiments, the present invention provides methods for extracting tocograph-like signals from biopotential data, i.e., data describing electrical potentials recorded at points on a person's skin through the use of skin-contacting objects, commonly referred to as electrodes. In some embodiments, the present invention provides methods for detecting uterine contractions from biopotential data. In some embodiments, the biopotential data is obtained through the use of non-contact electrodes placed at or near desired points on the person's body.

[0036] In some embodiments, the present invention provides a system for detecting, recording, and analyzing cardiac electrical activity data from a pregnant human subject. In some embodiments, a plurality of electrodes configured to detect fetal electrocardiogram signals are used to record the cardiac activity data. In some embodiments, a plurality of electrodes configured to detect fetal electrocardiogram signals and a plurality of acoustic sensors are used to record the cardiac activity data.

[0037] In some embodiments, a plurality of electrodes configured to detect fetal electrocardiogram signals are attached to the abdomen of a pregnant human subject. In some embodiments, a plurality of electrodes configured to detect fetal electrocardiogram signals are attached directly to the abdomen. In some embodiments, a plurality of electrodes configured to detect fetal electrocardiogram signals are incorporated into an article, such as a belt or patch, worn by or placed on a pregnant human subject. FIG. 3 illustrates an exemplary garment 300, which includes eight electrodes 310 incorporated into the garment 300 to be positioned around the abdomen of a pregnant human subject when the garment 300 is worn by the subject. In some embodiments, the garment 300 includes four acoustic sensors 320 incorporated into the garment 300 to be positioned around the abdomen of a pregnant human subject when the garment 300 is worn by the subject. In some embodiments, each of the acoustic sensors 320 is one of the acoustic sensors described in U.S. Patent No. 9,713,430. FIG. 4A is a front view of the positioning of eight electrodes 310 on a pregnant human subject's abdomen, according to some embodiments of the present invention. FIG. 4B is a side view of eight electrodes 310 on the abdomen of a pregnant woman, according to some embodiments of the present invention.

[0038] FIG. 2 is a flowchart of an exemplary method 200 of the present invention. In some embodiments, an exemplary computing device of the present invention programmed / configured according to method 200 is operable to receive as input raw biopotential data measured by a plurality of electrodes placed on the skin of a pregnant human subject and analyze such input to generate a tocograph-like signal. In some embodiments, the number of electrodes is between 2 and 10. In some embodiments, the number of electrodes is between 2 and 20. In some embodiments, the number of electrodes is between 2 and 30, inclusive. In some embodiments, the number of electrodes is between 2 and 40, inclusive. In some embodiments, the number of electrodes is between 4 and 10. In some embodiments, the number of electrodes is between 4 and 20. In some embodiments, the number of electrodes is between 4 and 30. In some embodiments, the number of electrodes is between 4 and 40. In some embodiments, the number of electrodes is between 6 and 10. In some embodiments, the number of electrodes is between 6 and 20. In some embodiments, the number of electrodes is between 6 and 30. In some embodiments, the number of electrodes is between 6 and 40. In some embodiments, the quantity of electrodes is between 8 and 10. In some embodiments, the quantity of electrodes is between 8 and 20. In some embodiments, the quantity of electrodes is between 8 and 30. In some embodiments, the quantity of electrodes is between 8 and 40. In some embodiments, the quantity of electrodes is 8. In some embodiments, an exemplary computing device of the present invention programmed / configured according to method 200 is operable to receive as input a maternal ECG signal that has already been extracted from the raw biopotential data (e.g., by separation from a fetal ECG signal that forms part of the same raw biopotential data). In some embodiments, an exemplary computing device of the present invention is programmed / configured according to method 200 via instructions stored on a non-transitory computer-readable medium. In some embodiments, an exemplary computing device of the present invention includes at least one computer processor that, when executing the instructions, becomes a specific computer processor programmed / configured according to method 200.

[0039] In some embodiments, an exemplary computing device of the present invention is programmed / configured to continuously perform one or more steps of method 200 along a moving time window. In some embodiments, the moving time window has a predefined length. In some embodiments, the predefined length is 60 seconds. In some embodiments, an exemplary computing device of the present invention is programmed / configured to continuously perform one or more steps of method 200 along a moving time window having a length between 1 second and 1 hour. In some embodiments, the length of the moving time window is between 30 seconds and 30 minutes. In some embodiments, the length of the moving time window is between 30 seconds and 10 minutes. In some embodiments, the length of the moving time window is between 30 seconds and 5 minutes. In some embodiments, the length of the moving time window is approximately 60 seconds. In some embodiments, the length of the moving time window is 60 seconds.

[0040] In step 210, an exemplary computing device of the present invention is programmed / configured to receive raw biopotential data as input and preprocess it. In some embodiments, the raw biopotential data is recorded through the use of at least two electrodes placed in proximity to the skin of the pregnant subject. In some embodiments, at least one of the electrodes is a signal electrode. In some embodiments, at least one of the electrodes is a reference electrode. In some embodiments, the reference electrode is placed at a point distant from the subject's uterus. In some embodiments, biopotential signals are recorded at each of multiple points around the pregnant subject's abdomen. In some embodiments, biopotential signals are recorded at each of eight points around the pregnant subject's abdomen. In some embodiments, the biopotential data is recorded at 1,000 samples per second. In some embodiments, the biopotential data is upsampled to 1,000 samples per second. In some embodiments, the biopotential data is recorded at a sampling rate between 100 and 10,000 samples per second. In some embodiments, the biopotential data is upsampled to a sampling rate between 100 and 10,000 samples per second. In some embodiments, pre-processing includes baseline removal (e.g., using a median filter and / or a moving average filter). In some embodiments, pre-processing includes low-pass filtering. In some embodiments, pre-processing includes low-pass filtering at 85 Hz. In some embodiments, pre-processing includes power line interference cancellation. FIG. 5 shows a portion of the raw biopotential data signal both before and after pre-processing.

[0041] In step 220, an exemplary computing device of the present invention is programmed / configured to detect maternal R-wave peaks in the preprocessed biopotential data resulting from the execution of step 210. In some embodiments, R-wave peaks are detected over 10-second segments of each data signal. In some embodiments, R-wave peak detection begins with derivative, threshold, and distance analysis. In some embodiments, detecting R-wave peaks in each data signal involves calculating the first derivative of the data signal over the 10-second segment, identifying R-wave peaks over the 10-second segment by identifying zero crossings of the first derivative, and excluding identified peaks having either (a) an absolute value less than a predetermined R-wave peak threshold absolute value or (b) a distance between adjacent identified R-wave peaks less than a predetermined R-wave peak threshold distance. In some embodiments, R-wave peak detection is performed in a manner similar to the electrocardiogram peak detection described in U.S. Pat. No. 9,392,952 (Patent Document 2), the entire contents of which are incorporated herein by reference. FIG. 6A illustrates a preprocessed biopotential data signal, with R-wave peaks detected as described above indicated by asterisks.

[0042] In some embodiments, the detection of R-wave peaks in step 220 continues with a peak re-detection process. In some embodiments, the peak re-detection process includes an automatic gain control ("AGC") analysis to detect windows with significantly different numbers of peaks. In some embodiments, the peak re-detection process includes a cross-correlation analysis. In some embodiments, the peak re-detection process includes an AGC analysis and a cross-correlation analysis. In some embodiments, the AGC analysis is appropriate for overcoming false negatives. In some embodiments, the cross-correlation analysis is appropriate for removing false positives. FIG. 6B shows the data signal after peak re-detection, with the R-wave peaks re-detected as described above indicated by asterisks. FIG. 6C shows an expanded view of a portion of the data signal of FIG. 6B.

[0043] In some embodiments, the detection of R-wave peaks in step 220 continues with the construction of a global peak array. In some embodiments, the global peak array is created from multiple channels of data (e.g., each channel corresponding to one or more of the electrodes 310). In some embodiments, the signal in each channel is given a quality score based on the relative energy of the peak. In some embodiments, the relative energy of the peak refers to the energy of the peak relative to the total energy of the signal being processed. In some embodiments, the energy of the peak is calculated by calculating the root mean square (“RMS”) of the QRS complex containing the R-wave peak, and the energy of the signal is calculated by calculating the RMS of the signal. In some embodiments, the relative energy of the peak is calculated by calculating the signal-to-noise ratio of the signal. In some embodiments, the channel with the highest quality score is considered the “best lead.” In some embodiments, the global peak array is constructed based on the best lead, and signals from other channels are also considered based on a voting mechanism. In some embodiments, after the global peak array is constructed based on the best lead, each of the remaining channels “votes” for each peak. A channel votes positively (e.g., gives a voting value of "1") for a given peak if it is included in the global peak array constructed based on the best read if it includes such a peak (e.g., as detected in the peak detection described above), and votes negatively (e.g., gives a voting value of "0") if it does not include such a peak. Peaks that receive more votes are considered to be higher quality peaks. In some embodiments, if a peak has a number of votes equal to or greater than a threshold, it is retained in the global peak array. In some embodiments, the vote number threshold is half the total number of channels. In some embodiments, if a peak has a number of votes less than the threshold, additional tests are performed on the peak. In some embodiments, the additional tests include calculating the correlation of the peak of the best read channel with a template calculated as the average of all peaks.In some embodiments, if the correlation is greater than a first threshold correlation value, the peak is retained in the global peak array. In some embodiments, the first threshold correlation value is 0.9. In some embodiments, if the correlation is less than the first threshold correlation value, further correlations are calculated for all reads with positive votes for the peak (i.e., not just the best read peak). In some embodiments, if the further correlation is greater than a second threshold correlation value, the peak is retained in the global peak array; if the further correlation is less than the second threshold correlation value, the peak is removed from the global peak array. In some embodiments, the second threshold correlation value is 0.85.

[0044] In some embodiments, once created, the global peak array is validated using physiological measurements. In some embodiments, the validation is performed by an exemplary inventive computing device such as that described in U.S. Pat. No. 9,392,952, the entire contents of which are incorporated herein. In some embodiments, the physiological parameters include RR intervals, mean, and standard deviation; and heart rate and heart rate variability. In some embodiments, the validation includes cross-correlation to overcome false negatives. FIG. 6D illustrates a data signal following the creation and validation of the global peak array as described above. In FIG. 6D, the peaks indicated by circled asterisks represent R-wave peaks previously detected (e.g., as shown in FIG. 6A), while the circles without asterisks represent R-wave peaks detected by cross-correlation to overcome false negatives as described above.

[0045] In some embodiments, if the initial step of R-wave detection fails (i.e., if no R-wave peaks are detected over a given sample), an independent component analysis ("ICA") algorithm is applied to the data samples, and the previous portions of step 220 are repeated. In some embodiments, an exemplary ICA algorithm is, for example, but not limited to, the FAST ICA algorithm. In some embodiments, the FAST ICA algorithm is utilized in connection with, for example, Hyvarinen et al., "Independent Component Analysis: Algorithms and Applications," Neural Networks 13(4-5):411-430 (2000).

[0046] Continuing with reference to FIG. 2 , at step 230, an exemplary computing device of the present invention is programmed / configured to extract a maternal ECG signal from a signal containing both maternal and fetal data. In some embodiments, if an exemplary computing device of the present invention programmed / configured to perform method 200 receives as input a maternal ECG signal after extraction from mixed maternal-fetal data, the exemplary computing device of the present invention is programmed / configured to skip step 230. FIG. 7A illustrates a portion of a signal in which R-wave peaks have been identified and which contains both maternal and fetal signals. While not intending to be limited to a particular theory, a primary challenge involved in the process of extracting a maternal ECG signal is that each maternal heartbeat is different from all other maternal heartbeats. In some embodiments, this challenge is addressed by using an adaptive reconstruction scheme to identify each maternal heartbeat. In some embodiments, the extraction process begins by segmenting the ECG signal into three source signals. In some embodiments, this segmentation includes using a curve length transform to locate the P wave, QRS complex, and T wave. In some embodiments, the curve length transform is as described in Zong et al., "A QT Interval Detection Algorithm Based On ECG Curve Length Transform," Computers In Cardiology 33:377-380 (October 2006). Figure 7B is an exemplary ECG signal including these portions.

[0047] Following the curve length transformation, step 230 continues by extracting the maternal signal using an adaptive template. In some embodiments, template adaptation is used to isolate the current beat. In some embodiments, extraction of the maternal signal using an adaptive template is performed as described in U.S. Pat. No. 9,392,952, the entire contents of which are incorporated herein by reference. In some embodiments, this process involves starting with the current template and adapting it using an iterative process to arrive at the current beat. In some embodiments, multipliers are defined for each portion of the signal (i.e., the P wave, the QRS complex, and the T wave) (referred to as P_mult, QRS_mult, and T_mult, respectively). In some embodiments, a shifting parameter is also defined. In some embodiments, the extraction uses a Levenberg-Marquardt nonlinear least mean squares algorithm, as shown below.

number

[0048] In some embodiments, the cost function is as follows:

number

[0049] In the above equation, φ represents the current beat ECG and φ represents the reconstructed ECG. In some embodiments, the method provides a local, stable, and repeatable solution. In some embodiments, the iterations proceed until the relative remaining energy reaches a threshold. In some embodiments, the threshold is between 0 db and -40 db. In some embodiments, the threshold is between -10 db and -40 db. In some embodiments, the threshold is between -20 db and -40 db. In some embodiments, the threshold is between -30 db and -40 db. In some embodiments, the threshold is between -10 db and -30 db. In some embodiments, the threshold is between -10 db and -20 db. In some embodiments, the threshold is between -20 db and -40 db. In some embodiments, the threshold is between -20 db and -30 db. In some embodiments, the threshold is between -30 db and -40 db. In some embodiments, the threshold is between -25 db and -35 db. In some embodiments, the threshold is about -20 db. In some embodiments, the threshold is about -20 db.

[0050] FIG. 7C shows an example signal including mixed maternal and fetal data. FIG. 7D shows a portion of the signal of FIG. 7C with an initial template for comparison. FIG. 7E shows a portion of the signal of FIG. 7C with an adapted template for comparison. FIG. 7F shows a portion of the signal of FIG. 7C, a current template, and the 0th iteration of adaptation. FIG. 7G shows a portion of the signal of FIG. 7C, the current template, the 0th iteration of adaptation, and the 1st iteration of adaptation. FIG. 7H shows a portion of the signal of FIG. 7C, the current template, and a reconstructed ECG signal (e.g., maternal ECG signal) based on the current template. FIG. 7I shows the progress of adaptation in terms of the logarithm of the error signal versus the number of iterations. FIG. 7J shows an extracted maternal ECG signal.

[0051] Continuing with reference to FIG. 2, at step 240, an exemplary computing device of the present invention is programmed / configured to perform signal cleanup on the maternal signal extracted at step 230. In some embodiments, the cleanup at step 240 includes filtering. In some embodiments, the filtering includes baseline removal using a moving average filter. In some embodiments, the filtering includes low-pass filtering. In some embodiments, the low-pass filtering is performed between 25 Hz and 125 Hz. In some embodiments, the low-pass filtering is performed between 50 Hz and 100 Hz. In some embodiments, the low-pass filtering is performed at 75 Hz. FIG. 8 shows a portion of an exemplary filtered maternal ECG following performance of step 240.

[0052] Continuing with reference to FIG. 2, in step 250, an exemplary computing device of the present invention is programmed / configured to calculate R-wave amplitudes of the filtered maternal ECG signal resulting from performance of step 240. In some embodiments, the R-wave amplitudes are calculated based on the maternal ECG peaks detected in step 220 and the maternal ECG signal extracted in step 230. In some embodiments, step 250 includes calculating the amplitudes of various R-waves. In some embodiments, the amplitudes are calculated as the value (e.g., signal amplitude) of the maternal ECG signal at the location of each detected peak. FIG. 9 shows an exemplary extracted maternal ECG signal with R-wave peaks annotated with circles.

[0053] Continuing with reference to FIG. 2 , in step 260, an exemplary computing device of the present invention is programmed / configured to create an R-wave amplitude signal over time based on the R-wave amplitudes calculated in step 250. In some embodiments, the calculated R-wave peaks are not sampled uniformly over time. Thus, in some embodiments, step 260 is performed to resample the R-wave amplitudes in a manner such that they are sampled uniformly over time (e.g., such that the time difference between every two adjacent samples is constant). In some embodiments, step 260 is performed by connecting the R-wave amplitude values ​​calculated in step 250 and resampling the connected R-wave amplitude values. In some embodiments, the resampling involves interpolation using query points defined in time. In some embodiments, the interpolation involves linear interpolation. In some embodiments, the interpolation involves spline interpolation. In some embodiments, the interpolation involves cubic interpolation. In some embodiments, the query points define the points in time at which the interpolation should occur. FIG. 10A illustrates an exemplary R-wave amplitude signal as created in step 260 based on the R-wave amplitudes from step 250. In Figure 10A, the maternal electrocardiogram is similar to that shown in Figure 8, with the detected R-wave peaks indicated by circles and the R-wave amplitude signal being the curve connecting the circles. Figure 10B shows the modulation of the R-wave amplitude signal over a larger time window.

[0054] Continuing with reference to FIG. 2 , at step 270, an exemplary computing device of the present invention is programmed / configured to clean up the R-wave amplitude signal by applying a moving average filter. In some embodiments, the moving average filter is applied to clean up high-frequency variations in the R-wave amplitude signal. In some embodiments, the moving average filter is applied over a predetermined time window. In some embodiments, the time window has a length between 1 second and 10 minutes. In some embodiments, the time window has a length between 1 second and 1 minute. In some embodiments, the time window has a length between 1 second and 30 seconds. In some embodiments, the time window has a length of 20 seconds. FIG. 11A shows the R-wave amplitude signal of FIG. 10B , with the signal resulting from application of the moving average filter indicated by the bold line along the center of the R-wave amplitude signal. As mentioned above, in some embodiments, multiple channels of data are considered as input to method 200. FIG. 11B shows a plot of filtered R-wave amplitude signals for multiple channels over the same time window.

[0055] Continuing with reference to FIG. 2 , at step 280, the exemplary computing device of the present invention is programmed / configured to calculate an average signal of all filtered R-wave signals (e.g., as shown in FIG. 11B ) per unit time. In some embodiments, one average signal is calculated for each time point for which there are samples. In some embodiments, the average signal is the 80th percentile of all signals at each time point. In some embodiments, the average signal is the 85th percentile of all signals at each time point. In some embodiments, the average signal is the 90th percentile of all signals at each time point. In some embodiments, the average signal is the 95th percentile of all signals at each time point. In some embodiments, the average signal is the 99th percentile of all signals at each time point. In some embodiments, the result of this averaging is a single signal with uniform sampling over time. At step 290, the exemplary computing device of the present invention is programmed / configured to normalize the signals calculated at step 280. In some embodiments, the signals are normalized by dividing by a constant factor. In some embodiments, the constant factor is between 2 and 1000 volts. In some embodiments, the constant factor is 50 volts. FIG. 12A shows an exemplary normalized electrical uterine signal following performance of steps 280 and 290. FIG. 12B shows a tocograph signal generated over the same time period, with contractions self-reported by the mother indicated by vertical lines. With reference to FIGS. 12A and 12B, it can be seen that the peaks in the exemplary normalized electrical uterine signal of FIG. 12A coincide with the self-reported contractions shown in FIG. 12B. Thus, in some embodiments, a normalized electrical uterine monitoring ("EUM") signal (e.g., the signal shown in FIG. 12A) generated through performance of exemplary method 200 is suitable for use in identifying contractions. In some embodiments, contractions are identified by identifying peaks in the EUM signal.

[0056] In some embodiments, the present invention is directed to a specifically programmed computer system including at least the following components: a non-transitory memory, the memory electronically storing computer-executable program code; and at least one computer processor, which, when executing the program code, becomes a specifically programmed computing processor configured to perform at least the following operations: receiving a plurality of biopotential signals collected at a plurality of locations on a pregnant mother's abdomen, detecting R-wave peaks in the biopotential signals, extracting maternal electrocardiogram ("ECG") signals from the biopotential signals, determining R-wave amplitudes of the maternal ECG signals, creating an R-wave amplitude signal for each of the maternal ECG signals, calculating an average of all the R-wave amplitude signals, and normalizing the average to generate an electrical uterine monitoring ("EUM") signal. In some embodiments, the operations also include identifying at least one uterine contraction based on a corresponding at least one peak in the EUM signal.

[0057] FIG. 13 is a flowchart of an exemplary inventive method 1300. In some embodiments, an exemplary inventive computing device programmed / configured according to method 1300 is operable to receive as input raw biopotential data measured by a plurality of electrodes placed on the skin of a pregnant human subject and analyze such input to generate a tocograph-like signal. In some embodiments, the number of electrodes is between 2 and 10. In some embodiments, the number of electrodes is between 2 and 20. In some embodiments, the number of electrodes is between 2 and 30, inclusive. In some embodiments, the number of electrodes is between 2 and 40, inclusive. In some embodiments, the number of electrodes is between 4 and 10. In some embodiments, the number of electrodes is between 4 and 20. In some embodiments, the number of electrodes is between 4 and 30. In some embodiments, the number of electrodes is between 4 and 40. In some embodiments, the number of electrodes is between 6 and 10. In some embodiments, the number of electrodes is between 6 and 20. In some embodiments, the number of electrodes is between 6 and 30. In some embodiments, the number of electrodes is between 6 and 40. In some embodiments, the quantity of electrodes is between 8 and 10. In some embodiments, the quantity of electrodes is between 8 and 20. In some embodiments, the quantity of electrodes is between 8 and 30. In some embodiments, the quantity of electrodes is between 8 and 40. In some embodiments, the quantity of electrodes is 8. In some embodiments, an exemplary computing device of the present invention programmed / configured according to method 1300 is operable to receive as input a maternal ECG signal that has already been extracted from the raw biopotential data (e.g., by separation from a fetal ECG signal that forms part of the same raw biopotential data). In some embodiments, an exemplary computing device of the present invention is programmed / configured according to method 1300 via instructions stored on a non-transitory computer-readable medium. In some embodiments, an exemplary computing device of the present invention includes at least one computer processor that, when executing the instructions, becomes a specific computer processor programmed / configured according to method 1300.

[0058] In some embodiments, an exemplary computing device of the present invention is programmed / configured to continuously perform one or more steps of method 1300 along a moving time window. In some embodiments, the moving time window has a predefined length. In some embodiments, the predefined length is 60 seconds. In some embodiments, an exemplary computing device of the present invention is programmed / configured to continuously perform one or more steps of method 1300 along a moving time window having a length between 1 second and 1 hour. In some embodiments, the length of the moving time window is between 30 seconds and 30 minutes. In some embodiments, the length of the moving time window is between 30 seconds and 10 minutes. In some embodiments, the length of the moving time window is between 30 seconds and 5 minutes. In some embodiments, the length of the moving time window is approximately 60 seconds. In some embodiments, the length of the moving time window is 60 seconds.

[0059] In step 1305, an exemplary computing device of the present invention is programmed / configured to receive raw biopotential data as input. Exemplary raw biopotential data is shown in FIGS. 18A, 18C, 18E, and 18G. In some embodiments, the raw biopotential data is recorded through the use of at least two electrodes placed in proximity to the skin of the pregnant subject. In some embodiments, at least one of the electrodes is a signal electrode. In some embodiments, at least one of the electrodes is a reference electrode. In some embodiments, the reference electrode is placed at a point distant from the subject's uterus. In some embodiments, biopotential signals are recorded at each of multiple points around the pregnant subject's abdomen. In some embodiments, biopotential signals are recorded at each of eight points around the pregnant subject's abdomen. In some embodiments, the biopotential data is recorded at 1,000 samples per second. In some embodiments, the biopotential data is upsampled to 1,000 samples per second. In some embodiments, the biopotential data is recorded at a sampling rate between 100 and 10,000 samples per second. In some embodiments, the biopotential data is upsampled to a sampling rate of between 100 and 10,000 samples per second. In some embodiments, the steps of method 1300 between receiving raw data and selecting a channel (i.e., steps 1310 through 1335) are performed for each of a plurality of signal channels, each signal channel being generated by the exemplary invention computing device as the difference between biopotential signals recorded by a particular pair of electrodes. In some embodiments in which method 1300 is performed through the use of data recorded with electrodes positioned as shown in Figures 4A and 4B, the channels are identified as follows: Channel 1: A1-A4 Channel 2: A2-A3 Channel 3: A2-A4 Channel 4: A4-A3 Channel 5: B1-B3 Channel 6: B1-B2 Channel 7: B3-B2 Channel 8: A1-A3

[0060] In step 1310, the exemplary computing device of the present invention is programmed / configured to preprocess the signal channels determined based on the raw biopotential data to generate multiple preprocessed signal channels. In some embodiments, the preprocessing includes one or more filters. In some embodiments, the preprocessing includes two or more filters. In some embodiments, the preprocessing includes a DC removal filter, a power line filter, and a high-pass filter. In some embodiments, the DC removal filter removes the average of the raw data for the current processing interval. In some embodiments, the power line filter includes a 10th-order bandstop infinite impulse response ("IIR") filter configured to minimize any noise at a preset frequency in the data. In some embodiments, the preset frequency is 50 Hz, and the power line filter includes cutoff frequencies of 49.5 Hz and 50.5 Hz. In some embodiments, the preset frequency is 60 Hz, and the power line filter includes cutoff frequencies of 59.5 Hz and 60.5 Hz. In some embodiments, the high-pass filtering is performed by subtracting a wandering baseline from the signal, where the baseline is calculated through a moving average window having a predetermined length. In some embodiments, the predetermined length is between 50 ms and 350 ms. In some embodiments, the predetermined length is between 100 ms and 300 ms. In some embodiments, the predetermined length is between 150 ms and 250 ms. In some embodiments, the predetermined length is between 175 ms and 225 ms. In some embodiments, the predetermined length is approximately 200 ms. In some embodiments, the predetermined length is 201 ms (i.e., 50 samples at a sampling rate of 250 samples / sec). In some embodiments, the baseline includes data with frequencies below 5 Hz, and therefore, the signal is high-pass filtered at approximately 5 Hz. Preprocessed data generated based on the raw biopotential data shown in FIGS. 18A, 18C, 18E, and 18G are shown in FIGS. 18B, 18D, 18F, and 18H, respectively.

[0061] Continuing with step 1310, in some embodiments, following application of the above-described filters, each data channel is checked for contact issues. In some embodiments, contact issues are identified in each data channel based on at least one of (a) the RMS of the data channel, (b) the signal-to-noise ratio ("SNR") of the data channel, and (c) the time variation of the peak relative energy of the data channel. In some embodiments, a data channel is identified as corrupted if it has an RMS value greater than a threshold RMS value. In some embodiments, the threshold RMS value is two local voltage units (e.g., a value of approximately 16.5 millivolts). In some embodiments, the threshold RMS value is between one local voltage unit and three local voltage units. Exemplary data channels identified as corrupted by this criterion are shown in FIGS. 18A and 18B. In some embodiments, a data channel is identified as corrupted if it has an SNR value less than a threshold SNR value. In some embodiments, the threshold SNR value is 50 dB. In some embodiments, the threshold SNR value is between 40 dB and 60 dB. In some embodiments, the threshold SNR value is between 30 dB and 70 dB. Exemplary data channels identified as corrupted by this criterion are shown in FIGS. 18C and 18D. In some embodiments, a data channel is identified as corrupted if the change in relative R-wave peak energy from one interval to another is greater than a threshold change amount. In some embodiments, the threshold change amount is 250%. In some embodiments, the threshold change amount is greater than or equal to 200% and less than or equal to 300%. In some embodiments, the threshold change amount is between 150% and 350%. Exemplary data channels identified as corrupted by this criterion are shown in FIGS. 18E and 18F. Exemplary data channels that were not identified as corrupted for any of the above reasons are shown in FIGS. 18G and 18H.

[0062] In step 1315, the exemplary computing device of the present invention is programmed / configured to extract R-wave peaks from the preprocessed signal channels and generate an R-wave peak data set. In some embodiments, step 1315 uses known maternal ECG peaks as input. In some embodiments, step 1315 uses maternal ECG peaks identified according to the techniques described in U.S. Patent No. 9,392,952 as input. In some embodiments, step 1315 includes refining the maternal ECG peak locations using the preprocessed data (e.g., generated by step 1310) and the known maternal ECG peaks. In some embodiments, refining the peak locations includes searching for the maximum absolute value in a window of samples before and after the known maternal ECG peak to ensure that the R-wave peak is located at the maximum point of the R wave for each one of the filtered signals. In some embodiments, the window includes plus or minus a predetermined length of time. In some embodiments, the predetermined length is between 50 milliseconds and 350 milliseconds. In some embodiments, the predetermined length is between 100 milliseconds and 300 milliseconds. In some embodiments, the predetermined length is between 150 and 250 milliseconds. In some embodiments, the predetermined length is between 175 and 225 milliseconds. In some embodiments, the predetermined length is approximately 200 milliseconds. In some embodiments, the window comprises plus or minus a number of samples ranging between 1 and 100 samples. Illustrations of known maternal ECG peaks and extracted R-wave peaks in an exemplary R-wave peak data set are shown in FIGS. 19A and 19B, respectively.

[0063] In step 1320, an exemplary computing device of the present invention is programmed / configured to remove electromyogram ("EMG") artifacts from the preprocessed data generated by step 1310 and data including the R-wave peaks extracted in step 1315. FIG. 20A shows exemplary preprocessed data used as input to step 1320. In some embodiments, EMG artifact removal is performed to correct for high amplitude peaks with increased high frequency energy, which usually, but not always, arise from maternal EMG activity. Other sources of such energy are high power line noise and high fetal activity. In some embodiments, EMG artifact removal involves finding corrupted peaks and replacing them with a median value. In some embodiments, finding corrupted peaks involves calculating a peak-to-peak RMS value based on the following formula:

[0064] The first step in correcting this artifact is to find the corrupted peaks, which requires calculating the peak-to-peak RMS value as follows:

number

[0065] In the above equation, the peak signal is the signal having the height of the R-peak (i.e., the amplitude of the peak of the R-wave), and the peak position is the signal having the time index of the R-peak found for each channel (i.e., the time index of each of the peaks of the R-wave). In some embodiments, there are two peak signal values ​​and two peak position values, one for the R-wave peak found using the filtered data and one found using the inverse signal (i.e., the signal obtained by multiplying the original signal data by −1 to result in a sign-inverted signal).

[0066] In some embodiments, finding corrupted peaks also includes finding peaks that indicate outliers in a maternal physical activity ("MPA") dataset. In some embodiments, such a signal (hereinafter referred to as an "envelope signal") is extracted as follows:

[0067] In some embodiments, the physical activity data is collected using a motion sensor. In some embodiments, the motion sensor includes a 3-axis accelerometer and a 3-axis gyroscope. In some embodiments, the motion sensor samples 50 times per second (50 sps). In some embodiments, the sensor is located on the same sensing device (e.g., a wearable device) that includes electrodes used to collect biopotential data for performance of method 1300 as the entire device (e.g., garment 300 shown in FIG. 3).

[0068] In some embodiments, the raw motion data is transformed. In some embodiments, the raw motion data is converted to g units in the case of accelerometer raw data or to degrees per second in the case of gyroscope raw data. In some embodiments, the transformed data is examined to distinguish between valid and invalid signals by determining whether the raw signal is saturated (e.g., has a certain maximum possible value). In some embodiments, the signal envelope is extracted as follows: First, in some embodiments, the data is checked for position changes. Because position changes are characterized by an increase in the accelerometer baseline, in some embodiments, a baseline filter is applied whenever a position change occurs. In some embodiments, the filtering is performed by employing a high-pass finite impulse response (“FIR”) filter. In some embodiments, the high-pass filter has a filter order of 400 and a frequency of 1 Hz. In some embodiments, a low-pass FIR filter is also applied to reject non-physiological motion. In some embodiments, the low-pass filter has a filter order of 400 and a frequency of 12 Hz. (Order 400, fc=12 Hz[1]) is applied as well. In some embodiments, following filtering, the magnitude of the accelerometer vector is calculated according to the following formula:

number

[0069] In this formula, AccMagnitudeVector(iSample) represents the square root of the sum of the squares of the three accelerometer axes (e.g., x, y, and z) for sample number iSample. In some embodiments, the magnitude vector of the gyroscope data is calculated according to the following formula:

number

[0070] In some embodiments, the peak in the MPA motor envelope is defined according to the following steps. Motion envelope peak = find(Motion envelope > P95%(Motion envelope)) Motion Envelope Peak Onset = Motion Envelope Peak - 2 Peak Width Motion Envelope Peak Offset = Motion Envelope Peak + 2 Peak Width

[0071] In the above, the peak width is defined as the distance between the peak and the first point where the envelope reaches 50% of the peak value, and P95%(x) is 95% of x. Figure 20B shows the data signal of Figure 20A and the corresponding motion envelope and peak-to-peak absolute sum (i.e., the sum of the absolute values ​​of all samples falling between adjacent peaks) calculated as above.

[0072] In some embodiments, a peak is deemed corrupted if: 1) Peaks with peak-to-peak RMS greater than 20 local voltage units 2) If the signal inspection stage concludes that there is a contact problem in the current processing section, the peak-to-peak RMS is higher than 8 local voltage units. 3) If during the signal inspection stage it is determined that there is a contact problem in the current processing section, but more than 50% of the points have a peak-to-peak effective voltage higher than 8 local voltage units, a threshold of 20 local voltage units is used. 4) Peaks located near the onset and offset of the motion envelope are suspected of being corrupted. The peak-to-peak RMS of these points should be greater than 6 to conclude that they are corrupted.

[0073] In some embodiments, if a peak is detected as a corrupted peak as described above, the amplitude of the peak is replaced by the median and the local median around the corrupted peak is calculated as follows:

number

[0074] In some embodiments, the corrupted data points themselves are excluded from the above calculations and replaced with a statistical value (e.g., global median, local median, mean, etc.). In some embodiments, if there are seven or fewer values ​​to use after exclusion, the global median is used as the local one, and the global median is calculated using standard techniques.

number

[0075] In some embodiments, if the absolute difference between the local median and the global median is greater than 0.1, the local median is used in place of the amplitude of the corrupted data point; otherwise, the global median is used as a replacement for the amplitude of the corrupted peak. Figure 20C shows the example data set of Figures 20A and 20B with the corrupted peak replaced as described above.

[0076] Continuing with reference to FIG. 13, in step 1325, an exemplary computing device of the present invention is programmed / configured to remove baseline artifacts from signals formed by R-wave peaks. In some embodiments, such artifacts are caused by sudden baseline or RMS changes. In some embodiments, such changes are often caused by changes in maternal position. FIG. 21A shows an exemplary data signal including a baseline artifact.

[0077] In some embodiments, such artifacts are found using the Grubbs test for outliers, which is a statistical test performed based on absolute deviation from the sample mean. In some embodiments, to correct such artifacts, the point of change should first be found. In some embodiments, the point of change is the point (e.g., data point) where the signal RMS or mean begins to change, and such a point should meet the following criteria: 1) Length (peak signal) - change point > 50 2) prctile(peak signal (change point: end), 10)>0.01 3a)

number

number

[0078] In some embodiments, the change point should meet the criteria above, and the peak signal up to this point is modified based on statistics defined below.

number

[0079] FIG. 21B shows the example data signal of FIG. 21A after baseline artifact removal has been performed according to step 1330.

[0080] Continuing with reference to FIG. 13, in step 1330, an exemplary computing device of the present invention is programmed / configured to remove outliers from an R-wave peak signal using an iterative process according to the Grubbs test for outliers. FIG. 22A shows an exemplary R-wave peak signal including outlier data points, as indicated by diamonds. In some embodiments, the iterative process of step 1330 stops when either of the following two conditions occurs: 1) Outlier P 95% / P 50% Ratio>1.5 2) Number of iterations > 4

[0081] In some embodiments, this process finds outliers at each iteration and trims the height of such outliers to the median of a local region around the outlier peak. In some embodiments, the local region is defined as a time window of a predetermined number of samples before and after the outlier peak. In some embodiments, the predetermined number of samples is between 0 and 20. In some embodiments, the predetermined number of samples is 10. FIG. 22B shows the example data signal of FIG. 22A following execution of step 1330. It may be seen that the outlier data points shown in FIG. 22A are no longer present in FIG. 22B. In some embodiments, following signal extraction, more outliers are identified and removed, as described in further detail below.

[0082] Continuing with reference to FIG. 13 , in step 1335, an exemplary computing device of the present invention is programmed / configured to interpolate and extract R-wave signal data from each of the R-wave peak signal datasets to generate an R-wave signal channel. In some embodiments, the peak signals output by step 1330 are temporally interpolated to provide a 4-sample / second signal. FIG. 23A shows an exemplary peak signal output by step 1330. In some embodiments, the interpolation is performed using cubic spline interpolation. In some embodiments, in cases where large gaps in the interpolated data exist, spurious high values ​​are present, and instead the interpolation method is shape-preserving piecewise cubic interpolation. In some embodiments, the shape-preserving piecewise cubic interpolation is piecewise cubic Hermitian interpolation polynomial (“PCHIP”) interpolation. In some embodiments, following the interpolation, the step of extracting the R-wave signal includes identifying additional outliers in the interpolated signal. In some embodiments, the additional outliers are identified in this step as one of the following: 1) Signal peaks with heights greater than one local voltage unit (i.e., peaks of the interpolated R-wave signal, not peaks of the raw biopotential signal) and their surroundings 2) A point lying between two consecutive R peaks separated by more than 10 seconds 3) If a severe contact problem is discovered during the data inspection phase (e.g., during steps 1320, 1335, and 1330),

[0083] In some embodiments, points identified as outliers based on meeting any of the three criteria above are discarded and replaced with a statistical value (e.g., either the local median or the global median) according to the process described above with reference to step 1320.

[0084] Continuing with step 1335, in some embodiments, following further outlier detection, signal statistics (e.g., median, minimum, and standard deviation) are calculated and a signal (e.g., a 1-minute signal time window for a given channel) is identified as corrupted if any of the following are true: 1) Even after removing outliers, the signal still has peaks with amplitudes greater than 1 local voltage unit and standard deviations greater than 0.1. 2) The median signal is greater than 0.65 and the minimum signal is less than 0.6. 3) More than 15% of the points that make up the signal are removed as outliers.

[0085] Continuing with step 1335, following identification of the corrupted signal, a sliding RMS window is applied to the signal. In some embodiments, the RMS window has a size ranging between 25 and 200 samples. In some embodiments, the RMS window has a size of 100 samples. In some embodiments, following application of the RMS window, a first-order polynomial function is fitted to the signal and then subtracted from the signal, thereby generating a clean version of the interpolated signal, which can be used in subsequent steps. FIG. 23B shows an exemplary R-wave signal following the interpolation of step 1330.

[0086] 13 , in step 1340, an exemplary computing device of the present invention is programmed / configured to perform channel selection, whereby a subset of exemplary R-wave signal channels are selected for use in generating the electrical uterine monitoring signal. In some embodiments, at the start of channel selection, all channels are considered eligible candidates, and channels are evaluated for possible elimination according to the following: 1) Exclude channels that have had contact problems for 10% or more of the processing intervals up to the current time. 2) If more than 50% of the channels are excluded based on the above, then instead exclude all channels that have contact issues for more than 15% of the processing interval.

[0087] If the above results in all channels being eliminated, instead, keep the channels that meet both of the following criteria and eliminate the remaining channels: 1) The standard deviation of the signal is between 0 and 0.1. 2) The signal width is within 0.2.

[0088] If the above still results in all channels being excluded, then only the first of the above conditions regarding standard deviation is used and the second of the above conditions regarding range is ignored. Figure 24A shows an example data set containing six data channels, with two data channels being excluded.

[0089] In some embodiments, after removing some channels as described above, the remaining channels are grouped into couples. In some embodiments, where the channels are defined as described above, a channel couple is any pair of the eight channels described above. In some embodiments, only couples that are independent of each other (i.e., couples that do not have a common electrode) are considered. In some embodiments, the possible couples are as follows: 1. Channels 1 and 2 (A1-A4 and A2-A3) 2. Channels 1 and 5 (A1-A4 and B1-B3) 3. Channels 1 and 6 (A1-A4 and B1-B2) 4. Channels 1 and 7 (A1-A4 and B3-B2) 5. Channels 2 and 5 (A2-A3 and B1-B3) 6. Channels 2 and 6 (A2-A3 and B1-B2) 7. Channels 2 and 7 (A2-A3 and B3-B2) 8. Channels 3 and 5 (A2-A4 and B1-B3) 9. Channels 3 and 6 (A2-A4 and B1-B2) 10. Channels 3 and 7 (A2-A and B3-B2) 11. Channels 3 and 8 (A2-A4 and A1-A3) 12. Channels 4 and 5 (A4-A3 and B1-B3) 13. Channels 4 and 6 (A4-A3 and B1-B2) 14. Channels 4 and 7 (A4-A3 and B3-B2) 15. Channels 5 and 8 (B1-B3 and A1-A3) 16. Channels 6 and 8 (B1-B2 and A1-A3) 17. Channels 7 and 8 (B3-B2 and A1-A3)

[0090] As can be seen, for each of the above channel pairs, the two channels that form the pair do not share common electrodes. In some embodiments, the Kendall rank correlation for each pair of channels is calculated using only significant points within the channels. In some embodiments, the Kendall correlation counts the matching rank signs for each pair of signals to test their statistical dependence.

[0091] In some embodiments, channels are then selected according to the following selection criteria: First, if the maximum Kendall correlation value is 0.7 or greater, the selected channels are any independent channels with a Kendall correlation value of 0.7 or greater. However, if all selected channels have previously been identified as corrupted, the output signal is identified as a corrupted signal. Furthermore, if any of the selected channels have previously been identified as corrupted or if any of the selected channels have a range greater than 0.3, such channels are excluded from the selected channels.

[0092] Second, if no channels were selected under the first criterion, and the maximum Kendall correlation value is greater than or equal to 0.5 but less than 0.7, the selected channel is any independent channel with a Kendall correlation value within this range. However, if all selected channels have previously been identified as corrupted, the output signal is identified as corrupted. Furthermore, if any of the selected channels have previously been identified as corrupted, or if any of the selected channels have a range greater than 0.3, such channels are excluded from the selected channels.

[0093] Third, if no channel is selected by the first or second criteria above, and the maximum Kendall correlation value is greater than zero but less than 0.5, all channels with a Kendall correlation value greater than zero are identified as selected channels. However, if the maximum correlation value is less than 0.3, the output signal is marked as corrupted, and all channels with a range greater than 0.3 are excluded.

[0094] Fourth, if no channels were selected under the first three criteria above, then all channels with a range greater than 0.3 and all channels with a deletion point greater than 15% are eliminated, and the remaining channels are selected and their output signals are identified as those that should be less sharpened, as described below with reference to step 1355.

[0095] Fifth, all channels not selected by any of the four criteria above are selected except for channels with severe contact problems. However, if the number of contact problems in a selected channel exceeds 15, the output signal is flagged as corrupted. Figure 24B shows an example data set following the channel selection of step 1340.

[0096] In some embodiments, rather than selecting channels in pairs based on their correlation values, the channels are selected individually.

[0097] 13 , in step 1345, an exemplary computing device of the present invention is programmed / configured to calculate a uterine activity signal (which may be referred to as an “electrical uterine monitoring” or “EUM” signal) based on the selected R-wave signal channel selected in step 1340. In some embodiments, for each sample (e.g., a set of data points at a given sampling time during the sampling interval of 4 samples per second for all selected channels), the 80th percentile of the signal for the selected channel is calculated according to:

number

[0098] FIG. 25A shows the 80% signal calculated based on the selected data channels shown in FIG. 24B. In some embodiments, the wandering baseline is then removed from the combined 80th percentile signal, as determined above, to generate the EUM signal. In some embodiments, a moving average window is considered to find the baseline. In some embodiments, the moving average window subtracts the average value during the window from the EUM signal. In some embodiments, the window length is between 0 and 20 minutes. In some embodiments, the window length is 10 minutes. FIG. 25B shows the example signal of FIG. 25A following baseline removal.

[0099] In step 1350, an exemplary inventive computing system is programmed / configured to normalize the EUM signal calculated in step 1345. In some embodiments, normalization comprises multiplying the EUM signal from step 1345 by a constant. In some embodiments, the constant is between 200 and 500. In some embodiments, the constant is between 250 and 450. In some embodiments, the constant is between 300 and 400. In some embodiments, the constant is between 325 and 375. In some embodiments, the constant is approximately 350. In some embodiments, the constant is 350. In some embodiments, the constant is 1, i.e., the original value of the extracted 80th percentile signal is maintained. FIG. 26 shows an exemplary data signal following normalization of the data signal of FIG. 25B according to step 1350.

[0100] In step 1355, the exemplary computing system of the present invention is programmed / configured to sharpen the normalized EUM signal generated by step 1350, thereby generating a sharpened EUM signal. In some embodiments, sharpening is performed only on signals that were not flagged as corrupted in the preceding steps; if all relevant signals are flagged as corrupted, the sharpening step is not performed. In some embodiments, the purpose of the sharpening step is to enhance all areas suspected of contraction. In some embodiments, sharpening proceeds as follows: First, if there is an EUM signal peak that exceeds a value of 200 local voltage units, the signal is marked as corrupted. Second, it is determined whether the signal was previously marked as corrupted. Third, the signal baseline is removed. In some embodiments, for baseline removal, if the signal duration is greater than 10 minutes, a 10-minute long moving average window is used to estimate the baseline; otherwise, the 10th percentile of the signal is used to estimate the baseline. In either case, the baseline is then subtracted from the EUM signal. Fourth, the signal baseline is defined as 30 visual voltage units. In some embodiments, the signal baseline thus defined following a normalization step provides an EUM signal that is in the range of 0 to 100 in a manner similar to the signal provided by an electrocardiograph.

[0101] Fifth, peaks are identified according to one of the following: If the signal is identified as not requiring much sharpening during step 1340, a peak is defined as having a minimum height of 35 visual voltage units and a minimum width of 300 samples. If the signal was not so identified, the peak was identified as having a minimum height of 35 visualization units and a minimum width of 220 samples.

[0102] In either case, the prominence of each peak is calculated according to the following formula:

number

[0103] After calculating the prominence of all peaks in the sample, peaks that fall into any of the following categories are excluded. The peak prominence is 12 or less and the height is 40 or less visual voltage units. The peak has a prominence that is less than 65% of the maximum prominence of all peaks in the sample.

[0104] In some embodiments, additional peaks are identified by identifying any further peaks (e.g., local maxima) that have a minimum height of 15 visual voltage units and a minimum width of 200 samples, and then removing all peaks that have a prominence greater than 20 visual voltage units.

[0105] In accordance with the above, sharpening is performed only if all of the following are true: (a) the signal is not corrupted (a "corrupted" signal is identified as above), (b) there are no deletion points in the signal, and (c) at least one peak was identified in the preceding part of this step. If sharpening is to be performed, peaks that meet any of the following conditions are removed prior to sharpening: Peaks have a prominence of less than 10 visible voltage units. The peak has a prominence of 35 or more visual voltage units. The peak width is 800 samples or more (200 seconds at 4 samples per second).

[0106] Peaks that meet the above conditions are excluded, and the following values ​​are calculated for each remaining peak.

number

number

number

[0107] After calculating these values, a mask is created with zero values ​​outside the peak and a Gaussian function inside the peak according to the following formula:

number

[0108] The mask is then smoothed with a moving average window having a predefined length. In some embodiments, the predefined length is between 10 and 50 seconds. In some embodiments, the predefined length is between 20 and 40 seconds. In some embodiments, the predefined length is between 25 and 35 seconds. In some embodiments, the predefined length is approximately 30 seconds. In some embodiments, the predefined length is 30 seconds. An exemplary EUM signal is shown in FIG. 27A, and an exemplary mask created in the above manner for the exemplary EUM signal of FIG. 27A is shown in FIG. 27B. The mask is then added to the existing EUM signal to generate a sharpened EUM signal. In some embodiments, the addition is performed using simple mathematical addition. An exemplary sharpened EUM signal created by adding the exemplary mask of FIG. 27B to the exemplary EUM signal of FIG. 27A is shown in FIG. 27C.

[0109] Returning to FIG. 13 , post-processing is performed in step 1360 to generate a post-processed EUM signal. In some embodiments, post-processing includes baseline removal. In some embodiments, baseline removal includes removing the signal baseline, as described above with reference to step 1355. In some embodiments, for baseline removal, if the signal duration is greater than 10 minutes, a 10-minute long moving average window is used to estimate the baseline; otherwise, 10% of the signal value is used to estimate the baseline. In either case, the baseline is then subtracted from the EUM signal, and the signal baseline is defined as 30 visual voltage units. Finally, all removed values ​​are set to a value of −1 visual voltage unit, and all values ​​greater than 100 visual voltage units are set to a value of 100 visual voltage units. FIG. 28 shows an exemplary post-processed signal generated by applying the post-processing of step 1360 to the exemplary sharpened signal of FIG. 27B .

[0110] Following step 1360, method 1300 is complete. As noted above, FIG. 28 illustrates an exemplary EUM signal calculated according to method 1300. FIG. 29 illustrates a representative tocograph signal obtained according to known techniques for the same subject during the same time period as the collection of biopotential data from which the EUM signal of FIG. 28 was calculated. It will be appreciated that FIGS. 28 and 29 are substantially similar to one another and include the same peaks that may be understood to represent contractions. Thus, the result of method 1300 may be appreciated as an EUM signal that can be used as a tocograph-like signal for monitoring maternal uterine activity, but that can be calculated based on non-invasively recorded biopotential signals.

[0111] Reference is now made to the following examples, which together with the above description illustrate, by way of non-limiting example, some embodiments of the present invention.

[0112] 14A-17B show further examples of comparisons between tochograph data and the output of exemplary method 200. In each of FIGS. 14A, 15A, 16A, and 17A, the tochograph signal is plotted against time, with the mother's self-reported contractions monitored by the tochograph indicated by vertical lines. In each of FIGS. 14B, 15B, 16B, and 17B, the filtered R-wave signals from each of the multiple channels are shown in different colors (e.g., similar to the plot shown in FIG. 11B), and the calculated normalized average signal is shown as a thick black line (e.g., similar to the plot shown in FIG. 12A). Each of FIGS. 14B, 15B, 16B, and 17B is shown adjacent to its counterpart in FIGS. 14A, 15A, 16A, and 17A for comparison (e.g., FIGS. 14A and 14B show different data recorded on the same mother over the same time interval, and FIGS. 15A-17B show such data). As discussed above with reference to FIGS. 12A and 12B, the peaks in the exemplary normalized uterine signal are seen to correspond to self-reported contractions.

[0113] A study was conducted to evaluate the effectiveness of exemplary embodiments. The study involved subjects with a BMI of 45 kg / m 2 This study compared EUM and TOCO recordings in pregnant women aged 18–50 years, carrying singleton fetuses at gestational ages >32+0 weeks, and without fetal abnormalities. EUM was calculated as described above for data samples measured for a minimum of 30 minutes. Analysis of a maternal myocardial R-wave amplitude-based uterine activity index, referred to herein as EUM, showed promising results as an innovative and reliable method for monitoring maternal uterine activity. EUM data were highly correlated with TOCO data. Therefore, EUM monitoring may provide useful data similar to TOCO data while overcoming the drawbacks of conventional tocodynamometry, such as discomfort.

[0114] 18A-27B show exemplary data present at various stages during the execution of exemplary method 1300. In particular, Figures 27A and 27B show a comparison of the output signal generated by exemplary method 1300 with a tochograph signal recorded during the same time interval.

[0115] 18A-18H illustrate exemplary raw data received as input to exemplary method 1300 (e.g., received at step 1305) and exemplary filtered raw data generated during exemplary method 1300 (e.g., generated by step 1310). In particular, FIGS. 18A, 18C, 18E, and 18G illustrate exemplary raw data, and FIGS. 18B, 18D, 18F, and 18H illustrate exemplary filtered data, respectively. While FIGS. 18A-18H represent single-channel raw and filtered biopotential data, those skilled in the art will appreciate that in an actual implementation of the above-described method 1300, data sets equivalent to those shown in FIGS. 18A-18H will exist for each channel of data. With reference to FIG. 18A, it may be seen that there is high power line noise around sample number 6000. With reference to FIG. 18B, it may be seen that the power line noise remains high. In some embodiments, this may result in this interval being flagged as having a severe contact problem due to a change in relative R-wave peak energy from one interval to another that is greater than the threshold described above with reference to step 1310 of exemplary method 1300. With reference to FIG. 18C, it may be seen that there is high power line noise around sample number 14000. With reference to FIG. 18D, it may be seen that the power line noise is still high. In some embodiments, this may result in this interval being flagged as having a severe contact problem because the signal RMS exceeds the threshold described above with reference to step 1310 of exemplary method 1300. With reference to FIG. 18E, it may be seen that there is high power line noise throughout the signal. With reference to FIG. 18F, it may be seen that the power line noise is still high. In some embodiments, this may result in this interval being flagged as having a severe contact problem because the SNR of this signal does not meet the threshold SNR described above with reference to step 1310 of exemplary method 1300. With reference to FIGS. 18G and 18H, a clear signal may be visible.In some embodiments, this may result in this interval not being flagged as having a contact problem.

[0116] 19A and 19B, the extraction of R-wave peaks according to step 1315 is illustrated. While FIGS. 19A and 19B illustrate R-wave peak extraction from a single channel, those skilled in the art will appreciate that in an actual implementation of the method 1300 described above, data sets equivalent to those shown in FIGS. 19A and 19B would exist for each channel of data. FIG. 19A illustrates filtered data (e.g., generated by step 1310) prior to execution of step 1315. In FIG. 19A, detected peak locations are represented by asterisks. FIG. 19B illustrates data with extracted peaks following execution of step 1315. In FIG. 19B, the peak locations are represented by asterisks. It may be noted that some of the peak locations indicated by asterisks in FIG. 19A are not located at the maxima of the peaks in the data, and such locations are correctly indicated by asterisks in FIG. 19B.

[0117] Referring now to FIGS. 20A-20C, the removal of EMG artifacts according to step 1320 is illustrated. While FIGS. 20A-20C illustrate the removal of EMG artifacts from a single channel, those skilled in the art will appreciate that in an actual implementation of the method 1300 described above, data sets equivalent to those shown in FIGS. 20A-20C would exist for each channel of data. FIG. 20A illustrates exemplary filtered data used in step 1320 (e.g., generated by step 1310). FIG. 20B illustrates the same filtered data of FIG. 20A, further including a representation of the motion envelope and peak-to-peak absolute sum. In FIG. 20B, suspected corrupted peaks are indicated by diamonds. FIG. 20C illustrates the corrected signal after EMG artifact correction, as generated by step 1320. In FIG. 20C, the suspected peaks have been removed, the corrected peaks are indicated by circles, and the original peak values ​​are shown in contrasting shades.

[0118] 21A and 21B, baseline artifact removal according to step 1325 is illustrated. While FIGS. 21A and 21B depict baseline artifact removal from a single channel, it will be apparent to one skilled in the art that in an actual implementation of the above-described method 1300, a data set equivalent to that shown in FIGS. 21A and 21B would exist for each channel of data. FIG. 21A illustrates exemplary data before baseline artifact removal that may be received as input to step 1325. In FIG. 21A, the baseline artifact is shown within a circle. In the data illustrated in FIG. 21A, the baseline ratio between the circled region and the remainder of the signal is less than 0.8. In some embodiments, dividing the remainder of the signal by this factor provides a corrected signal. FIG. 21B illustrates an exemplary corrected signal as may be produced by step 1325. In FIG. 21A, the baseline artifact region is shown within a circle. By comparing FIGS. 21A and 21B, it may be seen that the baseline artifact has been removed.

[0119] 22A and 22B, there is shown the trimming of outliers and gaps in accordance with step 1330. While FIGS. 22A and 22B illustrate the trimming of outliers and gaps from a single channel, it will be apparent to those skilled in the art that in an actual implementation of the method 1300 described above, a data set equivalent to that shown in FIGS. 22A and 22B would exist for each channel of data. FIG. 22A shows exemplary data that may be received as input to step 1330. It may be seen that the input data includes an outlier near sample 450, indicated by a diamond in FIG. 22A. FIG. 22B shows the exemplary data of FIG. 22A after execution of step 1330 to remove the outliers as described above. It may be seen that the outlier shown in FIG. 22A has been removed.

[0120] 23A and 23B, the interpolation and extraction of an R-wave peak signal in accordance with step 1330 is illustrated. While FIGS. 23A and 23B depict the extraction of an R-wave peak signal from a single channel, those skilled in the art will appreciate that in an actual implementation of the above-described method 1300, a data set equivalent to that shown in FIGS. 23A and 23B would exist for each channel of data. FIG. 23A illustrates an exemplary R-wave peak signal that may be provided as an output from step 1330 and received as an input to step 1335. FIG. 23B illustrates an exemplary clean, interpolated R-wave signal that may be produced by performance of step 1335.

[0121] 24A and 24B, channel selection according to step 1335 is shown. In the exemplary data set shown in FIGS. 24A and 24B, channels 3 and 8 were found to be ineligible for channel selection due to the presence of contact problems for more than 10% of the time interval. Therefore, only exemplary channels 1, 2, 4, 5, 6, and 7 are shown in FIGS. 24A and 24B. The independent channel pairs and corresponding Kendall correlation values ​​for the data shown in FIG. 24A are shown in the table below. [Table 1]

[0122] From the table above, it may be seen that the group consisting of channels 1, 2, 4, and 7 exhibits a moderate correlation (e.g., a correlation greater than 0.5 and less than 0.7). Therefore, in step 1340, channels 1, 2, 4, and 7 are selected. Figure 24B shows an exemplary data set output by step 1340, including selected channels 1, 2, 4, and 7.

[0123] 25A and 25B, there is shown the calculation of an EUM signal based on selected channels according to step 1345. The channel data shown in FIG. 24B is received as input to step 1345 to generate the output data shown in FIGS. 25A-25B. Referring to FIG. 25A, this shows the 80th percentile signal extracted from the signal shown in FIG. 24B. FIG. 25B shows the corrected signal obtained by applying wandering baseline removal to the signal shown in FIG. 25A.

[0124] Referring now to Figure 26, there is shown the calculation of a normalized EUM signal according to step 1350. The corrected data produced by step 1345 and shown in Figure 25B is received as input to step 1350 to generate a normalized EUM signal as shown in Figure 26. Figure 26 shows the normalized signal obtained by normalizing the signal shown in Figure 25B and setting the reference value to 30 visual voltage units. In Figure 26, it can be seen that there are three weak peaks in the signal.

[0125] 27A-27C, sharpening of an EUM signal according to step 1355 is shown. An exemplary normalized signal produced by step 1350, such as the exemplary normalized signal shown in FIG. 26, is received as input to step 1355 to produce a sharpened EUM signal. FIG. 27A shows an exemplary normalized EUM signal as produced by step 1350. FIG. 27B shows an exemplary enhancement mask produced according to step 1355. FIG. 27C shows an exemplary sharpened EUM signal produced by applying the normalized EUM signal of FIG. 27A to the mask of FIG. 27B.

[0126] Referring now to Figure 28, post-processing of the EUM signal according to step 1360 is shown. The sharpened EUM signal, such as that produced by step 1355, is received as input to step 1360 to produce a post-processed EUM signal. Figure 28 shows an exemplary post-processed EUM signal after removing the wandering baseline as described above with reference to step 1360. It may be seen that the three weak peaks shown in Figure 26 are more clearly visible in Figure 28 following the sharpening of step 1355 and the post-processing of step 1360.

[0127] Referring now to Figure 29, there is shown a tochograph signal corresponding to the exemplary EUM signal of Figure 28. As previously mentioned, the exemplary EUM signal of Figure 28 is generated according to method 1300. The exemplary tochograph signal of Figure 29 was acquired for the same subject during the same time interval as the data used to generate the exemplary EUM signal of Figure 28. It will be seen that Figures 28 and 29 are substantially consistent with each other and contain the same three peaks as each other.

[0128] In some embodiments, uterine monitoring is performed based on acoustic data collected using one or more acoustic sensors, such as acoustic sensor 320 described above with reference to FIG. 3 . In some embodiments, the process for uterine monitoring based on acoustic data is substantially similar to the process for uterine monitoring based on biopotential data described above with reference to method 1300 shown in FIG. 13 , except as described below. FIG. 30 shows an exemplary method 3000 for uterine monitoring based on acoustic data. In some embodiments, an exemplary computing device of the present invention is programmed / configured to perform method 3000. In some embodiments, the exemplary computing device of the present invention is programmed / configured according to method 3000 via instructions stored on a non-transitory computer-readable medium. In some embodiments, the exemplary computing device of the present invention includes at least one computer processor, which, when executing the instructions, becomes a specific computer processor programmed / configured according to method 3000. In some embodiments, the exemplary computing device of the present invention is specifically configured to solve the technical problems described below by performing method 3000.

[0129] In step 3005, the exemplary inventive computing device is specifically configured to receive raw acoustic data as input. In some embodiments, a set of raw acoustic data is received from each of a plurality of acoustic sensors positioned proximate the abdomen of the pregnant human subject. In some embodiments, the set of raw acoustic data is received from each of two, three, four, five, six, seven, eight, nine, ten, or more acoustic sensors. In one specific exemplary embodiment detailed in this description of method 3000, a set of raw acoustic data is received from each of four acoustic sensors, as shown in FIG. 3 .

[0130] In step 3010, an exemplary inventive computing device is specifically configured to preprocess the raw acoustic data to generate multiple channels of preprocessed acoustic data. In some embodiments, the preprocessing includes applying at least one filter (e.g., one filter, or two filters, or three filters, or four filters, or five filters, or six filters, or seven filters, or eight filters, or nine filters, or ten filters, or more) to the raw acoustic data, e.g., applying an amount X of filters to each of an amount Y of channels of raw data to generate an amount X times Y of preprocessed data channels. In some embodiments, the filter includes a bandpass filter. In some embodiments, the filter includes a DC filter. In some embodiments, the filter includes a finite impulse response filter, or an infinite impulse response (“IIR”) filter such as a Butterworth filter or a Chebyshev filter, or a combination thereof. In some embodiments, the filter includes a lowpass zero-phase lag IIR filter with a 50 Hz cutoff. In some embodiments, the filter includes a 12th-order Butterworth IIR filter, or a 3rd-order Butterworth IIR filter, or a 5th-order Butterworth IIR filter. In one exemplary embodiment, the filter includes five 12th-order Butterworth IIR filters having frequencies: 10-50 Hz, 15-50 Hz, 20-50 Hz, 25-50 Hz, and 30-50 Hz. In some embodiments, application of the five IIR filters to the four raw data channels produces 20 preprocessed data channels. Figure 31A shows exemplary preprocessed data channel data following step 3010. Figure 31B shows an expanded view of a small time window of the data shown in Figure 31A.

[0131] In step 3015, the exemplary computing device of the present invention is specifically configured to extract S1-S2 peaks from the pre-processed data channel. Those skilled in the art will recognize that S1 and S2 refer to the first and second sounds not associated with a cardiac cycle. In some embodiments, the term "S1-S2 peak" refers to the maximum point within a given S1-S2 complex. In some embodiments, the S1-S2 peak extraction of step 3015 is performed in a manner substantially similar to the R-wave peak extraction of step 1315 of method 1300 described above. FIG. 32A illustrates data from an exemplary data channel, showing the R-wave peaks with annotations. FIG. 32B is a close-up of a small time window of the data illustrated in FIG. 32A. FIG. 32C illustrates an exemplary S1-S2 amplitude signal determined based on the R-wave peaks as illustrated in FIG. 32A. FIG. 32D illustrates an exemplary R-wave amplitude signal over a larger time window.

[0132] In steps 3020, 3025, and 3030, the exemplary inventive computing device is specifically configured to remove artifacts and outliers from the dataset generated in step 3015 in the same manner as described above with reference to steps 1320, 1325, and 1330 of method 1300. Note that the acoustic data analyzed by exemplary method 3000 may not include electrical noise of the type described above with reference to step 1320, but rather may include movement-related noise that is typically recorded by an acoustic sensor. However, the process for removing such movement-related noise is substantially similar to the process for removing electrical noise described above. Figure 33 shows an exemplary dataset of an exemplary data channel following execution of steps 3020, 3025, and 3030.

[0133] In step 3035, the exemplary computing device of the present invention is specifically configured to interpolate and extract S1-S2 signal data from the data sets generated in step 3030 in a manner substantially similar to that described above with reference to step 1335 of method 1300. Figure 34 illustrates the extracted S1-S2 data sets for multiple channels calculated in step 3035.

[0134] In step 3040, the exemplary computing device of the present invention is specifically configured to perform channel selection in a manner substantially similar to that described above with reference to step 1340 of method 1300. However, the channel selection in step 3040 differs in one aspect from the channel selection in step 1340. As described above, some of the data channels used in step 1340 are not independent of each other due to the differential nature of the biopotential sensors, and as a result, only some of the data channels used in step 1340 may be combined with each other. In contrast, the acoustic sensors collecting the data used in method 3000 are single-ended, i.e., independent of each other. As a result, any two channels of data may be appropriately combined with each other in step 3040. Thus, for example, in an embodiment in which four raw data channels are processed with five different bandpass filters to generate 20 filtered data channels, there are 20 x 19, or 380, possible channel couples.

[0135] Following channel selection in step 3040, in step 3045, the exemplary computing device of the present invention is specifically configured to calculate an acoustic uterine activity signal in a manner substantially similar to that described above with reference to step 1345 of method 1300. In step 3050, the exemplary computing device of the present invention is specifically configured to normalize the acoustic uterine activity signal in a manner substantially similar to that described above with reference to step 1350 of method 1300. In step 3055, the exemplary computing device of the present invention is specifically configured to sharpen the normalized acoustic uterine activity signal in a manner substantially similar to that described above with reference to step 1355 of method 1300. In step 3060, the exemplary computing device of the present invention is specifically configured to post-process the sharpened acoustic uterine activity signal in a manner substantially similar to that described above with reference to step 1360 of method 1300.

[0136] In some embodiments, the output of exemplary method 3000 is an acoustic uterine monitoring signal determined non-invasively through analysis of data obtainable by acoustic sensors placed around the abdomen of a pregnant human subject. In some embodiments, the acoustic uterine monitoring signal produced by exemplary method 3000 provides uterine monitoring data similar to that produced by tocodynamometers and ultrasound transducers and can be used to monitor uterine activity, such as contractions.

[0137] Figures 35A-37B show examples of comparisons between tocograph data and the output of exemplary method 3000. In each of Figures 35A, 36A, and 37A, the tocograph signal is shown versus time. In each of Figures 35B, 36B, and 37B, the output of exemplary method 3000 using acoustic data recorded during the same time interval is shown. It can be seen that peaks in the exemplary acoustic-based uterine monitoring signal correspond to peaks in the tocograph data.

[0138] FIG. 38 illustrates a set of ECG-based EUM processed signals and PCG-based processed signals collected from biopotential sensors and acoustic sensors, respectively. In some embodiments, a uterine monitoring signal can be determined based on data collected from both biopotential sensors (e.g., as described above with reference to the methods shown in FIGS. 2 and 13 ) and acoustic sensors (e.g., as described above with reference to the method shown in FIG. 30 ). Example data collected from biopotential sensors or ECG-based processed EUM signals is shown in section 3801 (the first two rows represent eight channels), and example data collected from acoustic sensors or PCG-based processed signals is shown in 3802 (the last five rows represent twenty channels). Each column in section 3802 represents an acoustic channel, and each row in section 3802 represents one of five filters in a filter bank. Signals such as those shown in FIG. 38 can be combined or fused to generate a single uterine activity signal.

[0139] 39 shows a method 3900 for performing fusion processing to generate a uterine activity signal from an ECG-based processed EUM signal and a PCG-based processed signal. In some embodiments, the ECG-based signal can be received in parallel from N channels as shown at 3901, or sequentially from M channels at 3903. Machine learning techniques can then be performed at 3905 to determine optimal weights for each of channels M and N. The optimal weights determined at 3905 can be used to combine the received ECG-based processed EUM signal and the received PCG-based signal into a signal representative of uterine activity. Such a signal can be generated, for example, as a weighted average of the N and M channels, where each channel is associated with a signal.

[0140] In some embodiments, the machine learning channel weighting technique can be implemented as a gradient descent (GD) optimization process. For example, a weighting value can be assigned to each of the 28 channels depicted in Figure 38. The cost function for the gradient descent process can be defined as follows:

[0141] At each iteration of the gradient descent optimization process, a final signal output is determined based on the weighted average assigned to each of the 28 channels. Contractions are identified for these final signal outputs using a change-from-baseline detection algorithm, which defines the start and end time points of each contraction in the signal. For each identified contraction, a set of features is calculated. Such a feature set may include contraction rise time, contraction fall time, the ratio between contraction rise time and contraction fall time, SNR, contraction skewness, and other suitable features. The mean value of each feature is then calculated across all contractions in the final signal. An optimal target value is then determined for each feature (e.g., based on a standardized contraction dataset or an optimal contraction dataset). The cost function of the GD process may correspond to the difference between the optimal target value and the mean of the feature value.

[0142] In some embodiments, multiple instances of the gradient descent optimization process can be run simultaneously with different initial weights assigned to each channel. For example, in a first instance, all channels can be assigned equal or the same value. In a second example, weights can be assigned to channels based on the quality of the contraction features detected by such channels. For example, contractions can be identified through each channel, and for each contraction, a set of features, such as contraction rise time, contraction fall time, the ratio of contraction rise time to contraction fall time, SNR, contraction skewness, and other suitable features can be calculated. An average feature value can be calculated across all contractions identified in a channel. A weight inversely proportional to the difference between the average feature value and the optimal feature value can then be assigned to the channel. In a third example, a clustering algorithm can be used to assign weights to channels. For example, for each channel, the correlation between the SNR and all other channels can be determined. Clusters of channels can be defined according to their SNR and correlation with other channels. Each cluster can then be combined into a single channel, and each combined channel can be given a weight based on the quality of the contraction features detected by that channel. In some embodiments, an optimal result may be selected from the first, second, and third instances of the gradient descent optimization process described above. As noted above, in step 3907, a final output, i.e., a final uterine activity signal, is determined based on a weighted average of all selected channels. FIG. 44 illustrates a set of exemplary electrical uterine monitoring signals 4410 received according to step 3901, a set of exemplary acoustic uterine monitoring signals 4420 received according to step 3903, weights 4430 assigned to each of the signals according to step 3905 (e.g., determined by method 4000 described below (note that for clarity, only a portion of the weights 4430 are specifically pointed out in FIG. 44 )), and an exemplary final uterine activity signal 4440 determined according to step 3907.

[0143] In some embodiments, the process for determining channel weights (e.g., the process of step 3905) is performed according to exemplary method 4000 shown in FIG. 40. In step 4010, a set of electrical uterine activity signals and a set of acoustic-based / PCG-based uterine activity signals are received. In some embodiments, the set of electrical uterine activity signals is generated according to step 1335 of method 1300 shown in FIG. 13. In some embodiments, the set of acoustic uterine activity signals is generated according to step 3035 of method 3000 shown in FIG. 30.

[0144] In step 4020, multiple channel sets are initialized. In some embodiments, each channel set includes a different combination of channels (possibly overlapping with one another). In some embodiments, in each set of weights, certain channels are assigned a non-zero weight, while other channels are "zeroed out." For example, if one set of weights is defined as only a set of biopotential channels (e.g., electrical uterine activity channels), all channels originating from acoustic data are assigned a weight of 0, and channels originating from biopotential data are assigned weights based on their quality, as described in more detail below. After the initial selection of weights for each set, an optimization stage is performed using gradient descent and reinforcement, as described in more detail below with reference to subsequent steps of method 4000, and at the end of the process, an optimal set of optimized weights is selected and the data is weighted and averaged according to the selected set of weights. In some embodiments, the multiple weight sets initialized in step 4020 include four weight sets. In some embodiments, the multiple channel sets include:

[0145] 1. A "biological set" consisting of data issued only from biopotential channels.

[0146] 2. An "acoustic set" consisting of data emitted only from the acoustic channel.

[0147] 3. The "shrinkage basis set" that selects channels based on K-means clustering of the shrinkage features, which is described below.

[0148] 4. The "joined set" where all channels are considered.

[0149] In some embodiments, for the first two sets, non-zero weights are assigned to the channels based on data type, as described above. In some embodiments, for the contracted basis set, the initial weights are determined according to method 4100 described below with reference to FIG.

[0150] 41 is a flowchart of a method 4100 for initializing a contraction base set of channels. In step 4110, method 4100 takes as input a set of electrical uterine monitoring channels and a set of acoustic uterine monitoring channels, as described above with reference to step 4010 of method 4000. In step 4120, contractions are identified in each of the channels. In some embodiments, contractions are identified according to method 4200, described below with reference to FIG. 42.

[0151] In step 4130, multiple contraction features are determined for each channel. In some embodiments, six contraction features are determined for each channel. In some embodiments, the features determined for each channel include:

[0152] 1. Kurtosis of the signal during contraction (averaged across the contraction).

[0153] 2. Relative energy: the ratio of the sum of all values ​​during a contraction (sig(conts)) to the sum of all channel values, per channel:

number

[0154] 3. Relative Time: The combined time of all contractions divided by the time of the entire channel data.

[0155] 4. Differential energy: The ratio of the effective value of the first derivative of the signal during contraction to the effective value of the first derivative of the entire signal.

[0156] 5. Time skew: The ratio of the mean rise time to the mean fall time, calculated for each contraction as the difference between the onset and peak contraction amplitude and the peak amplitude and offset, respectively.

[0157] 6. Contraction SNR: Calculated as the average between two determined SNRs (global SNR and averaged contraction SNR). Global SNR is equal to the rms value of the derivatives of all contractile activity for a given channel divided by the rms value of the derivatives of all signals located outside the contraction. Averaged contraction SNR is equal to the averaged SNR over individual contractions, given by the RMS of the derivatives of contractile activity divided by the RMS of the derivatives of activity located immediately surrounding the particular contraction.

[0158] In some embodiments, the output of step 4130 is a feature matrix of size Nx6, where N is the number of channels considered and 6 is the number of features determined for each channel. In step 4140, the channels are clustered. In some embodiments, the channels are clustered by performing K-means clustering on the feature matrix output by step 4130. In some embodiments, different types of clustering methods (e.g., K-medox rastering, hierarchical clustering, etc.) are used to perform clustering on the feature matrix output by step 4130. In some embodiments, the cluster with the highest number of maxima across the features is then retained as the "best cluster."

[0159] In step 4150, the best cluster is refined. In some embodiments, in an iterative process, the best cluster of channels is refined by removing channels that may reduce the internal agreement between channels in the cluster. In some embodiments, for that purpose, an internal correlation matrix of all channel pairs in the cluster is calculated. In some embodiments, linear correlation (e.g., Pearson correlation) is applied to calculate the internal correlation matrix. In some embodiments, another correlation method is applied. In some embodiments, a candidate channel with the lowest correlation with others is preliminarily removed from the cluster, and the internal correlation matrix is ​​recalculated. In some embodiments, if the preliminarily removal of the candidate channel results in an improvement in the internal correlation beyond a predefined threshold, the candidate channel is removed. In some embodiments, as a second step, all channels that are not part of the best cluster are tested for cross-correlation with the averaged cluster signal, and those with sufficiently high cross-correlation and a small lag are added to the cluster. In some embodiments, the lag is calculated using a cross-correlation function (which provides the cross-correlation coefficient and lag as a one-dimensional array), and the final cross-correlation coefficient is taken as the maximum of the calculated cross-correlation coefficients, and the lag is taken as the lag value corresponding to the same array element for which the final correlation coefficient is calculated. In some embodiments, this calculation can be expressed in pseudocode as follows:

[0160] corr_coefs, lags = cross_correlation(signal1, signal2)

[0161] corr_coef, ind_of_corr_coef = max( corr_coefs )

[0162] lag = lags[ind_of_corr_coef]

[0163] In some embodiments, the cross-correlation is sufficiently high if ρ>0.8. In some embodiments, the lag is sufficiently low if the lag is less than a maximum lag threshold. In some embodiments, the maximum lag threshold is between 30 and 60 seconds, or between 30 and 40 seconds, or between 30 and 50 seconds, or between 40 and 60 seconds, or between 40 and 50 seconds, or between 50 and 60 seconds, or less than 60 seconds, or between 25 and 35 seconds, or approximately 30 seconds, or 30 seconds. In some embodiments, the result of step 4150 is a cluster of channels with very high inter-correlation for processing during the weight selection process below.

[0164] In step 4160, regions of excluded channels (e.g., channels not selected for inclusion in the best cluster after step 4150) are considered for inclusion in the best cluster if such regions belong to contractions with high quality. In some embodiments, this "regional" data inclusion is performed by finding suitable data points (e.g., data points with good contractile activity) in the non-included channels, zeroing out all other data points in the channel, and including those "processed" channels with good regions as part of the best cluster as well. In some embodiments, good regions in otherwise unincluded channels are identified as follows: If the contraction SNR of a particular channel (e.g., Feature #6 for each channel described above) exceeds a threshold, the contractions of that channel are tested for their individual SNR against the associated threshold. In some embodiments, both the SNR and the threshold are unitless values. In some embodiments, the threshold is any value greater than zero. In some embodiments, the threshold is between 0.1 and 5, or between 0.1 and 4, or between 0.1 and 3, or between 0.1 and 2, or between 0.1 and 1, or 0.5, or 1, or 1.5, or 2, or 2.5, or 3, or 3.5, or 4, or 4.5, or 5. In some embodiments, data from high SNR contractions are retained if they come from timepoints during which there is no contractile activity in the previously included "best cluster" channel and if they are equal to or greater than a minimum length. In some embodiments, the minimum length is between 1 and 9 minutes, or between 2 and 8 minutes, or between 3 and 7 minutes, or between 4 and 6 minutes, or about 5 minutes, or 5 minutes. In other words, contraction data is added to the pool of good channel data if the SNR is equal to or greater than the threshold and includes new timepoints that were not part of a previously included contraction, and these new timepoints are not too rare. In some embodiments, all other timepoints for such channels are then zeroed, and the remaining channels with "good" contractile activity are added to the best cluster.The best clusters are output in step 4170 for use as a reduced basis set for the channels to assign weights to.

[0165] Returning to Figure 40, in step 4030, a set of initial weights is defined for each of the channel sets defined in step 4020. In some embodiments, two weight subsets are defined for each channel. In some embodiments, the two weight subsets include: (1) a "channel voting" subset, as described below, and (2) an "equal by nature" subset, in which each channel in the channel set has an initial weight of 1 / N, where N is equal to the number of channels in the set.

[0166] In some embodiments, the initial weights in the channel voting subsets are determined as follows: The voting used to create the initial subset of weights in each set is a process in which, for each data point in a channel, the number of other channels that either feature a contraction or do not feature a contraction at that same data point is counted. In other words, every channel "votes" for each data point on the type of activity (e.g., contraction or not contraction) in every other channel. The votes across data points are then averaged, and a vote count metric is calculated for each specific channel, reflecting the degree of agreement across the entire set of channels with the contraction identified in the specific channel. The voting process is repeated, with each channel being voted for by every other channel.

[0167] In some embodiments, in addition to the votes, an average of two shrinkage scores is calculated for each channel, where the shrinkage scores are determined according to step 4280 of method 4200 described below. An initial weight is then calculated for each channel as the sum of the following three items:

[0168] 1. The ratio of the number of votes to the number of channels in the set.

[0169] 2. The ratio between the two contraction scores and a pre-set score threshold.

[0170] 3. Reliability of contractions, calculated as the ratio between the sum of signal values ​​during identified contractions and the sum of the entire signal, divided by the number of contractions.

[0171] The resulting weights are then normalized so that the sum of the weights across channels equals 1. In some embodiments, the uterine monitoring process analyzes the received data in a set of "frames," processes the data over the course of a given frame, and provides an output (e.g., a uterine monitoring signal) at the end of the frame. In some embodiments, a frame is 10 minutes long. In some embodiments, for any recorded frame that is not the first (e.g., beyond the first 10 minutes while a given patient is being monitored), the weights are averaged with the weights of the previous segment to mitigate sudden weight changes between processing segments. In some embodiments, the weighting for each channel is determined as described above as 0.6 times the channel's previous weight plus 0.4 times the channel's calculated current weight.

[0172] In step 4040, the weights are optimized. In some embodiments, the optimization is performed using a gradient descent process. In some embodiments, the gradient descent algorithm adjusts the weights by attempting to minimize a cost function in an iterative process. In some embodiments, the iterative process has a configurable maximum number of iterations. In some embodiments, the iterative process has a maximum of 20 iterations. In some embodiments, the iterative process has a maximum number of iterations of 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, or 19. In some embodiments, the cost function is calculated for each optimization iteration as follows: The signals of the relevant channels in a given set are averaged into a single time series using current weights (e.g., initial weights for the first iteration; weights determined for the previous iteration for each subsequent iteration) to produce an intermediate uterine activity trace. In some embodiments, the interim uterine activity trace is in a temporary form because weights have not yet been optimized and selected. In some embodiments, contractions are identified in the intermediate uterine activity trace using a contraction detection process described below with reference to method 4200 shown in FIG. 42. In some embodiments, the following equation is then applied to extract the cost function from the signal:

number

[0173] In the above equation, E_cont is the contraction energy calculated as the transient uterine signal value spanning two-thirds of the contraction width around its peak, summed across all contractions; E_tot is the sum over the entire transient signal; A_cont is the mean contraction amplitude calculated by averaging across contractions over one-third of the contraction width around its peak; R is the range (max-min) of baseline activity amplitude between contractions; and w identifies the given weight under consideration. Note again that all sets of weights have the same length, which is equal to the number of channels. Weights representing channels that by definition should not be included in the channel set (e.g., channels derived from biopotential signals relative to the acoustic channel set) are equal to zero.

[0174] In step 4050, the best subset of weights for each weight set is selected. In some embodiments, after optimization as described above with reference to step 4040, the two weight subsets within each set compete with each other, and the best subset of weights is selected to "represent" that set. At a later stage in method 4000, the weights from the different sets will compete among themselves towards selecting the final weights to be used in the fusion process. In some embodiments, the selection process between the two competing uterine activity traces within each set (e.g., the two weight subsets within each set) is based on the following strategy:

[0175] 1. Signal-to-noise ratio of uterine activity tracing.

[0176] 2. Cost function.

[0177] 3. Contraction reliability index (defined below).

[0178] 4. Difference index.

[0179] In some embodiments, the difference index quantifies the difference between the MUA signal generated so far throughout the session (i.e., from previously analyzed recording frames) and the signal that would have resulted if these previously recorded segments had been analyzed using only the current weight set. In some embodiments, the difference index is calculated as follows:

[0180]

number

[0181] In the above, r(S prev , S current ) is the Spearman correlation coefficient between the previously determined uterine activity trace and the uterine activity trace of the current set of selected channels (excluding the data frame currently being analyzed, which cannot be compared to previous data). In some embodiments, this index varies from 0 (identical traces) to 1 (any negative correlation between the two traces). This cannot be calculated for the first analyzed data segment, and therefore this step is omitted for the first analyzed data segment.

[0182] In some embodiments, the means of the four measures are compared between the two competing subsets and against a threshold, as well as their absolute values. In some embodiments, the comparison is based on a relative difference of 10%, or, if the relative difference is not met, an absolute difference greater than zero. If one of the two subsets exhibits both a better average value of the above four measures than the other and a better relative value relative to the threshold, it is selected for the channel. If the subset with the better average falls below the threshold, a decision tree is invoked, with each measurement having a different weight in the decision process.

[0183] In some embodiments, the decision tree is implemented based on a cost function and a confidence index for the two measurements as follows: → If conf_1 > conf2 →→ If cost_2 is not valid (e.g., has an invalid numeric value such as NaN or Inf), select method 1. →→→ If cost_1 is not valid, select method 2. →→→ If both cost_1 and cost_2 are valid, calculate the relative difference between cost_1 and cost_2. If the relative difference is in favor of measure 1 by a threshold (for example, a 10% threshold) (meaning that measure 1 has a lower cost function), select measure 1; otherwise, select measure 2.

[0184] where cost_1 is the cost function of the first subset in the comparison, cost_2 is the cost function of the second subset in the comparison, conf_1 is the contraction confidence index of the first subset in the comparison, and conf_2 is the contraction confidence index of the second subset in the comparison. In some embodiments, an updated uterine activity trace is finally calculated using the selected optimized weights, and contractions are redefined based on the updated uterine activity trace.

[0185] In step 4060, the signals in each channel set are enhanced. In some embodiments, two substeps are performed to enhance the signals. In some embodiments, the first substep is to enhance data in channels with significant contraction. In some embodiments, the first substep is performed by calculating, for each channel, a measure of similarity with the weighted average signal. To this end, three metrics are examined: (1) a correlation coefficient between the weighted average signal and the individual channels that exceeds a threshold (e.g., any value greater than 0.5, such as 0.6, 0.7, 0.75, 0.8, or 0.9); (2) a first parameter ("slope") of a first-order polynomial fit between the channel data and the weighted average; and (3) an estimated error (Δ) of said fit (e.g., the difference between the first-order polynomial fit and the channel). In some embodiments, the three metrics mentioned above are checked against thresholds (e.g., 0.55 for the first metric, 0.1 for the second metric, and 0.3 for the third metric, with all thresholds configurable as needed), and the weights associated with any channels that exceed the thresholds are retained. The remaining weights (e.g., those that do not exceed the associated thresholds) are zeroed. In some embodiments, following this substep, additional iterations of gradient difference optimization are performed on the resulting weights. Next, in some embodiments, the weights of traces among the selected weights that have higher power than the weighted average are further upscaled by a factor defined to minimize the Euclidean distance between the weighted channel and the weighted average. In some embodiments, contractions and their scores are then identified on the new weighted-averaged uterine activity traces for each channel.

[0186] The second substep considers weights from previous recording segments, if any. In some embodiments, the current recording segment is assigned a contribution weight (CW) between 0 and 1, and previous weights are assigned complementary contribution weights (1-CW). In some embodiments, the contribution weight CW assigned to a given segment, segment number N, is CW=1 / N. The more previous segments there are, the lower the CW assigned to the current segment: each additional recording segment adds 1 / segment number bits of information. The weights are then adjusted to balance between the current and previous sessions according to this weighting method.

[0187] In step 4070, a final set of weights is determined. As described above, prior to step 4070, a subset of weights was selected for each of the channel sets, thereby providing uterine activity candidates from which to select a final set of weights to generate the final uterine activity output. In some embodiments, the same four metrics used in step 4050 above are used to select the final set of weights. In contrast to step 4050, in step 4070, there are more than two candidates to compare and select (e.g., four channels each have a selected set of weights as described above). Thus, in step 4070, the selection is performed iteratively. The selection in step 4070 begins by comparing metrics for a first set of weights with metrics for a second set of weights and selecting the best one. The selected one of the weight sets is compared with a third set of weights, and the best one selected from that comparison is then compared with a fourth set of weights. The best set of weights from that comparison is then selected as the final set of weights to use in generating the maternal uterine activity signal. As mentioned above, in some embodiments, the final set of weights determined in step 4070 is used as an input in step 3905 of method 3900 to generate a weighted average for the data channel.

[0188] 42, which is a flowchart of a method 4200 for identifying contractions in a uterine activity signal. In some embodiments, method 4200 is applied during the performance of method 4000, described above. In some embodiments, method 4200 is applied to the electrical uterine activity signal generated by method 1300 and / or the acoustic uterine activity signal generated by method 3000. At step 4210, method 4200 receives as input a current uterine activity channel. While method 4200 is described with reference to a single uterine activity channel, it will be apparent to those skilled in the art that method 4200 can also be performed on multiple uterine activity channels, including either sequentially and / or simultaneously.

[0189] In step 4220, smoothed and enhanced versions of the input signal received in step 4210 are calculated. In some embodiments, the smoothed version is calculated by convolving the first derivative of the signal with a Hamming window and returning the cumulative sum of the results. In some embodiments, the convolution is performed after padding the signal with its right-left flipped version at both ends (as described in more detail below with reference to step 4230), creating a continuous padded signal such that there are no edge effects at the beginning or end of the original channel data trace. In some embodiments, the padding is discarded after the smoothing process. In some embodiments, the enhanced version is calculated by calculating the hyperbolic tangent of the z-score normalized smoothed signal. In some embodiments, the enhanced version produced in this manner is a smoothed, rounded time series in which transient modulations in heartbeat peak amplitude are easily apparent and detectable. FIG. 43 is a plot showing an example input uterine activity channel 4310, a smoothed channel 4320, and an enhanced channel 4330.

[0190] In step 4230, peaks are detected. In some embodiments, contraction peaks, their prominences, and their widths are detected on the enhanced channel 4330 after right-left flipping of the signal and adding padding to both ends. In some embodiments, padding the signal with its mirrored version allows for detection of incomplete contractions at the ends of the recording, as the mirrored signal "completes" half-contractions in its mirrored version. In some embodiments, the peaks, prominences, and widths are calculated as described above with reference to step 1355 of method 1300, as shown in FIG. 13 . In some embodiments, given the smoothness of the enhanced channel 4330 and the above parameters, the detected peaks are not cluttered. In some embodiments, a second iteration of peak detection is performed, increasing the sensitivity of the peak detector, since peaks located near the edges of the signal (not considering padding) may be missed. In some embodiments, sensitivity is increased in the second iteration by reducing the thresholds applied to peak widths and distances between peaks. In some embodiments, in the first iteration, the thresholds are 30 seconds for peak width and 60 seconds for peak-to-peak distance, and in the second iteration, the thresholds are 20 seconds for peak width and 50 seconds for peak-to-peak distance. It will be apparent to those skilled in the art that different magnitudes of threshold reduction are possible. In some embodiments, if any new peaks are detected in this second iteration, they are considered only if they are located 160 seconds or less from the signal edge. In some embodiments, in the third iteration, short contractions with high prominence that may have been missed but constitute physiologically valid contractions are identified. In some embodiments, in the third iteration, sensitivity is further improved by reducing the thresholds applied to peak width and peak-to-peak distance. In some embodiments, in the third iteration, the thresholds are 20 seconds for peak width and 40 seconds for peak-to-peak distance. In some embodiments, once the final contraction peak is determined, the width of each contraction is calculated using linear interpolation of left and right points that intercept the signal at half the peak's prominence.

[0191] In step 4240, outlier peaks are identified. In some embodiments, the Euclidean distance between each pair of peak prominences is calculated, and an error estimate for each peak is calculated as the sum of the distances to other peaks. In some embodiments, outlier errors are detected. In some embodiments, outlier error values ​​are those that are more than 3 scaled median absolute deviations (MAD) from the median. In some embodiments, the scaled MAD is calculated as K*MEDIAN(ABS(A-MEDIAN(A))), where A is the value to be evaluated and K is a scaling factor. In some embodiments, the scaling factor K is equal to approximately 1.5. In some embodiments, the scaling factor K is equal to approximately 1.4826. Additionally, in some embodiments, the contraction peak height is compared to a threshold determined based on the smoothed channel 4320 generated in step 4220. In some embodiments, the threshold is determined by normalizing the smoothed channel 4320 to its maximum value and is set to 0.2 after normalization. In some embodiments, outlier peaks and peaks with maximum values ​​less than a peak height threshold are discarded as peaks.

[0192] In step 4250, an incomplete contraction is detected. In some embodiments, an incomplete contraction is one that is still ongoing when the current data segment ends. In some embodiments, a contraction is detected as incomplete if the contraction peak and contraction offset are separated by less than a minimum required time and the activity levels before and after the contraction differ by more than a certain threshold. In some embodiments, the minimum required time is 1 minute. In some embodiments, the threshold is a relative threshold that is a 30% change between the value at the determined onset and the value at the determined offset. In some embodiments, if a contraction is identified as incomplete, it is marked as such for further use. For example, in some embodiments, a contraction marked as incomplete is given less weight when calculating the overall quality of the trace (e.g., weighted by a factor of 0.5 when determining the overall SNR of the trace). In some embodiments, the incomplete contraction is completed by considering data from subsequent segments.

[0193] In step 4260, a confidence index is calculated for each contraction. In some embodiments, three confidence indexes are calculated. In some embodiments, the confidence index is completed as follows:

[0194] 1. Contraction Relative Energy. This is calculated by first calculating the energy of a contraction as the sum of data points spanning two-thirds of the contraction width around its peak, and then, for each contraction, dividing that energy by the sum of the energies of all contractions.

[0195] 2. The ratio between the mean activity in the upper third (e.g., the mean amplitude of the third of the contraction width around its peak) and the range of values ​​between contractions. The range of values ​​between contractions is calculated by taking the distribution of all data points that do not fall within a contraction and finding the difference between the 5th and 95th percentiles of this distribution.

[0196] 3. The ratio between the range of values ​​during and between contractions. The range of values ​​between contractions is calculated as in item 2 immediately above. The range of values ​​for each contraction is calculated by taking the distribution of the data points that make up the contraction data and finding the difference between the 5th and 95th percentiles of this distribution.

[0197] In step 4270, noisy contractions and small contractions are removed based on the confidence index determined in step 4260. In some embodiments, the confidence index is compared to a predefined threshold and used to eliminate noisy contractions. In some embodiments, a noisy contraction is a contraction with a confidence index below a predefined threshold. In some embodiments, the predefined threshold is compared to 0.2, where activity between the contraction and the surrounding valley is removed, and small contractions are similarly removed. In some embodiments, small contractions are contractions with a normalized peak value less than 0.2 (e.g., less than 0.2 times the maximum peak value contraction).

[0198] In step 4280, a contraction score is determined for each contraction. In some embodiments, the contraction score is used in determining the initial weights, for example, as described above with reference to step 4030 of method 4000 shown in FIG. 40. In some embodiments, two contraction scores are determined. In some embodiments, a first contraction score is calculated as the difference in mean activity levels before and after a given contraction, normalized by the contraction peak amplitude. In some embodiments, a second contraction score is calculated as the normalized prominence, calculated as the difference between the contraction peak amplitude and the mean of activity around the contraction, divided by the peak amplitude.

[0199] As discussed herein, a technical problem in the field of maternal / fetal care is that existing solutions for monitoring uterine activity (e.g., contractions) through the use of tocodynamometers and ultrasound transducers require pregnant women to wear uncomfortable sensors that may produce unreliable data when worn by obese pregnant women (e.g., the sensors may not be sensitive enough to produce usable data). As discussed further herein, exemplary embodiments present a technical solution to this technical problem through the use of various sensors (e.g., biopotential and / or acoustic sensors) integrated into comfortable wearable devices and the analysis of data obtainable by such sensors (e.g., electrodes and / or acoustic sensors) to generate signals that can be utilized to monitor uterine activity. A further technical problem in the field of maternal / fetal care is that existing solutions for analysis based on signals obtainable by sensors (e.g., biopotential and / or acoustic sensors) that can be integrated into comfortable wearable devices are limited to analyzing such signals to extract cardiac data. As discussed herein, the exemplary embodiments provide a technical solution to this technical problem through analysis of biopotential and / or acoustic data to generate signals that can monitor uterine activity (e.g., contractions).

[0200] Publications cited throughout this document are incorporated herein by reference in their entirety. While various aspects of the present invention have been illustrated above with reference to examples and embodiments, it will be understood that the scope of the present invention is defined not by the foregoing description but rather by the following claims, appropriately interpreted under the principles of patent law. Moreover, many variations will be apparent to those skilled in the art, including that the various embodiments of the inventive methodologies, inventive systems, and inventive apparatus described herein can be used in any combination with one another. Furthermore, the various steps may be performed in any desired order (and any desired steps may be added and / or any undesired steps in a particular embodiment may be eliminated).

Claims

1. 1. A computer-implemented method comprising: by at least one computer processor, (i) a plurality of electrical contraction monitoring signal channels, the plurality of electrical contraction monitoring signal channels being different from the raw biopotential input; and (ii) a plurality of acoustic uterine contraction monitoring signal channels, the plurality of acoustic uterine contraction monitoring signal channels being different from the raw acoustic input; receiving the calculating, by the at least one computer processor, a plurality of channel weights for (i) the plurality of electrical uterine contraction monitoring signal channels and (ii) the plurality of acoustic uterine contraction monitoring signal channels using a machine learning algorithm that has the electrical uterine contraction monitoring signal channels as inputs and the acoustic uterine contraction monitoring signal channels as output, each of the channel weights corresponding to either (i) a particular one of the electrical uterine contraction monitoring signal channels or (ii) a particular one of the acoustic uterine contraction monitoring signal channels; calculating, by the at least one computer processor, a weighted average of (i) a plurality of electrical uterine contraction monitoring signal channels and (ii) a plurality of acoustic uterine contraction monitoring signal channels based on the plurality of channel weights; generating, by the at least one computer processor, a combined uterine contraction monitoring signal channel based on the weighted average; A method comprising:

2. the machine learning algorithm includes a gradient descent optimization process; The computer-implemented method of claim 1 .

3. The machine learning algorithm defining, by the at least one computer processor, a plurality of channel sets, each of the plurality of channel sets including at least some of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels; defining, by the at least one computer processor, a plurality of initial weight sets, each of the plurality of initial weight sets corresponding to a particular one of the plurality of channel sets; optimizing, by the at least one computer processor, the plurality of initial weight sets to generate a plurality of optimized weight sets, each of the plurality of optimized weight sets corresponding to a particular one of the plurality of channel sets; selecting, by the at least one computer processor, a particular one of the plurality of optimized weight sets as the plurality of channel weights; determined by a process including The computer-implemented method of claim 1 .

4. optimizing the plurality of initial weight sets comprises a gradient descent process; The computer-implemented method of claim 3 .

5. The step of selecting a particular one of the plurality of optimized weight sets comprises: generating, by the at least one computer processor, a plurality of intermediate uterine activity traces, each of the plurality of intermediate uterine activity traces corresponding to a particular one of the plurality of optimized weight sets; calculating, by the at least one computer processor, for each of the plurality of optimized weight sets, (a) a signal-to-noise ratio, (b) a cost function, (c) a contraction reliability index, and (d) a difference index for a particular one of the intermediate uterine activity traces corresponding to the particular one of the plurality of optimized weight sets; calculating, by the at least one computer processor, for each particular one of the plurality of optimized weight sets, an optimized weight set average that is an average of: (a) a signal-to-noise ratio of the particular one of the plurality of optimized weight sets; (b) a cost function of the particular one of the plurality of optimized weight sets; (c) a shrinkage confidence index of the particular one of the plurality of optimized weight sets; and (d) a difference index of the particular one of the plurality of optimized weight sets; selecting, by the at least one computer processor, one of the plurality of optimized weight sets as the particular one of the plurality of optimized weight sets; The computer-implemented method of claim 3 , wherein the method is performed by a process comprising:

6. Furthermore, generating, by the at least one computer processor, a first intermediate uterine activity trace and a second intermediate uterine activity trace corresponding to a particular one of the plurality of channel sets, the first intermediate uterine activity trace corresponding to a first one of the plurality of optimized weight sets for the particular one of the plurality of channel sets, and the second intermediate uterine activity trace corresponding to a second one of the plurality of optimized weight sets for the particular one of the plurality of channel sets; calculating, by the at least one computer processor, for a first of the plurality of optimized weight sets: (a) a signal-to-noise ratio of the first intermediate uterine activity trace; (b) a cost function; (c) a contraction reliability index; and (d) a difference index; calculating, by the at least one computer processor, for a second of the plurality of optimized weight sets: (a) a signal-to-noise ratio of the second intermediate uterine activity trace; (b) a cost function; (c) a contraction reliability index; and (d) a difference index; calculating, by the at least one computer processor, for a first of the plurality of optimized weight sets, a first average that is an average of (a) the signal-to-noise ratio of the first intermediate uterine activity trace, (b) a cost function of the first of the plurality of optimized weight sets, (c) a contraction confidence index of the first of the plurality of optimized weight sets, and (d) a difference index of the first of the plurality of optimized weight sets; calculating, by the at least one computer processor, for a second of the plurality of optimized weight sets, a second average that is the average of: (a) the signal-to-noise ratio of the second intermediate uterine activity trace; (b) the cost function of the second of the plurality of optimized weight sets; (c) the contraction confidence index of the second of the plurality of optimized weight sets; and (d) the difference index of the second of the plurality of optimized weight sets; selecting, by the at least one computer processor, a first one of the plurality of optimized weight sets as the best weight set for a particular one of the plurality of channel sets based on a determination that the first average is greater than the second average; selecting, by the at least one computer processor, a second one of the plurality of optimized weight sets as the best weight set for a particular one of the plurality of channel sets based on a determination that the second average is greater than the first average; The computer-implemented method of claim 3 , comprising:

7. defining the plurality of channel sets includes defining contraction-based channel sets; The contraction-based channel set comprises: identifying, by the at least one computer processor, a set of contractions in each of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels; extracting, by the at least one computer processor, contraction features for each of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels based on the set of contractions identified for each of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels; clustering, by the at least one computer processor, the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels into a plurality of clusters; selecting, by the at least one computer processor, one of the plurality of clusters as the contraction-based channel set; The computer-implemented method of claim 3 , wherein the value is determined by a process comprising:

8. The step of defining the plurality of channel sets further comprises: and refining, by the at least one computer processor, the one of the plurality of clusters by removing channels that reduce internal agreement between the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels in the one of the plurality of clusters. The computer-implemented method of claim 7, comprising:

9. The step of defining the plurality of channel sets further comprises: adding, by the at least one computer processor, to the one of the plurality of clusters, a portion of one of the plurality of electrical uterine contraction monitoring signal channels or one of the plurality of acoustic uterine contraction monitoring signal channels that is not included in the one of the plurality of clusters. The computer-implemented method of claim 7, comprising:

10. identifying a set of contractions in each of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels; for each of the plurality of electrical contraction monitoring signal channels and the plurality of acoustic contraction monitoring signal channels; generating, by the at least one computer processor, an enhanced version of one of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels by calculating the hyperbolic tangent of a smoothed signal normalized by a Z-score; detecting, by the at least one computer processor, a candidate set of contractions in an enhanced one of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels, the candidate set of contractions including a plurality of candidate contractions; calculating, by said at least one computer processor, a plurality of confidence measures for each candidate contraction; removing, by the at least one computer processor, at least one of the candidate contractions from the set of candidate contractions, the removal being based on the confidence measure corresponding to the removed at least one of the candidate contractions, whereby removal generates the set of contractions; The computer-implemented method of claim 7 , wherein the method is performed by a process comprising:

11. The step of defining, by the at least one computer processor, the plurality of initial weight sets comprises:

4. The computer-implemented method of claim 3, comprising generating, by the at least one computer processor, a set of channel voting weights and a set of natural equal weights for each of the channel sets.

12. receiving, by the at least one computer processor, the plurality of electrical contraction monitoring signal channels includes generating at least one of the plurality of electrical contraction monitoring signal channels; At least one of the plurality of electrical contraction monitoring signal channels receiving, by the at least one computer processor, the raw biopotential inputs, each of the raw biopotential inputs being received from a corresponding one of a plurality of electrodes, each of the plurality of electrodes positioned to measure a respective one of the raw biopotential inputs of the pregnant human subject; generating, by the at least one computer processor, a plurality of signal channels from the raw biopotential input, the plurality of signal channels from the raw biopotential input including at least three signal channels; pre-processing, by the at least one computer processor, respective signal channel data for each of the signal channels to create a plurality of pre-processed signal channels, each of the pre-processed signal channels including respective pre-processed signal channel data; extracting, by the at least one computer processor, a respective plurality of R-wave peaks from the pre-processed signal channel data for each of the pre-processed signal channels to create a plurality of R-wave peak data sets, each of the R-wave peak data sets including a respective plurality of R-wave peaks; removing, by the at least one computer processor, from the plurality of R-wave peak data sets, at least one of: (a) at least one signal artifact; or (b) at least one outlier data point, wherein the at least one signal artifact is one of an electromyogram artifact or a baseline artifact; replacing, by the at least one computer processor, the at least one signal artifact, the at least one outlier data point, or both, with at least one statistical value determined based on a corresponding one of the R-wave peak data sets from which the at least one signal artifact, the at least one outlier data point, or both have been removed, to create a plurality of interpolated R-wave peak data sets; generating, by the at least one computer processor, a respective R-wave signal data set for each R-wave signal channel at a predetermined sampling rate based on the respective interpolated R-wave peak data set to create a plurality of R-wave signal channels; selecting, by the at least one computer processor, at least one first selected R-wave signal channel and at least one second selected R-wave signal channel from the plurality of R-wave signal channels based on at least one correlation between (a) a respective R-wave signal data set of at least one first specific R-wave signal channel and (b) a respective R-wave signal data set of at least one second specific R-wave signal channel; generating, by the at least one computer processor, electrical contraction monitoring data representative of an electrical uterine contraction monitoring signal based on at least the respective R-wave signal data sets of the first selected R-wave signal channel and the respective R-wave signal data sets of the second selected R-wave signal channel, thereby creating the at least one electrical uterine contraction monitoring signal channel; The computer-implemented method of claim 1 , generated by a process comprising:

13. receiving, by the at least one computer processor, the plurality of acoustic uterine contraction monitoring signal channels includes generating at least one of the plurality of acoustic uterine contraction monitoring signal channels; At least one of the plurality of acoustic uterine contraction monitoring signal channels receiving, by the at least one computer processor, the raw acoustic inputs, each of the raw acoustic inputs being received from a corresponding one of a plurality of acoustic sensors, each of the plurality of acoustic sensors positioned to measure a respective one of the raw acoustic inputs of a pregnant human subject; generating, by the at least one computer processor, a plurality of signal channels from the raw acoustic input, the plurality of signal channels from the raw acoustic input including at least three signal channels; pre-processing, by the at least one computer processor, respective signal channel data for each of the signal channels to create a plurality of pre-processed signal channels, each of the pre-processed signal channels including respective pre-processed signal channel data; extracting, by the at least one computer processor, a respective plurality of S1-S2 peaks from the pre-processed signal channel data for each of the pre-processed signal channels to create a plurality of S1-S2 peak data sets, each of the S1-S2 peak data sets including a respective plurality of S1-S2 peaks; removing, by the at least one computer processor, from the plurality of S1-S2 peak data sets, at least one of: (a) at least one signal artifact; or (b) at least one outlier data point, wherein the at least one signal artifact is one of a motion-related artifact or a baseline artifact; replacing, with the at least one computer processor, the at least one signal artifact, the at least one outlier data point, or both, with at least one statistical value determined based on a corresponding one of the S1-S2 peak data sets from which the at least one signal artifact, the at least one outlier data point, or both have been removed, to create a plurality of interpolated S1-S2 peak data sets; generating, by the at least one computer processor, a respective S1-S2 signal data set for each S1-S2 signal channel at a predetermined sampling rate based on the respective interpolated S1-S2 peak data set to create a plurality of S1-S2 signal channels; selecting, by the at least one computer processor, at least one first selected S1-S2 signal channel and at least one second selected S1-S2 signal channel from the plurality of S1-S2 signal channels based on at least one correlation between (a) a respective S1-S2 signal data set of at least one first specific S1-S2 signal channel and (b) a respective S1-S2 signal data set of at least one second specific S1-S2 signal channel; generating, by the at least one computer processor, acoustic contraction monitoring data representative of an acoustic uterine contraction monitoring signal based at least on the respective S1-S2 signal data sets of the first selected S1-S2 signal channel and the respective S1-S2 signal data sets of the second selected S1-S2 signal channel, thereby creating at least one of the plurality of acoustic uterine contraction monitoring signal channels; The computer-implemented method of claim 1 , generated by a process comprising:

Images