Fusion signal processing for maternal uterine activity detection

Abstract

A computer-implemented method includes providing, by at least one computer processor, a plurality of signal channels, wherein the plurality of signal channels includes 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, wherein each channel weight corresponds 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 of each signal channel.

Description

[0001] Cross-references to related applications

[0002] This application is an international (PCT) patent application that relates to and claims the benefit of co-owned, co-pending U.S. Provisional Patent Application No. 62 / 970,585, filed February 5, 2020, entitled “Fusion Signal Processing for Maternal Uterine Activity Detection,” the contents of which are incorporated herein by reference in their entirety. Technical Field

[0003] This invention generally relates to the monitoring of pregnant women. More specifically, this invention relates to the analysis of sensed bioelectrical and / or acoustic data to generate a computational representation of uterine activity, such as uterine contractions. Background Technology

[0004] Uterine contractions are a temporary process during which the uterine muscles shorten and the spacing between muscle cells decreases. These structural changes in the muscles lead to increased pressure within the uterine cavity, allowing the fetus to be pushed down into a lower position for delivery. During uterine contractions, the structure of the myometrial cells (i.e., uterine cells) changes, and the uterine wall thickens. Figure 1A An illustration of a relaxed uterus is shown, in which the uterine muscle walls are relaxed. Figure 1B An illustration shows a contracting uterus, where the uterine muscle walls contract and push the fetus toward the cervix.

[0005] Uterine contractions are monitored to assess the progress of labor. Typically, this is done using two sensors: a dynamometer, a strain gauge-based sensor located on the pregnant woman's abdomen; and an ultrasound transducer, also located on the abdomen. The dynamometer's signal is used to provide a labor map (“TOCO”), which is analyzed to identify uterine contractions, while the ultrasound transducer's signal is used to detect fetal heart rate, maternal heart rate, and fetal movement. However, these sensors can be uncomfortable to wear and may produce unreliable data when worn by obese pregnant women. Summary of the Invention

[0006] In some embodiments, the present application provides a specially programmed computer system comprising at least the following components: a non-transitory memory electrically storing computer executable program code; and at least one computer processor that, when executing the program code, becomes a specially programmed computing processor configured to at least: receive a plurality of biopotential signals collected at a plurality of locations on a pregnant woman's abdomen; detect R-wave peaks in the biopotential signals; extract maternal electrocardiogram ("ECG") signals from the biopotential signals; determine R-wave amplitudes in the maternal ECG signals; create an R-wave amplitude signal for each maternal ECG signal; calculate a mean of all R-wave amplitude signals; and normalize the mean to produce an electrical uterine monitoring ("EUM") signal. In some embodiments, the operations further include identifying at least one uterine contraction based on at least one corresponding peak in the EUM signal.

[0007] In some embodiments, the present application provides a method comprising: receiving a plurality of biopotential signals collected at a plurality of locations on a pregnant woman's abdomen; detecting R-wave peaks in the biopotential signals; extracting maternal electrocardiogram ("ECG") signals from the biopotential signals; determining R-wave amplitudes in the maternal ECG signals; creating an R-wave amplitude signal for each maternal ECG signal; calculating a mean of all R-wave amplitude signals; and normalizing the mean to produce an EUM signal. In some embodiments, the method further includes identifying at least one uterine contraction based on at least one corresponding peak in the EUM signal.

[0008] In one embodiment, a computer-implemented method, comprising: receiving, by at least one computer processor, a plurality of raw biopotential inputs, wherein each raw biopotential input is received from a corresponding one of a plurality of electrodes, wherein each of the plurality of electrodes is positioned so as to measure a respective one of the raw biopotential inputs of a pregnant human subject; generating, by the at least one computer processor, a plurality of signal channels from the plurality of raw biopotential inputs, wherein the plurality of signal channels comprises at least three signal channels; pre-processing, by the at least one computer processor, respective signal channel data of each signal channel to produce a plurality of pre-processed signal channels, wherein each pre-processed signal channel comprises respective pre-processed signal channel data; extracting, by the at least one computer processor, a respective plurality of R-peak from the pre-processed signal channel data of each pre-processed signal channel to produce a plurality of R-peak datasets, wherein each R-peak dataset comprises a respective plurality of R-peaks; removing, by the at least one computer processor, from the plurality of R-peak datasets at least one of: (a) at least one signal artifact and (b) at least one outlier data point, wherein the at least one signal artifact is one of an electromyography artifact and 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 respective one of the R-peak datasets from which the at least one signal artifact, the at least one outlier data point, or both, were removed; generating, by the at least one computer processor, a respective R-wave signal dataset for a respective R-wave signal channel at a predetermined sampling rate based on each respective interpolated R-peak dataset to produce 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) the respective R-wave signal dataset of the at least one first particular R-wave signal channel and (b) the respective R-wave signal dataset of the at least one second particular R-wave signal channel; generating, by the at least one computer processor, electrical uterine monitoring data representing an electrical uterine monitoring signal based on at least the respective R-wave signal dataset of the first selected R-wave signal channel and the respective R-wave signal dataset of the second selected R-wave signal channel.

[0009] In one embodiment, the computer-implemented method further includes sharpening, by the at least one computer processor, the electrical uterine monitoring data to produce a sharpened electrical uterine monitoring signal. In one embodiment, the sharpening step is omitted if the electrical uterine monitoring data is computed based on a selected one of the electrical uterine monitoring signal channels (i.e., the corrupted electrical uterine signal monitoring channel). In one embodiment, the computer-implemented method further includes post-processing the sharpened electrical monitoring signal data to produce a post-processed electrical uterine monitoring signal. In one embodiment, the sharpening step includes identifying a set of peaks in the electrical uterine monitoring signal data; determining a prominence of each peak; removing peaks having a prominence less than at least one threshold prominence value from the set of peaks; computing a mask based on the remaining peaks of the set of peaks; smoothing the mask based on a moving average window to produce a smoothed mask; and adding the smoothed mask to the electrical uterine monitoring signal data to produce the sharpened electrical uterine monitoring signal data. In one embodiment, the at least one threshold prominence value includes at least one threshold prominence value selected from the group consisting of an absolute prominence value and a relative prominence value computed based on a maximum prominence of the peaks in the set of peaks. In one embodiment, the mask includes zero values outside of a region of the remaining peaks and non-zero values inside of the region of the remaining peaks, wherein the non-zero values are computed based on a Gaussian function.

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

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

[0012] In one embodiment, the step of removing at least one of signal artifacts and outliers data points includes removing at least one electromyography artifact by a process including: 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 an inter-peak root mean square value greater than a threshold value; and replacing the corrupted peak with a median value, wherein the median value is a local median value or a global median value.

[0013] In one embodiment, the step of removing at least one of a signal artifact and an outlier data point includes removing at least one baseline artifact by a process including: identifying a change point in an R-peak in one of the plurality of R-peak data sets; subdividing the one of the plurality of R-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 of the first portion; determining a second root mean square value of the second portion; determining a balancing 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- peaks in the first portion by the balancing factor.

[0014] In one embodiment, the step of removing at least one of a signal artifact and an outlier includes removing at least one outlier according to a Grubbs test for outliers.

[0015] In one embodiment, the step of generating a respective R-wave data set based on each respective R-peak data set includes interpolating between R-peaks of each respective R-peak data set, and wherein interpolating between R-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.

[0016] In one 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 percentiles of previous intervals in which each R-wave signal channel experienced a contact issue; grouping the selected candidate R-wave signal channels into a plurality of pairs, wherein each pair includes two selected candidate R-wave channels independent of each other; calculating a correlation value for each pair; and selecting the candidate R-wave signal channels of at least one pair having a correlation value exceeding a threshold correlation value as the selected at least one first R-wave signal channel and the selected at least one second R-wave signal channel.

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

[0018] In one embodiment, the statistical value is one of a local median, a global median, and a mean.

[0019] In some embodiments, a computer-implemented method comprises: providing, by at least one computer processor, a plurality of signal channels, wherein the plurality of signal channels comprises 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, wherein each channel weight corresponds 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 of each signal channel.

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

[0021] In some embodiments, the plurality of channel weights are determined by a process comprising: defining, by the at least one computer processor, a plurality of channel sets, each of the plurality of channel sets comprising at least some of the plurality of signal channels; defining, by the at least one computer processor, a plurality of initial weight sets, wherein each of the plurality of initial weight sets corresponds 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, wherein each of the plurality of optimized weight sets corresponds 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.

[0022] In some embodiments, the step of optimizing the plurality of initial weight sets comprises a gradient descent process.

[0023] In some embodiments, the step of selecting a best one of the plurality of optimized weight sets is performed by a process comprising: generating, by the at least one computer processor, a plurality of temporary uterine activity traces, wherein each of the plurality of temporary uterine activity traces corresponds 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 of the particular one temporary uterine activity trace corresponding to the particular one of the plurality of optimized weight sets, (b) a cost function, (c) a contraction confidence measure, and (d) a difference index; calculating, by the at least one computer processor, for each particular one of the plurality of optimized weight sets, an optimized weight set mean that is a mean of (a) the signal-to-noise ratio of the particular one of the plurality of optimized weight sets, (b) the cost function of the particular one of the plurality of optimized weight sets, (c) the contraction confidence measure of the particular one of the plurality of optimized weight sets, and (d) the difference index of the particular one of the plurality of optimized weight sets; and selecting, by the at least one computer processor, a one of the plurality of optimized weight sets having a best optimized weight set mean as the best one of the plurality of optimized weight sets.

[0024] In some embodiments, the computer-implemented method further comprises: generating, by the at least one computer processor, a first provisional uterine activity trace and a second provisional uterine activity trace corresponding to a particular one of the plurality of channel sets, wherein the first provisional uterine activity trace corresponds to a first one of the plurality of sets of optimization weights for the particular one of the plurality of channel sets, and wherein the second provisional uterine activity trace corresponds to a second one of the plurality of sets of optimization weights 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 sets of optimization weights, (a) a signal-to-noise ratio of the first provisional uterine activity trace, (b) a cost function, (c) a contraction confidence measure, and (d) a difference index; calculating, by the at least one computer processor, for the second one of the plurality of sets of optimization weights, (a) a signal-to-noise ratio of the second provisional uterine activity trace, (b) a cost function, (c) a contraction confidence measure, and (d) a difference index; calculating, by the at least one computer processor, a first mean value for the first one of the plurality of sets of optimization weights, the first mean value being a mean of (a) the signal-to-noise ratio of the first provisional uterine activity trace, (b) the cost function for the first one of the plurality of sets of optimization weights, (c) the contraction confidence measure for the first one of the plurality of sets of optimization weights, and (d) the difference index for the first one of the plurality of sets of optimization weights; calculating, by the at least one computer processor, a second mean value for the second one of the plurality of sets of optimization weights, the second mean value being a mean of (a) the signal-to-noise ratio of the second provisional uterine activity trace, (b) the cost function for the second one of the plurality of sets of optimization weights, (c) the contraction confidence measure for the second one of the plurality of sets of optimization weights, and (d) the difference index for the second one of the plurality of sets of optimization weights; selecting, by the at least one computer processor, the first one of the plurality of sets of weights as the best set of weights for the particular one of the plurality of channel sets based on a determination that the first mean value is better than the second mean value; and selecting, by the at least one computer processor, the second one of the plurality of sets of weights as the best set of weights for the particular one of the plurality of channel sets based on a determination that the second mean value is better than the first mean value. In some embodiments, the computer-implemented method further comprises: enhancing, by the at least one computer processor, the plurality of channel sets prior to the step of selecting, by the at least one computer processor, the best one of the plurality of sets of optimization weights as the plurality of channel weights.

[0025] In some embodiments, the step of defining a plurality of channel sets comprises defining a contraction-based channel set, and wherein the contraction-based channel set is determined by a process comprising: identifying, by the at least one computer processor, a set of contractions in each of the plurality of signal channels; extracting, by the 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 the at least one computer processor, the plurality of signal channels into a plurality of clusters; and selecting, by the at least one computer processor, a best one of the plurality of clusters as the contraction-based channel set. In some embodiments, the step of defining a plurality of channel sets further comprises refining, by the at least one computer processor, the best one of the plurality of clusters. In some embodiments, the step of defining a plurality of channel sets further comprises adding, by the at least one computer processor, to the best one of the plurality of clusters, a portion of one signal channel that is not included in the best one of the plurality of clusters. In some embodiments, the step of identifying a set of contractions in each of the plurality of signal channels is performed by a process comprising, for each of the plurality of signal channels: generating, by the at least one computer processor, an enhanced version of the one of the plurality of signal channels; detecting, by the at least one computer processor, a set of candidate contractions in the enhanced one of the plurality of signal channels, wherein the set of candidate contractions comprises a plurality of candidate contractions; computing, by the at least one computer processor, a plurality of confidence measures for each of the candidate contractions; and removing, by the at least one computer processor, at least one of the candidate contractions from the set of candidate contractions based on the confidence measures for the candidate contractions corresponding to the at least one of the candidate contractions being eliminated, thereby producing the set of contractions.

[0026] In some embodiments, the step of defining a plurality of initial weight sets by the at least one computer processor comprises generating, by the at least one computer processor, for each channel set, a channel vote weight set and a born-equal weight set.

[0027] In some embodiments, the step of providing a plurality of signal channels includes generating at least one of the plurality of electrical uterine monitoring signal channels, and wherein the 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, wherein each raw biopotential input is received from a corresponding one of a plurality of electrodes, wherein each of the plurality of electrodes is positioned so as to measure a respective one of the raw biopotential inputs of a pregnant human subject; generating, by the at least one computer processor, a plurality of signal channels from the plurality of raw biopotential inputs, wherein the plurality of signal channels includes at least three signal channels; pre-processing, by the at least one computer processor, respective signal channel data of each signal channel to produce a plurality of pre-processed signal channels, wherein each pre-processed signal channel includes 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 of each pre-processed signal channel to produce a plurality of R-wave peak data sets, wherein each R-wave peak data set includes 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 and (b) at least one outlier data point, wherein the at least one signal artifact is one of an electromyography artifact and 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 to produce a plurality of interpolated R-wave peak data sets, the at least one statistical value being determined based on a respective 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, were removed; generating, by the at least one computer processor, a respective R-wave signal data set for a respective R-wave signal channel at a predetermined sampling rate based on each respective interpolated R-wave peak data set to produce 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) the respective R-wave signal data set of the at least one first particular R-wave signal channel and (b) the respective R-wave signal data set of the at least one second particular R-wave signal channel; and generating, by the at least one computer processor, electrical uterine monitoring data representing an electrical uterine monitoring signal based on at least 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, thereby producing the at least one electrical uterine monitoring signal channel.

[0028] In some embodiments, the step of providing a plurality of signal channels comprises generating at least one acoustic uterine monitoring signal channel of a plurality of acoustic uterine monitoring signal channels, and wherein at least one acoustic uterine monitoring signal channel of the plurality of acoustic uterine monitoring signal channels is generated by a process comprising: receiving, by the at least one computer processor, a plurality of raw acoustic inputs, wherein each raw acoustic input is received from a corresponding acoustic sensor of a plurality of acoustic sensors, wherein each of the plurality of acoustic sensors is positioned so as 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 plurality of raw acoustic inputs, wherein the plurality of signal channels comprises at least three signal channels; pre-processing, by the at least one computer processor, respective signal channel data of each signal channel to produce a plurality of pre-processed signal channels, wherein each pre-processed signal channel comprises 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 of each pre-processed signal channel to produce a plurality of S1-S2 peak data sets, wherein each S1-S2 peak data set comprises 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 and (b) at least one outlier data point, wherein the at least one signal artifact is one of a motion-related artifact and 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 respective 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, were removed, to produce 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 a respective S1-S2 signal channel based on each respective interpolated S1-S2 peak data set at a predetermined sampling rate to produce 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) the respective S1-S2 signal data set of the at least one first particular S1-S2 signal channel and (b) the respective S1-S2 signal data set of the at least one second particular S1-S2 signal channel; and generating, by the at least one computer processor, acoustic uterine monitoring data representing an acoustic uterine monitoring signal based on at least the respective S1-S2 signal data set of the first selected S1-S2 signal channel and the respective S1-S2 signal data set of the second selected S1-S2 signal channel, thereby producing the at least one acoustic uterine monitoring signal channel. BRIEF DESCRIPTION OF DRAWINGS

[0029] Figure 1A A representative uterus is shown in a non-contracted state.

[0030] Figure 1B A representative uterus is shown in a contracted state.

[0031] Figure 2 A flowchart of an example method is shown.

[0032] Figure 3 An example garment including a plurality of biopotential sensors useful for sensing data to be analyzed is shown according to Figure 2 an example method.

[0033] Figure 4A A front view of the location of a pair of ECG sensors on the abdomen of a pregnant woman is shown according to some embodiments of the present application.

[0034] Figure 4B A side view of the location of a pair of ECG sensors on the abdomen of a pregnant woman is shown according to some embodiments of the present application.

[0035] Figure 5 An example biopotential signal is shown before and after pre-processing.

[0036] Figure 6A An example biopotential signal is shown after pre-processing, with detected R-wave peaks indicated.

[0037] Figure 6B An example biopotential signal is shown after peak re-detection Figure 6A .

[0038] Figure 6C An enlarged view of a portion of the signal is shown Figure 6B .

[0039] Figure 6D An example biopotential signal is shown after checking detected peaks Figure 6B .

[0040] Figure 7A A portion of an example biopotential signal is shown including identified R-wave peaks.

[0041] Figure 7B A portion of an example biopotential signal is shown with P-waves, QRS complexes, and T-waves identified.

[0042] Figure 7C An example biopotential signal is shown including mixed maternal and fetal data.

[0043] Figure 7D A portion of the signal is shown Figure 7C , along with an initial template.

[0044] Figure 7E It shows Figure 7C A portion of the signals and the adaptation template.

[0045] Figure 7F It shows Figure 7C A portion of the signals, as well as the current template and the 0th iteration of the adaptation.

[0046] Figure 7G It shows Figure 7C A portion of the signals, as well as the current template, the 0th iteration of the adaptation, and the 1st iteration of the adaptation.

[0047] Figure 7H It shows Figure 7C A portion of the signal, as well as the current template and the parent ECG signal reconstructed based on the current template.

[0048] Figure 7I The adaptive process is illustrated as a logarithm of the error signal plotted relative to the number of iterations.

[0049] Figure 7J The extracted maternal ECG signal is shown.

[0050] Figure 8 An example of a filtered parent ECG signal is shown.

[0051] Figure 9 An exemplary parent ECG signal with R-wave peaks is shown.

[0052] Figure 10A An exemplary R-wave amplitude signal is shown.

[0053] Figure 10B An exemplary modulated R-wave amplitude signal is shown.

[0054] Figure 11A An exemplary modulated R-wave amplitude signal and the result of applying a moving average filter to it are shown.

[0055] Figure 11B An exemplary filtered R-wave amplitude signal for multiple channels within the same time window is shown.

[0056] Figure 12A It shows the basis Figure 11B The first exemplary normalized electrical uterine signal generated by the exemplary filtered R-wave amplitude signal shown is an example of this.

[0057] Figure 12B The image shows a first labor map signal, in which self-reported contractions are marked, recorded during the same time period as the exemplary normalized electrical uterine signal.

[0058] Figure 13 A flowchart of the second exemplary method is shown.

[0059] Figure 14A The second labor diagram signal, marked with self-reported contractions, is shown.

[0060] Figure 14B It shows from and Figure 14A A second exemplary electrical uterine signal was derived from bioelectrical potential data recorded during the same time period shown.

[0061] Figure 15A The diagram shows a third labor signal marked with self-reported contractions.

[0062] Figure 15B It shows from and Figure 15A A third exemplary electrical uterine signal derived from bioelectrical potential data recorded during the same time period shown.

[0063] Figure 16A The fourth labor diagram signal, marked with self-reported contractions, is shown.

[0064] Figure 16B It shows from and Figure 16A The fourth exemplary electrical uterine signal is derived from bioelectrical potential data recorded during the same time period shown.

[0065] Figure 17A The fifth labor diagram signal, marked with self-reported contractions, is shown.

[0066] Figure 17B It shows from and Figure 17A The fifth exemplary electrical uterine signal is derived from bioelectrical potential data recorded during the same time period shown.

[0067] Figure 18A An exemplary raw bioelectric potential dataset is shown.

[0068] Figure 18B It shows the basis Figure 18A An example of a filtered dataset is an example of an example of a raw dataset.

[0069] Figure 18C An exemplary raw bioelectric potential dataset is shown.

[0070] Figure 18D It shows the basis Figure 18C An example of a filtered dataset is an example of an example of a raw dataset.

[0071] Figure 18E An exemplary raw bioelectric potential dataset is shown.

[0072] Figure 18FAn exemplary filtered data set based on the exemplary raw data set of Figure 18E is shown.

[0073] Figure 18G An exemplary raw bioelectric potential data set is shown.

[0074] Figure 18H An exemplary filtered data set based on the exemplary raw data set of Figure 18G is shown.

[0075] Figure 19A An exemplary filtered data set with input peak locations is shown.

[0076] Figure 19B An exemplary filtered data set of Figure 19A with extracted peak locations is shown.

[0077] Figure 20A An exemplary filtered data set is shown.

[0078] Figure 20B An exemplary filtered data set of Figure 20A with representations of the corresponding maternal motion envelope and inter-peak absolute sum is shown.

[0079] Figure 20C An exemplary corrected data set produced by removing electromyographic artifacts from the filtered data set of Figure 20A is shown.

[0080] Figure 21A An exemplary corrected data set including baseline artifacts is shown.

[0081] Figure 21B An exemplary corrected data set of Figure 21A after removal of baseline artifacts is shown.

[0082] Figure 22A An exemplary corrected data set including outlier data points is shown.

[0083] Figure 22B An exemplary corrected data set of Figure 22A after removal of outlier data points is shown.

[0084] Figure 23A An exemplary R-wave peak signal is shown.

[0085] Figure 23B An exemplary R-wave signal generated based on the exemplary R-wave peak signal of Figure 23A is shown.

[0086] Figure 24A A set of exemplary candidate R-wave signal channels is shown.

[0087] Figure 24B A set of example selected signal channels generated based on the set of example candidate R-wave signal channels shown in Figure 24A

[0088] Figure 25A An example electrical uterine monitoring signal generated based on the set of selected signal channels shown in Figure 24B

[0089] Figure 25B An example corrected electrical uterine monitoring signal generated by applying a drift baseline removal to the example electrical uterine monitoring signal shown in Figure 25A

[0090] Figure 26 An example normalized electrical uterine monitoring signal generated based on the example corrected electrical uterine monitoring signal shown in Figure 25B

[0091] Figure 27A An example normalized electrical uterine monitoring signal shown in

[0092] Figure 27B An example sharpened mask generated based on the example normalized electrical uterine monitoring signal shown in Figure 27A

[0093] Figure 27C An example sharpened electrical uterine monitoring signal generated based on the example normalized electrical uterine monitoring signal shown in Figure 27A Figure 27B

[0094] Figure 28 An example post-processed electrical uterine monitoring signal shown in

[0095] Figure 29 A labor chart signal corresponding to the example post-processed electrical uterine monitoring signal shown in Figure 28

[0096] A flowchart of a third example method shown in Figure 30

[0097] An example pre-processed data set shown in Figure 31A

[0098] An enlarged view of a portion of the example pre-processed data set shown in Figure 31B Figure 31A An example R-wave peak extracted in the example pre-processed data set shown in

[0099] Figure 32A

[0100] Figure 32B ​​​​​​​​​An enlarged view of an R-wave peak extracted in an exemplary pre-processed data set is shown.

[0101] Figure 32C An exemplary R-wave amplitude signal is shown.

[0102] Figure 32D An exemplary R-wave amplitude signal is shown over a larger time window.

[0103] Figure 33 A filtered exemplary R-wave amplitude signal is shown.

[0104] Figure 34 Four data channels of exemplary R-wave data are shown.

[0105] Figure 35A A sixth tocogram signal is shown.

[0106] Figure 35B A first exemplary acoustic uterine signal derived from acoustic data recorded during the same time period as shown in Figure 35A A second exemplary acoustic uterine signal derived from acoustic data recorded during the same time period as shown in

[0107] Figure 36A A seventh tocogram signal is shown.

[0108] Figure 36B A second exemplary acoustic uterine signal derived from acoustic data recorded during the same time period as shown in Figure 36A A third exemplary acoustic uterine signal derived from acoustic data recorded during the same time period as shown in

[0109] Figure 37A An eighth tocogram signal is shown.

[0110] Figure 37B A third exemplary acoustic uterine signal derived from acoustic data recorded during the same time period as shown in Figure 37A A third exemplary acoustic uterine signal derived from acoustic data recorded during the same time period as shown in

[0111] Figure 38 A set of ECG-based EUM processing signals and PCG-based processing signals is shown.

[0112] Figure 39 A fusion process to generate a uterine activity signal from the ECG-based EUM processing signals and PCG-based processing signals is shown.

[0113] Figure 40 A process to generate a set of weights for generating a uterine activity signal based on electrical uterine monitoring data and acoustic uterine monitoring data is shown.

[0114] Figure 41 A process to define an initial set of channels in a method for Figure 40

[0115] ​ Figure 42 The process of shrinking within the identifier dataset is shown.

[0116] Figure 43 As shown by Figure 42 The generated exemplary uterine monitoring signal and exemplary smoothed and enhanced signal.

[0117] Figure 44 It shows the result of Figure 39 The method generates exemplary electrical uterine monitoring data signals, exemplary acoustic uterine monitoring data signals, and exemplary output uterine monitoring signals. Detailed Implementation

[0118] Other objects and advantages of the invention will become apparent from the following description taken in conjunction with the accompanying drawings, amidst the disclosed benefits and improvements. Detailed embodiments of the invention are disclosed herein; however, it should be understood that the disclosed embodiments are merely illustrative and the invention can be practiced in various forms. Furthermore, each example given in connection with the various embodiments of the invention is illustrative and not restrictive.

[0119] Throughout the specification and claims, the following terms take their explicitly associated meanings unless the context clearly specifies otherwise. The phrases “in one embodiment,” “in an embodiment,” and “in some embodiments” as used herein do not necessarily refer to the same embodiment, although they may. Furthermore, the phrases “in another embodiment” and “in some other embodiments” as used herein do not necessarily refer to different embodiments, although they may refer to different embodiments. Therefore, as described below, various embodiments of the invention can be readily combined without departing from the scope or spirit of the invention.

[0120] As used herein, the term "based on" is not exclusive and allows for consideration based on additional factors not described unless the context explicitly states otherwise. Furthermore, throughout the specification, the meanings of "a," "an," and "the" include plural references. The meaning of "in" includes both "in" and "on." The scope discussed herein is inclusive (e.g., the scope of "between 0 and 2" includes the values ​​0 and 2, and all values ​​in between).

[0121] As used herein, the term “contact area” includes the contact area between skin-to-skin contacts in a pregnant human subject, i.e., the surface area through which an electric current can pass between skin-to-skin contacts in a pregnant human subject.

[0122] In some embodiments, the present application provides a method for extracting labor-like signals from bio-potential data, i.e., data describing electrical potentials recorded at points on the skin of a human body using skin contacts, commonly referred to as electrodes. In some embodiments, the present application provides a method for detecting uterine contractions from bio-potential data. In some embodiments, the bio-potential data is obtained using non-contact electrodes positioned on or near desired points on the human body.

[0123] In some embodiments, the present application 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.

[0124] 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, for example, a plurality of electrodes configured to detect fetal electrocardiogram signals are included into an article, such as a belt, a patch, and the like, and the article is worn or placed on the body of the pregnant human subject. Figure 3 An exemplary garment 300 is shown, which includes eight electrodes 310 included into the garment 300 so that when the garment 300 is worn by a subject, the electrodes are positioned around the abdomen of a pregnant human subject. In some embodiments, the garment 300 includes four acoustic sensors 320 included into the garment 300 so that when the garment 300 is worn by a subject, the acoustic sensors are positioned around the abdomen of a pregnant human subject. In some embodiments, each acoustic sensor 320 is one of the acoustic sensors described in U.S. Patent No. 9,713,430. Figure 4A A front view of the positions of eight electrodes 310 on the abdomen of a pregnant woman is shown, in accordance with some embodiments of the present application. Figure 4B A side view of the positions of eight electrodes 310 on the abdomen of a pregnant woman is shown, in accordance with some embodiments of the present application.

[0125] Figure 2A flowchart showing an example inventive method 200 is shown. In some embodiments, an example inventive computing device programmed / configured according to method 200 is operable to receive as input raw biopotential data measured by a plurality of electrodes positioned on the skin of a pregnant human subject, and analyze such input to produce a laborogram-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. In some embodiments, the number of electrodes is between 2 and 40. 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 number of electrodes is between 8 and 10. In some embodiments, the number of electrodes is between 8 and 20. In some embodiments, the number of electrodes is between 8 and 30. In some embodiments, the number of electrodes is between 8 and 40. In some embodiments, the number of electrodes is 8. In some embodiments, an example inventive computing device programmed / configured according to method 200 is operable to receive as input a maternal ECG signal that has been extracted from the raw biopotential data (e.g., by separating from a fetal ECG signal that forms part of the same raw biopotential data). In some embodiments, an example inventive computing device is programmed / configured according to method 200 via instructions stored in a non-transitory computer readable medium. In some embodiments, an example inventive computing device includes at least one computer processor that, when executing the instructions, becomes a specially programmed computer processor programmed / configured according to method 200.

[0126] In some embodiments, the exemplary inventive computing device is programmed / configured to continuously perform one or more steps of the method 200 along a moving time window. In some embodiments, the moving time window has a predetermined length. In some embodiments, the predetermined length is 60 seconds. In some embodiments, the exemplary inventive computing device is programmed / configured to continuously perform one or more steps of the 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.

[0127] In step 210, the exemplary inventive computing device is programmed / configured to receive raw biopotential data as input and to pre-process it. In some embodiments, the raw biopotential data is recorded by using at least two electrodes located in the vicinity of 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 located at a point away from the uterus of the subject. In some embodiments, a biopotential signal is recorded at each of several points around the abdomen of the pregnant subject. In some embodiments, a biopotential signal is recorded at each of eight points around the abdomen of the pregnant subject. In some embodiments, the biopotential data is recorded at a rate of 1000 samples per second. In some embodiments, the biopotential data is up-sampled at a rate of 1000 samples per second. In some embodiments, the biopotential data is recorded at a sampling rate between 100 and 10000 samples per second. In some embodiments, the biopotential data is up-sampled at a sampling rate between 100 and 10000 samples per second. In some embodiments, the pre-processing includes baseline removal (e.g., using a median filter and / or a moving average filter). In some embodiments, the pre-processing includes low-pass filtering. In some embodiments, the pre-processing includes low-pass filtering at 85 Hz. In some embodiments, the pre-processing includes power-line interference removal. Figure 5 A portion of the raw biopotential data signal is shown before and after pre-processing.

[0128] In step 220, the example inventive computing device is programmed / configured to detect maternal R-peak in the pre-processed biopotential data resulting from performing step 210. In some embodiments, R-peak is detected on a 10 second segment of each data signal. In some embodiments, detection of R-peak begins by analyzing derivatives, thresholds, and distances. In some embodiments, detecting R-peak in each data signal includes computing a first derivative of the data signal in the 10 second segment, identifying R-peak in the 10 second segment by identifying zero-crossings of the first derivative, and excluding identified peaks having any of: (a) an absolute value less than a predetermined R-peak threshold absolute value, and (b) a distance between adjacent identified R-peak less than a predetermined R-peak threshold distance. In some embodiments, detection of R-peak is performed in a manner similar to electrocardiogram peak detection described in U.S. Patent No. 9,392,952, the contents of which are incorporated by reference herein in their entirety. Figure 6A A pre-processed biopotential data signal is shown, in which R-peaks detected as described above are indicated with asterisks.

[0129] In some embodiments, detection of R-peak of step 220 continues a peak re-detection process. In some embodiments, the peak re-detection process includes automatic gain control (“AGC”) analysis to detect windows in which the number of peaks is significantly different. In some embodiments, the peak re-detection process includes cross-correlation analysis. In some embodiments, the peak re-detection process includes AGC analysis and cross-correlation analysis. In some embodiments, the AGC analysis is adapted to overcome false negatives. In some embodiments, the cross-correlation analysis is adapted to remove false positives. Figure 6B A data signal after peak re-detection is shown, in which R-peaks re-detected as described above are indicated with asterisks. Figure 6C An enlarged view of a portion of the data signal of Figure 6B is shown.

[0130] In some embodiments, the R-wave peak detection of step 220 continues to build a global peak array. In some embodiments, the global peak array is created from multiple data channels (e.g., each channel corresponds to one or more electrodes 310). In some embodiments, the signal of each channel is given a quality score based on the relative energy of the peaks. In some embodiments, the relative energy of a peak refers to the energy of the peak relative to the total energy of the signal being processed. In some embodiments, the energy of a peak is calculated by computing the root mean square (“RMS”) of the QRS complex containing the R-wave peak, and the energy of the signal is calculated by computing the RMS of the signal. In some embodiments, the relative energy of a peak is calculated by computing 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, with signals from other channels also considered based on a voting mechanism. In some embodiments, after the global peak array is constructed based on the best lead, each remaining channel votes on each peak. If a given peak is contained in the global peak array constructed based on the best lead (e.g., as detected in the peak detection described above), the channel votes for the given peak (e.g., gives a vote value of “1”), and if the peak is not contained, the channel votes against the peak (e.g., gives a vote value of “0”). Peaks that receive more votes are considered to be higher quality peaks. In some embodiments, a peak is retained in the global peak array if the number of votes for the peak is greater than a threshold value. In some embodiments, the threshold number of votes is half the total number of channels. In some embodiments, if the number of votes for a peak is less than the threshold value, an additional test is performed on the peak. In some embodiments, the additional test includes computing the correlation of the peak in the best lead channel to a template computed 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, a further correlation is computed for all leads that voted for the peak (i.e., not just the best lead 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, and if the further correlation is less than the second threshold correlation value, the peak is excluded from the global peak array. In some embodiments, the second threshold correlation value is 0.85.

[0131] In some embodiments, once created, the global peak array is checked using physiologic measurements. In some embodiments, the checking is performed by the exemplary inventive computing device described in U.S. Patent No. 9,392,952, the contents of which are incorporated by reference herein in their entirety. In some embodiments, the physiologic parameters include R-R intervals, mean and standard deviation; and heart rate and heart rate variability. In some embodiments, the checking includes cross-correlation to overcome false negatives. Figure 6D Data signals are shown after the global peak array is created and checked as described above. In Figure 6D the peak denoted by a circled asterisk represents a previously detected R-wave peak (e.g., as shown in Figure 6A the peak denoted by a circled asterisk represents a previously detected R-wave peak (e.g., as shown in

[0132] In some embodiments, if the initial step of R-wave detection is unsuccessful (i.e., if no R-wave peak is detected on a given sample), an independent component analysis (“ICA”) algorithm is applied to the data samples, and the preceding portion of step 220 is repeated. In some embodiments, the exemplary ICA algorithm is, for example, but not limited to, the FAST ICA algorithm. In some embodiments, the FAST ICA algorithm is utilized, for example, in accordance with Hyvarinen et al., “Independent component analysis: Algorithms and applications,” Neural Networks 13(4-5): 411-430 (2000).

[0133] With continued reference to Figure 2 In step 230, the exemplary inventive computing device is programmed / configured to extract a maternal ECG signal from the signal comprising maternal and fetal data. In some embodiments, when the exemplary inventive computing device programmed / configured to perform method 200 receives the maternal ECG signal as input after extraction from the mixed maternal-fetal data, the exemplary inventive computing device is programmed / configured to skip step 230. Figure 7AA portion of the signal is shown in which the R-wave peaks have been identified and includes both maternal and fetal signals. Without intending to be bound by any particular theory, a major challenge involved in extracting the maternal ECG signal is that each maternal heartbeat is different from all other maternal heartbeats. In some embodiments, this challenge is addressed by identifying each maternal heartbeat using an adaptive reconstruction scheme. In some embodiments, the extraction process begins with segmenting the ECG signal into signals of three origins. In some embodiments, this segmentation includes using a curve length transform to find 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 An exemplary ECG signal including these segments is shown.

[0134] Following the curve length transform, step 230 continues by extracting the maternal signal using an adaptive template. In some embodiments, the template is adapted for isolating the current beat. In some embodiments, extracting the maternal signal using an adaptive template is performed as described in U.S. Patent No. 9,392,952, the contents of which are incorporated by reference herein in their entirety. In some embodiments, the process includes starting with a current template and using an iterative process to adjust the current template to the current beat. In some embodiments, for each segment of the signal (i.e., P-wave, QRS complex, and T-wave), a multiplier is defined (referred to as P_mult, QRS_mult, and T_mult, respectively). In some embodiments, a shift parameter is also defined. In some embodiments, the extraction uses a Levenberg-Marquardt non-linear least mean square algorithm as follows:

[0135]

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

[0137] E = || φ m - φ c || 2

[0138] In the above expression, φ m represents the current beat ECG, φ crepresents a reconstructed ECG. In some embodiments, the method provides a local, stable, and repeatable solution. In some embodiments, the iteration continues until the relative residual 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 approximately -20 db. In some embodiments, the threshold is approximately -20 db.

[0139] Figure 7C An exemplary signal including mixed maternal and fetal data is shown. Figure 7D A portion of the signal of Figure 7C is shown along with an initial template for comparison. Figure 7E A portion of the signal of Figure 7C is shown along with an adapted template for comparison. Figure 7F A portion of the signal of Figure 7C is shown along with a current template and a 0th iteration of adaptation. Figure 7G A portion of the signal of Figure 7C is shown along with a current template, a 0th iteration of adaptation, and a 1st iteration of adaptation. Figure 7H A portion of the signal of Figure 7C is shown along with a current template and a reconstructed ECG signal (e.g., a maternal ECG signal) based on the current template. Figure 7I A log of error signals versus iterations of adaptation is shown. Figure 7J An extracted maternal ECG signal is shown.

[0140] With continued reference to Figure 2 , in step 240, the exemplary inventive computing device is programmed / configured to perform signal cleaning on the maternal signal extracted in step 230. In some embodiments, the cleaning of step 240 includes filtering. In some embodiments, the filtering includes removing a baseline 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. Figure 8 A portion of an exemplary filtered maternal ECG is shown after performance of step 240.

[0141] With continued reference to Figure 2 In step 250, the exemplary inventive computing device is programmed / configured to compute R-wave amplitudes of the filtered maternal ECG signal resulting from performance of step 240. In some embodiments, the R-wave amplitudes are computed 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 computing the amplitudes of various R-waves. In some embodiments, the amplitudes are computed as the values (e.g., signal amplitudes) of the maternal ECG signal at each detected peak location. Figure 9 An exemplary extracted maternal ECG signal is shown, in which R-wave peaks are annotated with circles.

[0142] With continued reference to Figure 2 In step 260, the exemplary inventive computing device is programmed / configured to create an R-wave amplitude signal over time based on the R-wave amplitudes computed in step 250. In some embodiments, the computed R-wave peaks are not uniformly sampled over time. Thus, in some embodiments, step 260 is performed in order to resample the R-wave amplitudes in such a way that these amplitudes will be uniformly sampled 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 computed in step 250 and resampling the connected R-wave amplitude values. In some embodiments, the resampling includes interpolation with defined query time points. In some embodiments, the interpolation includes linear interpolation. In some embodiments, the interpolation includes spline interpolation. In some embodiments, the interpolation includes cubic interpolation. In some embodiments, the query points define the time points at which interpolation should occur. Figure 10A An exemplary R-wave amplitude signal created in step 260 based on the R-wave amplitudes from step 250 is shown. In Figure 10A In step 260, the exemplary inventive computing device is programmed / configured to create an R-wave amplitude signal over time based on the R-wave amplitudes computed in step 250. In some embodiments, the computed R-wave peaks are not uniformly sampled over time. Thus, in some embodiments, step 260 is performed in order to resample the R-wave amplitudes in such a way that these amplitudes will be uniformly sampled 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 computed in step 250 and resampling the connected R-wave amplitude values. In some embodiments, the resampling includes interpolation with defined query time points. In some embodiments, the interpolation includes linear interpolation. In some embodiments, the interpolation includes spline interpolation. In some embodiments, the interpolation includes cubic interpolation. In some embodiments, the query points define the time points at which interpolation should occur. Figure 8 In step 260, the exemplary inventive computing device is programmed / configured to create an R-wave amplitude signal over time based on the R-wave amplitudes computed in step 250. In some embodiments, the computed R-wave peaks are not uniformly sampled over time. Thus, in some embodiments, step 260 is performed in order to resample the R-wave amplitudes in such a way that these amplitudes will be uniformly sampled 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 computed in step 250 and resampling the connected R-wave amplitude values. In some embodiments, the resampling includes interpolation with defined query time points. In some embodiments, the interpolation includes linear interpolation. In some embodiments, the interpolation includes spline interpolation. In some embodiments, the interpolation includes cubic interpolation. In some embodiments, the query points define the time points at which interpolation should occur. Figure 10B Modulation of the R-wave amplitude signal over a larger time window is shown.

[0143] With continued reference to Figure 2At step 270, the example inventive computing device is programmed / configured to clean the R-wave amplitude signal by applying a moving average filter. In some embodiments, the moving average filter is applied to clean the R-wave amplitude signal of high frequency variations. In some embodiments, the moving average filter is applied over a predetermined time window. In some embodiments, the length of the time window is between 1 second and 10 minutes. In some embodiments, the length of the time window is between 1 second and 1 minute. In some embodiments, the length of the time window is between 1 second and 30 seconds. In some embodiments, the length of the time window is 20 seconds. Figure 11A The R-wave amplitude signal of Figure 10B is shown, where the signal resulting from the application of the moving average filter is shown in bold along the middle of the R-wave amplitude signal. As described above, in some embodiments, multiple data channels are considered inputs to the method 200. Figure 11B A plot of the filtered R-wave amplitude signals of multiple channels over the same time window is shown.

[0144] Continuing with Figure 2 , at step 280, the example inventive computing device is programmed / configured to compute an average signal of all filtered R-wave signals (e.g., as shown in Figure 11B ) per unit time. In some embodiments, a single average signal is computed at each time point at which samples exist. 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 that is uniformly sampled over time. At step 290, the example inventive computing device is programmed / configured to normalize the signal computed at step 280. In some embodiments, the signal is normalized by dividing by a constant factor. In some embodiments, the constant factor is between 2 volts and 1000 volts. In some embodiments, the constant factor is 50 volts. Figure 12A An example normalized electrical uterine signal after steps 280 and 290 are performed is shown. Figure 12B A labor graph signal generated over an example inventive computing system time period is shown, where contractions self-reported by the mother are indicated by vertical lines. Reference is made to Figure 12A and Figure 12B It can be seen that Figure 12A the peaks of the example normalized electrical uterine signal in Figure 12Bself-reported contractions shown in FIG. 2. Thus, in some embodiments, a normalized electrical uterine monitoring (“EUM”) signal (e.g., the signal shown in FIG. 2) generated by performing the exemplary method 200 is suitable for use in identifying contractions. In some embodiments, a contraction is identified by identifying a peak in the EUM signal. Figure 12A

[0145] In some embodiments, the present disclosure relates to a specially programmed computer system comprising at least the following components: a non-transitory memory electrically storing computer executable program code; and at least one computer processor that, when executing the program code, becomes a specially programmed computing processor configured to at least: receive a plurality of biopotential signals collected at a plurality of locations on the abdomen of a pregnant woman; detect R-wave peaks in the biopotential signals; extract a maternal electrocardiogram (“ECG”) signal from the biopotential signals; determine R-wave amplitudes in the maternal ECG signal; create an R-wave amplitude signal for each maternal ECG signal; calculate a mean of all R-wave amplitude signals; and normalize the mean to produce an electrical uterine monitoring (“EUM”) signal. In some embodiments, the operations further comprise identifying at least one uterine contraction based on at least one corresponding peak in the EUM signal.

[0146] Figure 13 ​A flowchart showing an example inventive method 1300 is shown. In some embodiments, an example inventive computing device programmed / configured according to method 1300 is operable to receive as input raw biopotential data measured by a plurality of electrodes positioned on the skin of a pregnant human subject, and analyze such input to produce a laborogram-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. In some embodiments, the number of electrodes is between 2 and 40. 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 number of electrodes is between 8 and 10. In some embodiments, the number of electrodes is between 8 and 20. In some embodiments, the number of electrodes is between 8 and 30. In some embodiments, the number of electrodes is between 8 and 40. In some embodiments, the number of electrodes is 8. In some embodiments, an example inventive computing device programmed / configured according to method 1300 is operable to receive as input a maternal ECG signal that has been extracted from the raw biopotential data (e.g., by separating from a fetal ECG signal that forms part of the same raw biopotential data). In some embodiments, an example inventive computing device is programmed / configured according to method 1300 via instructions stored in a non-transitory computer readable medium. In some embodiments, an example inventive computing device includes at least one computer processor that, when executing the instructions, becomes a specially programmed computer processor programmed / configured according to method 1300.

[0147] In some embodiments, the exemplary inventive computing device is programmed / configured to continuously perform one or more steps of the method 1300 along a moving time window. In some embodiments, the moving time window has a predetermined length. In some embodiments, the predetermined length is 60 seconds. In some embodiments, the exemplary inventive computing device is programmed / configured to continuously perform one or more steps of the 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.

[0148] In step 1305, the exemplary inventive computing device is programmed / configured to receive raw biopotential data as input. Exemplary raw biopotential data is shown in FIGS. 1301-1304. Figure 18A , Figure 18C , Figure 18E and Figure 18G In some embodiments, the raw biopotential data is recorded using at least two electrodes positioned proximate to the skin of a 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 positioned at a point distal from the subject’s uterus. In some embodiments, a biopotential signal is recorded at each of several points around the pregnant subject’s abdomen. In some embodiments, a biopotential signal is recorded at each of eight points around the pregnant subject’s abdomen. In some embodiments, the biopotential data is recorded at a rate of 1000 samples per second. In some embodiments, the biopotential data is upsampled at a rate of 1000 samples per second. In some embodiments, the biopotential data is recorded at a sampling rate between 100 and 10000 samples per second. In some embodiments, the biopotential data is upsampled at a sampling rate between 100 and 10000 samples per second. In some embodiments, the steps of the method 1300 between the reception of raw data and the channel selection (i.e., steps 1310-1335) are performed for each of a plurality of signal channels, wherein each signal channel is generated by the exemplary inventive computing device as a difference between biopotential signals recorded by a particular pair of electrodes. In some embodiments, wherein the method 1300 is performed using data recorded at electrodes positioned as shown in FIGS. 1301-1304, the channels are identified as follows: Figure 4A and Figure 4B

[0149] • Channel 1: Al-A4

[0150] ​• Channel 2: A2-A3

[0151] • Channel 3: A2-A4

[0152] • Channel 4: A4-A3

[0153] • Channel 5: B1-B3

[0154] • Channel 6: B1-B2

[0155] • Channel 7: B3-B2

[0156] • Channel 8: A1-A3

[0157] In step 1310, the example inventive computing device is programmed / configured to pre-process the signal channels determined based on the raw biopotential data to produce a plurality of pre-processed signal channels. In some embodiments, the pre-processing includes one or more filters. In some embodiments, the pre-processing includes more than one filter. In some embodiments, the pre-processing includes a DC removal filter, a power line filter, and a high-pass filter. In some embodiments, the DC removal filter removes the mean of the raw data for the current processing interval. In some embodiments, the power line filter includes a 10thorder band-stop infinite impulse response ("IIR") filter configured to minimize any noise at a preconfigured frequency in the data. In some embodiments, the preconfigured frequency is 50 Hz, and the power line filter includes cutoff frequencies of 49.5 Hz and 50.5 Hz. In some embodiments, the preconfigured 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 drift baseline from the signal, where the baseline is computed by a moving average window of a predetermined length. 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 milliseconds and 250 milliseconds. In some embodiments, the predetermined length is between 175 milliseconds and 225 milliseconds. In some embodiments, the predetermined length is approximately 200 milliseconds. In some embodiments, the predetermined length is 201 milliseconds (i.e., 50 samples at a sampling rate of 250 samples per second) long. In some embodiments, the baseline includes data from frequencies below 5 Hz, so the signal is high-pass filtered at approximately 5 Hz. The pre-processed data generated based on the raw biopotential data shown in Figure 18A , Figure 18C , Figure 18E and Figure 18G are shown in Figure 18B , Figure 18D , Figure 18F and Figure 18H , respectively.

[0158] With continued reference to step 1310, in some embodiments, after applying the above-described filters, each data channel is checked for contact problems. In some embodiments, a contact problem in each data channel is identified 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 in the peak relative energy of the data channel. In some embodiments, a data channel is identified as corrupted if the RMS value of the data channel is 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 and three local voltage units. Figure 18A and Figure 18B An exemplary data channel that is identified as corrupted on this basis is shown. In some embodiments, a data channel is identified as corrupted if the SNR value of the data channel is 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. Figure 18C and Figure 18D An exemplary data channel that is identified as corrupted on this basis is shown. In some embodiments, a data channel is identified as corrupted if the relative R-wave peak energy of the data channel varies from one interval to another by more than a threshold variation amount. In some embodiments, the threshold variation amount is 250%. In some embodiments, the threshold variation amount is between 200% and 300%. In some embodiments, the threshold variation amount is between 150% and 350%. Figure 18E and Figure 18F An exemplary data channel that is identified as corrupted on this basis is shown. In some embodiments, a data channel is identified as corrupted if the relative R-wave peak energy of the data channel varies from one interval to another by more than a threshold variation amount. In some embodiments, the threshold variation amount is 250%. In some embodiments, the threshold variation amount is between 200% and 300%. In some embodiments, the threshold variation amount is between 150% and 350%. Figure 18G and Figure 18H An exemplary data channel that is not identified as corrupted for any of the above reasons is shown in FIG. 7.

[0159] In step 1315, the example inventive computing device is programmed / configured to extract R-wave peaks from the pre-processed signal channels to produce an R-wave peak dataset. In some embodiments, step 1315 uses known maternal ECG peaks as input. In some embodiments, step 1315 uses known 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 maternal ECG peak locations using pre-processed data (e.g., as produced by step 1310) and known maternal ECG peaks. In some embodiments, peak location refinement 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 at the maximum point of the R-wave of each filtered signal. In some embodiments, the window includes a positive or negative 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 milliseconds and 250 milliseconds. In some embodiments, the predetermined length is between 175 milliseconds and 225 milliseconds. In some embodiments, the predetermined length is approximately 200 milliseconds. In some embodiments, the window includes a positive or negative number of samples in a range between 1 sample and 100 samples. Figure 19A and Figure 19B FIGS. 13A and 13B respectively show plots of a known maternal ECG peak and an extracted R-wave peak in an example R-wave peak dataset.

[0160] In step 1320, the example inventive computing device is programmed / configured to remove electromyography (“EMG”) artifacts from data including the pre-processed data produced by step 1310 and the R-wave peaks extracted in step 1315. Figure 20A FIG. 14 shows example pre-processed data used as input for step 1320. In some embodiments, removal of EMG artifacts is performed in order to correct peaks with high amplitude where high frequency energy is increased, which is often but not always derived from maternal EMG activity. Other sources of such energy are high power line noise and high fetal activity. In some embodiments, removal of EMG artifacts includes finding corrupted peaks and replacing these peaks with a median. In some embodiments, finding corrupted peaks includes calculating an inter-peak RMS value based on the following equation:

[0161] A first step in correcting the artifacts is finding the corrupted peaks. To do so, an inter-peak root mean square value is calculated, so:

[0162] Inter-peak RMS(iPeak) = RMS(Signal(peak position(iPeak)+1 : peak position(iPeak+1)-1))

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

[0164] In some embodiments, finding corrupted peaks further includes finding abnormal peaks in maternal physical activity ("MPA") data sets. In some embodiments, such a signal (hereinafter referred to as an "envelope signal") is extracted as follows:

[0165] In some embodiments, physical activity data is collected using a motion sensor. In some embodiments, the motion sensor includes a three-axis accelerometer and a three-axis gyroscope. In some embodiments, the motion sensor samples at 50 samples per second (50sps). In some embodiments, the sensor is located on the same sensing device (e.g., a wearable device) that contains electrodes for collecting biopotential data for performing method 1300 as a whole (e.g., a garment 300 as shown). Figure 3

[0166] In some embodiments, raw motion data is converted. In some embodiments, in the case of accelerometer raw data, the raw motion data is converted to g units, and in the case of gyroscope raw data, to degrees per second. In some embodiments, the converted data is examined to distinguish between valid and invalid signals by determining whether the raw signal is saturated (e.g., the raw signal has a constant maximum possible value). In some embodiments, a signal envelope is extracted as follows. First, in some embodiments, the position change of the data is examined. Because position change is 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 filter order of the high-pass filter is 400, and the frequency is 1 hertz. In some embodiments, to eliminate any non-physiological motion, a low-pass FIR filter is also applied. In some embodiments, the filter order of the low-pass filter is 400, and the frequency is 12 hertz. Also applied is a (400 order, fc=12 Hz [1]). In some embodiments, after filtering, the magnitude of the accelerometer vector is calculated according to the following equation:

[0167]

[0168] ​In this equation, 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 for the gyroscope data is calculated according to the following equation:

[0169]

[0170] In this equation, GyroManitudeVector(iSample) represents the square root of the sum of the squares of the three gyroscope axes (e.g., x, y, and z) for sample number iSample. In some embodiments, after the accelerometer magnitude vector and the gyroscope magnitude vector are calculated, the envelopes of the gyroscope magnitude vector and the accelerometer magnitude vector are extracted by applying an RMS window to the gyroscope magnitude vector and the accelerometer magnitude vector, respectively. In some embodiments, the length of the RMS window is 50 samples. In some embodiments, after the envelopes of the gyroscope magnitude vector and the accelerometer magnitude vector are extracted, the two envelopes are averaged (e.g., mean, median, etc.) to produce the MPA motion envelope.

[0171] In some embodiments, the peaks in the MPA motion envelope are defined according to the following steps:

[0172] MotionEnvelopePeak = find (MotionEnvelope > P 95％ (MotionEnvelope))

[0173] MotionEnvelopePeakStart = MotionEnvelopePeak - 2*PeakWidth

[0174] MotionEnvelopePeakOffset = MotionEnvelopePeak + 2*PeakWidth

[0175] In the above, the PeakWidth is defined as the distance between the peak and the first point where the envelope reaches 50% of the peak value, and P 95％ (x) is the 95th percentile of x. Figure 20B Data signals of Figure 20A and the corresponding motion envelope and inter-peak absolute sum (i.e., the sum of the absolute values of all samples falling between adjacent peaks) calculated according to the above are shown.

[0176] In some embodiments, a peak is determined to be corrupted if the peak is as follows:

[0177] 1) a peak with an inter-peak RMS higher than 20 local voltage units

[0178] 2) a peak with an inter-peak RMS higher than 8 local voltage units if the signal check phase concludes that there is a contact problem in the current processing interval

[0179] 3) If the signal check phase concludes that there is a contact problem in the current processing interval, but more than 50% of the points have a peak-to-peak RMS above 8 local voltage units, use a 20 local voltage unit threshold

[0180] 4) Peaks located around the start and offset of the motion envelope are suspected to be damaged. The peak-to-peak RMS of these points should exceed 6 for the peak to be concluded as damaged.

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

[0182] Local Median = Median (Peak Signal (Damaged Peak - 10 : Damaged Peak + 10))

[0183] In some embodiments, the damaged data points themselves are excluded from the above calculation and replaced by a statistical value (e.g., global median, local median, mean, etc.). In some embodiments, if 7 or fewer values are to be used after exclusion, the global median is used as the local median, where the global median is calculated using standard techniques:

[0184] Local Median = Global Median = Median (Signal)

[0185] In some embodiments, if the absolute difference between the local median and the global median exceeds 0.1, the local median is used instead of the amplitude of the damaged data point, otherwise the global median is used instead of the amplitude of the damaged peak. Figure 20C Exemplary data sets are shown Figure 20A and Figure 20B where damaged peaks as described above are replaced.

[0186] With continued reference to Figure 13 , in step 1325, the exemplary inventive computing device is programmed / configured to remove baseline artifacts from the signal formed by the R-wave peaks. In some embodiments, such artifacts are caused by sudden baseline or RMS changes. In some embodiments, such changes are typically caused by maternal position changes. Figure 21A An exemplary data signal is shown that includes baseline artifacts.

[0187] In some embodiments, such artifacts are found using the Grubbs test for outliers, which is a statistical test performed based on the absolute deviation from the sample mean. In some embodiments, to correct for such artifacts, the change point should first be found. In some embodiments, the change point is the point (e.g., data point) at which the signal RMS or mean starts to change; such a point should satisfy the following criteria:

[0188] 1) Length (Peak Signal) - Change Point > 50

[0189] 2) prctile(peak signal (change point: end), 10) > 0.01

[0190] 3a) (where P 10％ (x) is the 10th percentile of x)

[0191] or

[0192] 3b)

[0193] In some embodiments, if the change point satisfies the above criteria, the peak signal up to that point is changed based on a statistical value defined as follows:

[0194]

[0195] Figure 21B An exemplary data signal of Figure 21A after performing baseline artifact removal according to step 1330 is shown.

[0196] With continued reference to Figure 13 , in step 1330, the exemplary inventive computing device is programmed / configured to remove outliers from the R-wave peak signal using an iterative process according to the Grubbs test for outliers. Figure 22A An exemplary R-wave peak signal including outlier data points is shown, as indicated by the diamond. In some embodiments, the iterative process of step 1330 stops when either of the following two conditions occurs:

[0197] 1)

[0198] 2) number of iterations > 4

[0199] In some embodiments, the process finds an outlier point in each iteration, and the height of such outlier point is trimmed to the median of the 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. Figure 22B An exemplary data signal of Figure 22A after performing step 1330 is shown. As can be seen, Figure 22A the outlier data points shown in Figure 22B are no longer present in In some embodiments, more outliers are revealed and removed after signal extraction, which will be described in further detail below.

[0200] With continued reference to Figure 13In step 1335, the example inventive computing device is programmed / configured to interpolate and extract R-wave signal data from each R-peak signal dataset to produce an R-wave signal channel. In some embodiments, the peak signal output from step 1330 is time interpolated to provide a signal of 4 samples per second. Figure 23A Example peak signals output from step 1330 are shown. In some embodiments, cubic spline interpolation is used to accomplish the interpolation. In some embodiments, in the case of large gaps in the interpolated data, there are high values of error, and instead, the interpolation method is a conforming piecewise cubic interpolation. In some embodiments, the conforming piecewise cubic interpolation is a piecewise cubic Hermite interpolating polynomial (“PCHIP”) interpolation. In some embodiments, after interpolation, the process of extracting the R-wave signal includes identifying other outliers in the interpolated signal. In some embodiments, in this step, further outliers are identified as one of the following:

[0201] 1) a signal peak (i.e., a peak in the interpolated R-wave signal, not a peak in the original biopotential signal) with a height greater than 1 local voltage unit and its surroundings

[0202] 2) a point located between two consecutive R-peaks that are more than 10 seconds apart

[0203] 3) the number of minutes in which a severe contact problem was found during the data checking phase (e.g., during steps 1320, 1335, and 1330)

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

[0205] Continuing with the description of step 1335, in some embodiments, after further outlier detection, signal statistics (e.g., median, minimum, and standard deviation) are computed, and a signal (e.g., a one-minute signal time window for a given channel) is identified as a corrupted signal if any of the following conditions are met:

[0206] 1) the signal still has a peak with an amplitude greater than one local voltage unit and a standard deviation greater than 0.1 after the outliers are excluded

[0207] 2) the median of the signal is greater than 0.65 and the minimum is less than 0.6

[0208] 3) more than 15% of the points comprising the signal have been deleted as outliers

[0209] Continuing with step 1335, after identifying the corrupted signal, a sliding RMS window is applied to the signal. In some embodiments, the size of the RMS window is in the range between 25 and 200 samples. In some embodiments, the size of the RMS window is 100 samples. In some embodiments, after applying the RMS window, a first order polynomial function is fit to the signal, which is then subtracted from the signal, resulting in a clean version of the interpolated signal, which can be used in subsequent steps. Figure 23B An exemplary R-wave signal after interpolation at step 1330 is shown.

[0210] Continuing with reference to Figure 13 At step 1340, the exemplary inventive computing device is programmed / configured to perform channel selection, selecting a subset of the exemplary R-wave signal channels for use in generating the electrical uterine monitoring signal. In some embodiments, at the beginning of channel selection, all channels are considered eligible candidates, and possible exclusion of channels is assessed according to:

[0211] 1) Exclude any channel that has had contact issues in more than 10% of the processing intervals so far

[0212] 2) If more than 50% of the channels are excluded based on the above, then instead exclude all channels that have had contact issues in more than 15% of the processing intervals

[0213] If the above results in exclusion of all channels, then instead, retain any channel that meets both of the following criteria, and exclude the remaining channels:

[0214] 1) The signal has a standard deviation between 0 and 0.1

[0215] 2) The signal range is less than 0.2

[0216] If the above still results in exclusion of all channels, then use only the above first condition related to standard deviation, and ignore the above second condition related to range. Figure 24A An exemplary data set including six data channels, with two data channels excluded, is shown.

[0217] In some embodiments, after removing some channels as described above, the remaining channels are grouped into pairs. In some embodiments defining channels as described above, the channel pairs are any of the above eight channels. In some embodiments, only pairs that are independent of each other (i.e., pairs that do not have a common electrode) are considered. In some embodiments, the possible pairs are as follows:

[0218] 1. Channels 1 and 2 (A1-A4 and A2-A3)

[0219] 2. Channels 1 and 5 (A1-A4 and B1-B3)

[0220] 3. Channels 1 and 6 (A1-A4 and B1-B2)

[0221] 4. Channels 1 and 7 (A1-A4 and B3-B2)

[0222] 5. Channels 2 and 5 (A2-A3 and B1-B3)

[0223] 6. Channels 2 and 6 (A2-A3 and B1-B2)

[0224] 7. Channels 2 and 7 (A2-A3 and B3-B2)

[0225] 8. Channels 3 and 5 (A2-A4 and B1-B3)

[0226] 9. Channels 3 and 6 (A2-A4 and B1-B2)

[0227] 10. Channels 3 and 7 (A2-A4 and B3-B2)

[0228] 11. Channels 3 and 8 (A2-A4 and A1-A3)

[0229] 12. Channels 4 and 5 (A4-A3 and B1-B3)

[0230] 13. Channels 4 and 6 (A4-A3 and B1-B2)

[0231] 14. Channels 4 and 7 (A4-A3 and B3-B2)

[0232] 15. Channels 5 and 8 (B1-B3 and A1-A3)

[0233] 16. Channels 6 and 8 (B1-B2 and A1-A3)

[0234] 17. Channels 7 and 8 (B3-B2 and A1-A3)

[0235] It can be seen that for each channel pair listed above, the two channels forming the pair do not share a common electrode. In some embodiments, only the valid points within a channel are used to calculate the Kendall rank correlation for each pair of channels. In some embodiments, the Kendall correlation counts the number of concordant rank signs for each pair of signals to test their statistical correlation.

[0236] In some embodiments, the channels are then selected by the following selection criteria. First, if the maximum Kendall's tau value is greater than or equal to 0.7, the selected channels are any individual channel with a Kendall's tau value in this range. However, if all selected channels were previously identified as damaged, the output signal is identified as a damaged signal. In addition, if any selected channel was previously identified as damaged, or if the range of any selected channel is greater than 0.3, any such channel is excluded from the selected channels.

[0237] Second, if no channels are selected under the first criteria above, if the maximum Kendall's tau value is greater than or equal to 0.5 but less than 0.7, the selected channels are any individual channel with a Kendall's tau value in this range. However, if all selected channels were previously identified as damaged, the output signal is identified as a damaged signal. In addition, if any selected channel was previously identified as damaged, or if the range of any selected channel is greater than 0.3, any such channel is excluded from the selected channels.

[0238] Third, if no channels are selected under the first or second criteria above, if the maximum Kendall's tau value is greater than 0 but less than 0.5, all channels with a Kendall's tau value greater than 0 are identified as selected channels. However, if the maximum correlation value is less than 0.3, the output signal is flagged as damaged and all channels with a range greater than 0.3 are excluded.

[0239] Fourth, if no channels are selected under the first three criteria above, all channels with a range greater than 0.3 and all channels with more than 15% of the deletion points are excluded, the remaining channels are selected, and the output signal is identified as a less sharpened signal, as will be discussed below with reference to step 1355.

[0240] Fifth, if no channels are selected under any of the four criteria above, all channels except those with severe contact problems are selected. However, if the number of contact problems in the selected channels exceeds 15, the output signal is flagged as damaged. Figure 24B An exemplary data set following the channel selection of step 1340 is shown.

[0241] In some embodiments, rather than selecting channels in pairs based on the correlation value of the pair, the channels are selected individually.

[0242] Continuing with reference to Figure 13In step 1345, the exemplary inventive computing device 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 for a given sampling time during a sampling interval of four samples per second for all selected channels), the 80th percentile of the signal of the selected channel is calculated according to the following formula:

[0243] Combined signal (iSample) = P 80％ (Interpolated peak signal (selected channel, iSample))

[0244] Figure 25A It shows the basis Figure 24B The selected data channel is shown to calculate the 80th percentile signal. In some embodiments, the drift baseline is then removed from the 80th percentile signal of the combination determined above to produce the EUM signal. In some embodiments, a moving average window is considered to find the baseline. In some embodiments, the moving average window is the mean during the window period subtracted from the EUM signal. In some embodiments, the window length is between 0 minutes and 20 minutes. In some embodiments, the window length is 10 minutes. Figure 25B Showing the result after removing the baseline Figure 25A An example signal.

[0245] In step 1350, the exemplary inventive computing system is programmed / configured to normalize the EUM signal calculated in step 1345. In some embodiments, normalization includes 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., preserving the original value of the extracted 80th percentile signal. Figure 26 The normalization according to step 1350 is shown. Figure 25B An example data signal following the data signal.

[0246] In step 1355, the exemplary inventive computing system is programmed / configured to sharpen the normalized EUM signal produced in step 1350, thereby producing a sharpened EUM signal. In some embodiments, sharpening is performed only on signals that were not marked as damaged in a previous step; if all relevant signals are marked as damaged, the sharpening step is not performed. In some embodiments, the goal of the sharpening step is to enhance all regions with suspicious contractions. In some embodiments, sharpening is performed as follows. First, if there is any peak in the EUM signal with a value that exceeds 200 local voltage units, the signal is marked as damaged. Second, it is determined whether the signal was previously marked as damaged. Third, the signal baseline is removed. In some embodiments, for baseline removal, if the signal duration exceeds 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 visualization voltage units. In some embodiments, the signal baseline defined in this way after the normalization step provides an EUM signal in the 0-100 range in a manner similar to the signal provided by a heart rate monitor.

[0247] Fifth, peaks are identified according to one of the following:

[0248] • If the signal was identified during step 1340 as a signal that requires less sharpening, a peak is defined as having a minimum height of 35 visualization voltage units and a minimum width of 300 samples.

[0249] • If the signal was not identified, a peak is identified as having a minimum height of 35 visualization units and a minimum width of 220 samples.

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

[0251] Peak Prominence = Peak Height - P 10％ (EUM signal)

[0252] After calculating the prominence of all peaks in the sample, each peak is eliminated if it satisfies any of the following:

[0253] • The peak has a prominence less than 12 and a height less than 40 visualization voltage units

[0254] • The peak has a prominence less than 65% of the maximum prominence of all peaks in the sample

[0255] In some embodiments, additional peaks are identified by identifying any other peaks (e.g., local maxima) with a minimum height of 15 visualization voltage units and a minimum width of 200 samples, and then eliminating all peaks with a prominence higher than 20 visualization voltage units.

[0256] As noted above, sharpening is only performed if all of the following conditions are met: (a) the signal is not corrupted (as noted above, a "corrupted" signal is identified); (b) there are no deleted points in the signal; and (c) at least one peak is identified in the preceding portion of the step. If sharpening is to be performed, then prior to sharpening, each peak is eliminated if it meets any of the following conditions:

[0257] • the prominence of the peak is less than 10 visualization voltage units

[0258] • the prominence of the peak is more than 35 visualization voltage units

[0259] • the width of the peak is more than 800 samples (i.e., 4 samples per second, 200 seconds)

[0260] After any peaks that meet one of the above conditions are eliminated, the following values are calculated for each remaining peak:

[0261] μ = mean (peak start, peak end)

[0262]

[0263] t = peak start: peak end

[0264] Once these values are calculated, a mask of zeros outside the peak region and a Gaussian function inside the peak region is created according to the following formula:

[0265]

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

[0267] Referring again to Figure 13In step 1360, post-processing is performed 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 exceeds 10 minutes, a 10-minute 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, 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 ​​higher than 100 visual voltage units are set to a value of 100 visual voltage units. Figure 28 This demonstrates applying the post-processing of step 1360 to... Figure 27B An exemplary post-processing signal generated from an exemplary sharpening signal.

[0268] After step 1360, method 1300 is completed. As described above, Figure 28 An exemplary EUM signal calculated according to method 1300 is shown. Figure 29 This illustrates representative labor map signals obtained from the same subject and during the same time period as the bioelectric data collection, based on known techniques, and the calculation of the bioelectric data. Figure 28 The EUM signal. It can be seen that... Figure 28 and Figure 29 They are essentially similar to each other and include the same peaks, which can be interpreted as indicating contractions. Therefore, it can be seen that the result of method 1300 is an EUM signal, which can be used as a signal similar to a labor map to monitor maternal uterine activity, but it can be calculated based on non-invasively recorded bioelectric potential signals.

[0269] The following examples, together with the foregoing description, illustrate some embodiments of the invention in a non-limiting manner.

[0270] Example

[0271] Figures 14A-17B Other comparative embodiments are shown between birth chart data and the output of exemplary method 200. Figure 14A , Figure 15A , Figure 16A and Figure 17A Each of the diagrams shows the labor graph signal relative to time, with contractions self-reported by the mother and monitored by the labor graph, represented by vertical lines. Figure 14B , Figure 15B , Figure 16B and Figure 17B In each of the multiple channels, the filtered R-wave signal is shown in a different color (e.g., similar to...). Figure 11Bthe calculated normalized average signal (e.g., similar to Figure 12A For comparison, Figure 14B , Figure 15B , Figure 16B and Figure 17B are shown adjacent to the corresponding one of Figure 14A , Figure 15A , Figure 16A and Figure 17A (e.g., the peak in the exemplary normalized uterine signal corresponds to the self-reported contraction). Figure 14A and Figure 14B show different data recorded for the same mother during the same time interval, as does Figures 15A to 17B As discussed above with reference to Figure 12A and Figure 12B , it can be seen that the peaks in the exemplary normalized uterine signal correspond to self-reported contractions.

[0272] A study was conducted to evaluate the effectiveness of exemplary embodiments. The study involved comparing EUM and TOCO recordings for single pregnancies with no fetal abnormalities, gestational age > 32+0 weeks, age 18-50 years, BMI < 45 kg / m 2 As described above, EUM was calculated by measuring a data sample of at least 30 minutes. Analysis of the maternal heart R-wave amplitude-based uterine activity index, referred to herein as EUM, showed promising results as an innovative and reliable method of monitoring maternal uterine activity. EUM data was highly correlated with TOCO data. Thus, EUM monitoring can provide data that is as useful as TOCO data while overcoming the drawbacks of traditional tocometry, e.g., discomfort.

[0273] Figures 18A to 27B Exemplary data present at various stages during the performance of exemplary method 1300 is shown. In particular, Figure 27A and Figure 27B show a comparison of the output signal shown by exemplary method 1300 to a tocogram signal recorded during the same time interval.

[0274] FIGS. 18A-18H Exemplary raw data received as input to exemplary method 1300 (e.g., received in step 1305) and exemplary filtered raw data produced during exemplary method 1300 (e.g., produced by step 1310) are shown. In particular, FIG. 18A , FIG. 18C , FIG. 18E and FIG. 18G show exemplary raw data, while FIG. 18B , FIG. 18D , FIG. 18F and FIG. 18HExemplary filtered data is shown. It will be apparent to those skilled in the art that, FIGS. 18A-18H represents raw and filtered biopotential data for a single channel, and in actual implementations of the method 1300 as described above, for each data channel, a data set comparable to that shown in FIGS. 18A-18H will be presented. Referring to FIG. 18A , it can be seen that there is high power line noise near sample number 6000. Referring to FIG. 18B , it can be seen that the power line noise is still high; in some embodiments, this can result in this interval being flagged as having a severe contact problem due to the relative R-wave peak energy changing from one interval to another being greater than the threshold value discussed above with reference to step 1310 of the exemplary method 1300. Referring to FIG. 18C , it can be seen that there is high power line noise near sample number 14000. Referring to FIG. 18D , it can be seen that the power line noise is still high; in some embodiments, this can result in this interval being flagged as having a severe contact problem due to the signal RMS exceeding the threshold value discussed above with reference to step 1310 of the exemplary method 1300. Referring to FIG. 18E , it can be seen that there is high power line noise throughout the signal. Referring to FIG. 18F , it can be seen that the power line noise is still high; in some embodiments, this can result in this interval being flagged as having a severe contact problem due to the SNR of this signal failing to meet the threshold SNR discussed above with reference to step 1310 of the exemplary method 1300. Referring to FIG. 18G and FIG. 18H , it can be seen that a clean signal is visible; in some embodiments, this can result in this interval not being flagged as having a contact problem.

[0275] Referring now to FIG. 19A and FIG. 19B , R-wave peaks are shown extracted according to step 1315. It will be apparent to those skilled in the art that, FIG. 19A and FIG. 19B represent R-wave peaks extracted from a single channel, and in actual implementations of the method 1300 as described above, for each data channel, a data set comparable to that shown in FIG. 19A and FIG. 19B will be presented. FIG. 19A shows filtered data (e.g., produced by step 1310) prior to performing step 1315. In FIG. 19A , the detected peak locations are indicated with asterisks. FIG. 19B shows peak data extracted after performing step 1315. In FIG. 19B , the peak locations are indicated with asterisks. It can be seen that, in FIG. 19ASome of the peak positions denoted by asterisks are not at the maximum value of the peaks in the data, and such positions are correctly denoted by the asterisks in FIG. 19B .

[0276] Reference is now made to FIGS. 20A-20C , which shows removal of EMG artifacts according to step 1320. It will be apparent to those skilled in the art that FIGS. 20A-20C denotes removal of EMG artifacts from a single channel, and in actual implementations of method 1300 as described above, a data set comparable to that shown in FIGS. 20A-20C will be presented for each data channel. FIG. 20A Exemplary filtered data (e.g., produced by step 1310) used in step 1320 is shown. FIG. 20B The same filtered data as in FIG. 20A is shown, and also includes representations of motion envelope and inter-peak sum. In FIG. 20B , peaks suspected of being corrupted are denoted by diamonds. FIG. 20C The EMG artifact-corrected corrected signal produced by step 1320 is shown. In FIG. 20C , the suspect peaks have been removed, and the corrected peaks are shown in circles, with the original peak values shown in contrasted shading.

[0277] Reference is now made to FIG. 21A and FIG. 21B , which shows removal of baseline artifacts according to step 1325. It will be apparent to those skilled in the art that FIG. 21A and FIG. 21B denote removal of baseline artifacts from a single channel, and in actual implementations of method 1300 as described above, a data set comparable to that shown in FIG. 21A and FIG. 21B will be presented for each data channel. FIG. 21A Exemplary data prior to baseline artifact removal, which can be received as input to step 1325, is shown. In FIG. 21A , baseline artifacts are denoted by circles. In FIG. 21A the data shown, the baseline ratio between the circled areas and the remaining signal is less than 0.8. In some embodiments, the corrected signal is provided by dividing the remaining signal by this factor. FIG. 21B Exemplary corrected signal such as can be produced by step 1325 is shown. In FIG. 21A , the baseline artifact areas are denoted within circles. By comparing FIG. 21A and FIG. 21B it can be seen that the baseline artifacts have been removed.

[0278] Reference is now made to FIG. 22A and FIG. 22B, showing the trimming of outliers and gaps according to step 1330. It will be apparent to those skilled in the art that, FIG. 22A and FIG. 22B represent the trimmed data set from the individual channels of outliers and gaps, and in actual implementations of the method 1300 as described above, will present a data set comparable to that shown in FIG. 22A and FIG. 22B for each data channel. FIG. 22A shows exemplary data that can be received as input to step 1330. As can be seen, the input data includes outliers near sample 450, which are represented in FIG. 22A with diamond shapes. FIG. 22B shows exemplary data of FIG. 22A after step 1330 has been performed to remove the outliers, as described above. As can be seen, the outliers shown in FIG. 22A have been removed.

[0279] Reference is now made to FIG. 23A and FIG. 23B , showing the interpolation and extraction of R-peak signals according to step 1330. It will be apparent to those skilled in the art that, FIG. 23A and FIG. 23B represent the extraction of R-peak signals from the individual channels, and in actual implementations of the method 1300 as described above, will present a data set comparable to that shown in FIG. 23A and FIG. 23B for each data channel. FIG. 23A shows exemplary R-peak signals that can be provided as output from step 1330, and received as input to step 1335. FIG. 23B shows exemplary clean interpolated R-wave signals that can be generated by performing step 1335.

[0280] Reference is now made to FIG. 24A and FIG. 24B , showing the channel selection according to step 1335. In the exemplary data sets shown in FIG. 24A and FIG. 24B , channels 3 and 8 were found to be ineligible for channel selection due to contact problems in more than 10% of the time intervals. Accordingly, in FIG. 24A and FIG. 24B only exemplary channels 1, 2, 4, 5, 6 and 7 are shown. FIG. 24A The independent channel pairs of the data shown in

[0281]

[0282]

[0283] From the above table, it can be seen that the group consisting of channels 1, 2, 4, and 7 shows a moderate correlation (e.g., a correlation greater than 0.5 but less than 0.7). Accordingly, channels 1, 2, 4, and 7 are selected in step 1340. FIG. 24B An exemplary data set output by step 1340 is shown, including the selected channels 1, 2, 4, and 7.

[0284] Referring now to FIG. 25A and FIG. 25B , the calculation of EUM signals based on the selected channels according to step 1345 is shown. The channel data shown in FIG. 24B is received as input to step 1345 in order to produce the output data shown in FIGS. 25A-25B . Referring to FIG. 25A , this figure shows the 80thpercentile signal extracted from the signal shown in FIG. 24B . FIG. 25B A corrected signal obtained by applying a drift baseline removal to the signal shown in FIG. 25A is shown.

[0285] Referring now to FIG. 26 , the calculation of normalized EUM signals according to step 1350 is shown. The corrected data produced by step 1345 and as shown in FIG. 25B is received as input to step 1350 in order to produce the normalized EUM signals as shown in FIG. 26 . FIG. 26 A normalized signal obtained by normalizing the signal shown in FIG. 25B and setting the baseline value to 30 visualized voltage units is shown. As can be seen in FIG. 26 , there are three weak peaks in the signal.

[0286] Referring now to FIGS. 27A-27C , the sharpening of the 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 in order to produce the sharpened EUM signal. FIG. 27A An exemplary normalized EUM signal produced by step 1350 is shown. FIG. 27B An exemplary enhancement mask generated according to step 1355 is shown. FIG. 27C An exemplary sharpened EUM signal produced by adding the normalized EUM signal of FIG. 27A to the mask of FIG. 27B is shown.

[0287] Referring now to FIG. 28This illustrates the post-processing of the EUM signal according to step 1360. The sharpened EUM signal generated in step 1355 is received as input to step 1360 to generate the post-processed EUM signal. FIG. 28 An exemplary post-processed EUM signal is shown after removing the drift baseline as described above with reference to step 1360. It can be seen that after the sharpening in step 1355 and the post-processing in step 1360, FIG. 26 The three weak peaks shown are in FIG. 28 It is more clearly visible in the middle.

[0288] Now for reference FIG. 29 , showing the corresponding FIG. 28 The exemplary EUM signal is a labor map signal. As previously described, it is generated according to method 1300. FIG. 28 An example EUM signal. In conjunction with the signal used to generate... FIG. 28 The exemplary EUM signal data was captured for the same subject during the same time interval. FIG. 29 The representative labor pattern signals. It can be seen that... FIG. 28 and FIG. 29 They are basically matched with each other and include three peaks that are the same in each other.

[0289] In some implementations, this is based on the use of one or more acoustic sensors (such as those referenced above). FIG. 3 The acoustic data collected by the acoustic sensor 320 is used to perform uterine monitoring. In some embodiments, the uterine monitoring process based on acoustic data is substantially similar to that described above. FIG. 13 The method 1300 described herein is a uterine monitoring process based on bioelectrical potential data, except as described below. FIG. 30 An exemplary method 3000 for uterine monitoring based on acoustic data is illustrated. In some embodiments, an exemplary inventive computing device is programmed / configured to perform method 3000. In some embodiments, the exemplary inventive computing device is programmed / configured according to method 3000 via instructions stored in a non-transitory computer-readable medium. In some embodiments, the exemplary inventive computing device includes at least one computer processor that, when the instructions are executed, becomes a specially programmed computer processor programmed / configured according to method 3000. In some embodiments, the exemplary inventive computing device is specifically configured to solve the technical problems discussed below by performing method 3000.

[0290] In step 3005, the example 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 to the abdomen of a pregnant human subject. In some embodiments, a set of raw acoustic data is received from each of two, or three, or four, or five, or six, or seven, or eight, or nine, or ten, or more number of acoustic sensors. In one specific example embodiment to be discussed in detail in the present specification of method 3000, a set of raw acoustic data is received from each of four acoustic sensors, as shown in FIG. 3

[0291] In step 3010, the example inventive computing device is specifically configured to pre-process the raw acoustic data to produce a plurality of channels of pre-processed acoustic data. In some embodiments, the pre-processing 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 a greater number of filters) to the raw acoustic data, e.g., applying each of a number X of filters to each of a number Y of channels of the raw data to produce a number X times Y of channels of pre-processed data. In some embodiments, the filters include bandpass filters. In some embodiments, the filters include DC filters. In some embodiments, the filters include finite impulse response filters or infinite impulse response (“IIR”) filters, such as Butterworth filters or Chebyshev filters or combinations thereof. In some embodiments, the filters include low-pass zero-phase lag IIR filters with a 50 Hz cutoff. In some embodiments, the filters include twelve-order Butterworth IIR filters, three-order Butterworth IIR filters, or five-order Butterworth IIR filters. In one example embodiment, the filters include five twelve-order Butterworth IIR filters with frequencies of 10-50 Hz, 15-50 Hz, 20-50 Hz, 25-50 Hz, and 30-50 Hz. In some embodiments, the five IIR filters are applied to the four raw data channels, resulting in twenty (20) channels of pre-processed data. FIG. 31A Data in example pre-processed data channels following step 3010 is shown. FIG. 31B An enlarged view of a time window of data shown in FIG. 31A

[0292] ​​In step 3015, the example inventive computing device is specifically configured to extract S1-S2 peaks from the pre-processed data channel. Those skilled in the art will appreciate that S1 and S2 refer to the first and second sounds of 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 substantially similar fashion to the R-wave peak extraction of step 1315 of method 1300 as described above. FIG. 32A Data from an example data channel is shown, with annotated R-wave peaks shown. FIG. 32B An enlarged view of a small time window of the data shown. FIG. 32A An enlarged view of a small time window of the data shown. FIG. 32C An example S1-S2 amplitude signal is shown, based on the R-wave peaks shown. FIG. 32A An example S1-S2 amplitude signal is shown, based on the R-wave peaks shown. FIG. 32D An example R-wave amplitude signal is shown, over a larger time window.

[0293] In steps 3020, 3025, and 3030, the example inventive computing device is specifically configured to remove artifacts and outliers from the data set produced in step 3015 in substantially similar fashion to that described above with reference to steps 1320, 1325, and 1330 of method 1300. It should be noted that the acoustic data analyzed by example method 3000 does not include electrical noise of the type discussed above with reference to step 1320, but can instead generally include motion-related noise recorded by the acoustic sensor. However, the process of removing such motion-related noise is substantially similar to the process of removing electrical noise described above. FIG. 33 An example data set of the example data channel after performance of steps 3020, 3025, and 3030 is shown.

[0294] In step 3035, the example inventive computing device is specifically configured to interpolate and extract S1-S2 signal data from the data set produced in step 3030 in substantially similar fashion to that described above with reference to step 1335 of method 1300. FIG. 34 Extracted S1-S2 data sets for multiple channels computed in step 3035 are shown.

[0295] In step 3040, the example inventive computing device is particularly 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 of step 3040 differs from step 1340 in one respect. As described above, some of the data channels used in step 1340 are not independent of one another due to the different nature of the biopotential sensors, and thus only some of the data channels used in step 1340 can be coupled to one another. In contrast, the acoustic sensors that collect the data used in method 3000 are single-ended, i.e., independent of one another. Thus, in step 3040, any two data channels can be appropriately coupled to one another. Thus, for example, in an implementation in which four raw data channels are processed using five different bandpass filters to produce twenty filtered data channels, there are twenty times nineteen (i.e., 380) possible channel pairs.

[0296] Following the channel selection of step 3040, in step 3045, the example inventive computing device is particularly configured to compute 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 example inventive computing device is particularly 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 example inventive computing device is particularly 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 example inventive computing device is particularly 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.

[0297] In some implementations, the output of example method 3000 is an acoustic uterine monitoring signal that is non-invasively determined by analyzing data obtainable by acoustic sensors positioned around the abdomen of a pregnant human subject. In some implementations, the acoustic uterine monitoring signal generated by example method 3000 provides similar uterine monitoring data to that generated by a tocometer and an ultrasound transducer, and can be used to monitor uterine activity, such as contractions.

[0298] FIGS. 35A-37B Examples of comparisons between tocometer data and the output of example method 3000 are shown. In FIG. 35A , FIG. 36A and FIG. 37A In each of FIG. 35B , FIG. 36B and FIG. 37BIn each example, the output of method 3000, which uses acoustic data recorded during the same time interval, is shown. It can be seen that the peaks in the exemplary acoustic-based uterine monitoring signal correspond to the peaks in the labor map data.

[0299] FIG. 38 A set of ECG-based EUM processed signals and PCG-based processed signals collected from a biopotential sensor and an acoustic sensor, respectively, are shown. In some embodiments, the processing can be based on signals collected from a biopotential sensor (e.g., as referenced above). FIG. 2 and FIG. 13 The methods discussed herein) and acoustic sensors (e.g., as referenced above) FIG. 30 The method described herein uses data collected to determine uterine monitoring signals. Example data collected from bioelectrical potential sensors or ECG-based processed EUM signals is shown in section 3801 (the first two rows represent 8 channels), and example data collected from acoustic sensors or PCG-based processed signals is shown in section 3802 (the last five rows represent 20 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. Such as FIG. 38 The signals shown can be combined or fused to generate a single uterine activity signal.

[0300] FIG. 39 A method 3900 for implementing a fusion process to generate a uterine activity signal from ECG-based processed EUM signals and PCG-based processed signals is illustrated. In some embodiments, as shown at 3901, ECG-based signals can be received in parallel from N channels, or at 3903, PCG-based signals can be received sequentially from M channels. Subsequently, machine learning techniques can be performed at 3905 to determine the 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 signals and the received PCG-based signals into a signal representing uterine activity. For example, such a signal can be generated as a weighted average of N and M channels, where each channel is associated with one signal.

[0301] In some implementations, machine learning channel weighting techniques can be implemented as gradient descent (GD) optimization processes. For example, it can be used to... FIG. 38 Each of the 28 channels represented in the diagram is assigned a weight value. The cost function for the gradient descent process can be defined as follows.

[0302] At each iteration of the gradient descent optimization process, a final signal output is determined based on a weighted average provided to each of the 28 channels. For such a final signal output, a detection algorithm based on baseline variation is used to identify contractions, which algorithm defines a start time point and an end time point for each contraction in the signal. For each identified contraction, a set of features is computed. Such a set of features can include a contraction rise time, a contraction fall time, a ratio between the contraction rise time and the contraction fall time, an SNR, a skew of the contraction, and other suitable features. Thereafter, an average value for each feature is computed across all contractions in the final signal. Then, for each feature, an optimal target value is determined (e.g., based on a normalized or optimal contraction dataset). The cost function of the GD process can correspond to a difference between the optimal target value and the average value of the feature.

[0303] In some embodiments, multiple instances of the gradient descent optimization process can be performed simultaneously, with different initial weights assigned to each channel. For example, in a first instance, all channels can be assigned an equal value or the same value. In a second instance, weights can be assigned to channels based on a quality of contraction features detected by such channels. For example, contractions can be identified by each channel, and for each contraction, a set of features can be computed, e.g., a contraction rise time, a contraction fall time, a ratio between the contraction rise time and the contraction fall time, an SNR, a skew of the contraction, and other suitable features. An average feature value can be computed across all contractions identified in a channel. Then, a weight can be assigned to the channel inversely proportional to a difference between the average feature and an optimal feature value. In a third instance, a clustering algorithm can be used to assign weights to channels. For example, for each channel, an SNR and its correlation with all other channels can be determined. Channel clusters can be defined according to the SNR of a channel and its correlation with other channels. Thereafter, each cluster can be combined into a single channel, and each combined channel can be assigned a weight based on a quality of contraction features detected by such channel. In some embodiments, the best result can be selected from the first, second, and third instances of the gradient descent optimization process described above. As described 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 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 signal according to step 3905 (e.g., determined by the method 4000 to be described below; for clarity, only some of the weights 4430 are specifically annotated in FIG. 44), and an exemplary final uterine activity signal 4440 determined according to step 3907 are shown. FIG. 44

[0304] ​In some embodiments, the process of determining channel weights is performed according to the exemplary method 4000 (e.g., the process of step 3905) shown below. FIG. 40 Step 4010 of the method 4000 shown below generates a set of electrical uterine activity signals. In some embodiments, the set of electrical uterine activity signals is generated according to step 1335 of the method 1300 shown above. FIG. 13 Step 4010 of the method 4000 shown below generates a set of electrical uterine activity signals. In some embodiments, the set of electrical uterine activity signals is generated according to step 1335 of the method 1300 shown above. FIG. 30 Step 3035 of the method 3000 shown below generates a set of acoustic uterine activity signals.

[0305] In step 4020, a plurality of channel sets is initialized. In some embodiments, each channel set includes a different combination of channels (possibly overlapping with each other). 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 a set of weights is defined only for a set of biopotential channels (e.g., electrical uterine activity channels), then all channels emitted from acoustic data will be assigned a weight of 0, and channels emitted from biopotential data will be assigned a weight based on their quality, which will be described in detail below. As will be discussed in further detail below with reference to the subsequent steps of the method 4000, after an initial selection of weights for each set, an optimization phase is performed using gradient descent and augmentation, and at the end of this process, a set of optimal, optimized weights will be selected, and the data will be weighted and averaged according to the selected set of weights. In some embodiments, the plurality of weight sets initialized in step 4020 includes four (4) weight sets. In some embodiments, the plurality of channel sets includes:

[0306] 1. “Biopotential set”, consisting only of data emitted by biopotential channels.

[0307] 2. “Acoustic set”, consisting only of data emitted by acoustic channels.

[0308] 3. “Contraction-based set”, in which channels are selected based on K-means clustering of contraction features, which will be discussed below.

[0309] 4. “Combined set”, in which all channels are considered.

[0310] In some embodiments, for the first two sets, non-zero weights are assigned to channels based on data type, as described above. In some embodiments, for the contraction-based set, initial weights are determined according to the method 4100, which will be described below with reference to FIG. 41 Method 4100.

[0311] FIG. 41A flowchart showing a method 4100 for initializing a set of contraction-based channels is shown. In step 4110, the 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 channel. In some embodiments, the contractions are identified according to method 4200, which will be described below with reference to FIG. 42 Method 4200 is described.

[0312] In step 4130, a plurality of contraction features are determined for each channel. In some embodiments, six (6) contraction features are determined for each channel. In some embodiments, the features determined for each channel include the following:

[0313] 1. Kurtosis of the signal during contractions, averaged across contractions.

[0314] 2. Relative energy: ratio of the sum of all values during contractions (sig(conts)) to the sum of all channel values, per channel:

[0315] 3. Relative time: total duration of all contractions together divided by the duration of the entire channel data.

[0316] 4. Derivative energy: ratio between the RMS of the first derivative of the signal during contractions and the RMS of the first derivative of the entire signal.

[0317] 5. Time skew: ratio between the average rise time and the average fall time, which in turn is calculated for each contraction as the difference between the start peak contraction amplitude and the peak amplitude offset.

[0318] 6. Contraction SNR: calculated as the average of two determined SNRs: global SNR and average contraction SNR. For a given channel, the global SNR is equal to the RMS of the derivative of all contraction activity divided by the RMS of the derivative of all signal outside of contractions. The average contraction SNR is equal to the average SNR across individual contractions, given by the RMS of the derivative of contraction activity divided by the RMS of the activity derivative located around a particular contraction.

[0319] In some embodiments, the output of step 4130 is a feature matrix of size N by 6, 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, a different type of clustering method (e.g., K-medoids clustering, hierarchical clustering, etc.) is used to perform clustering on the feature matrix output by step 4130. In some embodiments, the cluster with the maximum number of maxima across features is kept as the “best cluster”.

[0320] In step 4150, the best cluster is improved. In some embodiments, in an iterative process, the best cluster of channels is improved by removing channels that can reduce the internal consistency between the channels in the cluster. In some embodiments, for this purpose, an internal correlation matrix of all pairs of channels within the cluster is computed. In some embodiments, a linear correlation (e.g., Pearson correlation) is applied to compute the internal correlation matrix. In some embodiments, another correlation method is applied. In some embodiments, as a first step, candidate channels with the lowest correlation to other channels are preliminarily removed from the cluster, and the internal correlation matrix is recomputed. In some embodiments, as a result of the preliminary removal of candidate channels, a candidate channel is removed if the internal correlation improves above a predefined threshold. In some embodiments, as a second step, all channels that are not part of the best cluster are tested for cross-correlation with the average cluster signal, and if the cross-correlation is sufficiently high and has a small lag, these channels are added to the cluster. In some embodiments, the lag is computed using a cross-correlation function (which provides a cross-correlation coefficient and a lag as a one-dimensional array), where the final cross-correlation coefficient is taken as the maximum of the computed cross-correlation coefficients, and the lag is taken as the lag value corresponding to the same array element at which the final correlation coefficient is computed. In some embodiments, this computation can be represented in pseudo-code as follows:

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

[0322] corr_coef, ind_of_corr_coef = max(corr_coefs)

[0323] lag = lags[ind_of_corr_coef]

[0324] In some embodiments, the cross-correlation is sufficiently high if p > 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 seconds and 60 seconds, or between 30 seconds and 40 seconds, or between 30 seconds and 50 seconds, or between 40 seconds and 60 seconds, or between 50 seconds and 60 seconds, or less than 60 seconds, or between 25 seconds 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 internal correlation for processing during the following weight selection.

[0325] In step 4160, if an excluded channel region (e.g., a channel not selected for inclusion in the best cluster after step 4150) belongs to a high quality contraction, then such a region is considered for inclusion in the best cluster. In some embodiments, this “regional” data inclusion is performed by finding suitable data points (e.g., data points with good contraction activity) within the excluded channel, zeroing all other data points of that channel, and also including those “processed” channels with good regions as part of the best cluster. In some embodiments, good regions in other excluded channels are identified as follows: if the contraction SNR (e.g., feature #6 for each channel as described above) of a particular channel is above a threshold, then the individual SNR of the contraction of that channel is tested against a correlation threshold. In some embodiments, both the SNR and the threshold are unitless values. In some embodiments, the threshold is an arbitrary 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, if data from a high SNR contraction is from a time point with no contraction activity in a previously included “best cluster” channel, and if that data exceeds a minimum length, then that data is retained. In some embodiments, the minimum length is between 1 minute and 9 minutes, or between 2 minutes and 8 minutes, or between 3 minutes and 7 minutes, or between 4 minutes and 6 minutes, or approximately 5 minutes, or 5 minutes. In other words, if the SNR is above the threshold, and if new time points are included that are not part of a previously included contraction, and if those new time points are not too sparse, then contraction data is added to the pool of good channel data. In some embodiments, all other time points of such channels are then zeroed, and the channels with remaining “good” contraction activity are added to the best cluster. The best cluster is output in step 4170 for use as the contraction-based channel set to which weights are assigned.

[0326] Referring again to FIG. 40 In step 4030, a set of initial weights is defined for each channel set defined in step 4020. In some embodiments, two subsets of weights are defined for each channel. In some embodiments, the two subsets of weights include (1) a“channel voting” subset, as will be described below, and (2) a“Born equation” subset, in which the initial weight for each channel within a channel set is 1 / N, N equal to the number of channels within the set.

[0327] In some embodiments, the initial weights in the channel voting subset are determined as follows. The voting used to create the first subset of weights within each set is a process by which, for each data point within a channel, the number of other channels that have a contraction or do not have a contraction feature at that same data point is counted. In other words, all channels vote data point by data point on the type of activity (e.g., contraction or non-contraction) in all other channels. The votes across data points are then averaged to calculate a voting count metric for each particular channel, reflecting the degree of agreement across the set of channels with the contraction identified on that particular channel. The voting process is repeated, with each channel being voted on by all other channels.

[0328] In some embodiments, in addition to the voting, the average of two contraction scores for each channel is also calculated, where the contraction scores are determined according to step 4280 of method 4200, which will be described below. The initial weight for each channel is then calculated as the sum of the following three items:

[0329] 1. The ratio of the voting count to the number of channels in the set.

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

[0331] 3. The reliability of the contraction, calculated as the ratio between the sum of signal values during the identified contractions and the sum of the entire signal, divided by the number of contractions.

[0332] The resulting weights are then normalized so that the sum of the weights across channels is 1. In some embodiments, the uterine monitoring process analyzes received data in“frames” of the set, processing the data during a given frame and providing an output (e.g., a uterine monitoring signal) at the end of the frame. In some embodiments, the length of a frame is 10 minutes. In some embodiments, for any recorded frame that is not the first recorded frame (e.g., beyond the first 10 minutes of monitoring a given patient), the weights are averaged with the weights of the previous segment to mitigate sudden weight changes between processed segments. In some embodiments, as described above, the weight for each channel is determined to be 0.6 times the previous weight for that channel plus 0.4 times the calculated current weight for that channel.

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

[0334]

[0335] In the above expressions, E cont is the contraction energy, calculated as the sum of all temporary uterine signal values across two-thirds of the contraction width around its peak, and summed across all contractions; E tot is the sum of the entire temporary signal; A cont is the average contraction amplitude, averaged across contractions, calculated over one-third of the contraction width around its peak; R is the range (max-min) of the baseline activity amplitude between contractions; w denotes the given weight under consideration. Again, note that all weight sets have the same length, which is equal to the number of channels. The weights of channels that are not included in the channel set according to their definition (e.g., channels that originate from bio-potential signals with respect to the acoustic channel set) are equal to zero.

[0336] 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, two subsets of weights within each set compete with each other, and the best subset of weights is selected to “represent” the set. In later stages of the method 4000, the weights of different sets will compete among themselves to select the best 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 subsets of weights within each set) is based on the following metrics:

[0337] 1. Signal-to-noise ratio of the uterine activity trace.​

[0338] 2. Cost function.

[0339] 3. Shrinkage confidence measure, defined below.

[0340] 4. "Difference index".

[0341] In some embodiments, the difference index quantifies the difference between the MUA signal that has been generated so far throughout the session (i.e., from previously analyzed recording frames) and the signal that would have been generated if only the current weight set was used to analyze these previous recording segments. In some embodiments, the difference index is calculated as follows:

[0342] Difference index = 1 - max{0, r(S prev , S current )}

[0343] 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 currently selected set of channels (excluding the currently analyzed data frame, which cannot be compared to the previous data). In some embodiments, the index varies from 0 (identical traces) to 1 (any negative correlation between the two traces). This cannot be computed for the first analyzed data segment, so this step is omitted for the first analyzed data segment.

[0344] In some embodiments, the mean of the four measures between the two competing subsets is compared to a threshold and its absolute value. In some embodiments, the comparison is based on a relative difference of 10%, or, if the relative difference is not met, on an absolute difference greater than zero. If one of the two subsets exhibits a greater mean of the above four measures and relative to the threshold than the other, this subset is selected for the channel. If the subset with the better mean is below the threshold, a decision tree is activated, where each measure has a different importance in the decision process.

[0345] In some embodiments, the decision tree is based on the cost function and the confidence measure of the two measures and is executed as follows:

[0346] → If conf_1 > conf2

[0347] →→ If cost_2 is invalid (e.g., has an invalid value such as NaN or Inf), select method 1.

[0348] →→→ Else, if cost_1 is invalid, select method 2.

[0349] →→ If both cost_1 and cost_2 are valid, then the relative difference between cost_1 and cost_2 is calculated. If the relative difference supports metric 1 (meaning the cost function of method 1 is lower) by a threshold (e.g., a 10% threshold), then metric 1 is selected, otherwise, metric 2 is selected.

[0350] In the above, 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 metric of the first subset in the comparison, and conf_2 is the contraction confidence metric of the second subset in the comparison. In some embodiments, the updated uterine activity trace is finally computed using the selected optimization weights, and contractions are redefined based on the updated uterine activity trace.

[0351] In step 4060, the signals in each set of channels are enhanced. In some embodiments, two sub-steps are performed to enhance the signals. In some embodiments, the first sub-step is to enhance the data in the channels that have labeled contractions. In some embodiments, the first sub-step is performed by computing a similarity metric to the weighted average signal for each channel. To this end, three metrics are examined: (1) the correlation coefficient between the weighted average signal and each channel 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) the first parameter of a first order polynomial fit (‘slope’) between the channel data and the weighted average; and (3) the estimated error (delta) of the above fit (e.g., the difference between the first order polynomial fit and the channel). In some embodiments, the above three metrics are examined against thresholds (e.g., a threshold of 0.55 for the first metric, a threshold of 0.1 for the second metric, and a threshold of 0.3 for the third metric, all thresholds can be configured as needed), and the weights associated with any channels that exceed the thresholds are retained. The remaining weights (e.g., weights that do not cross the correlation threshold) are zeroed. The remaining weights are then scaled to sum to 1. In some embodiments, after this sub-step, an additional iteration of gradient differential optimization is run on the resulting weights. Next, in some embodiments, the weights of the traces that have higher energy than the weighted average are further amplified by a factor that is defined as the minimization of the Euclidean distance between the weighted channels and the weighted average. In some embodiments, contractions and their scores are then identified on the new weighted average uterine activity trace of each channel.

[0352] The second sub-step takes into account the weights from previous recording segments, if present. In some embodiments, the current recording segment is assigned a contribution weight (CW) from 0 to 1, and the previous weights are assigned a complementary contribution weight (1-CW). In some embodiments, the contribution weight CW assigned to a given segment with segment number N is CW = 1 / N. As more previous segments exist, the current segment will be assigned a lower CW: each additional recording segment adds 1 / segment number bits of information. The weights are then adjusted according to this weighting method to maintain a balance between the current session and previous sessions.

[0353] In step 4070, the most weighted set is determined. As described above, prior to step 4070, a weighted subset has been selected for each channel set, providing uterine activity candidates from which the most weighted set is selected to produce a final uterine activity output. In some embodiments, to select the most weighted set, the same four metrics are used as in step 4050 above. In step 4070, in contrast to step 4050, there are more than two candidates to compare and select from (e.g., four channels all have selected weighted sets as described above). Accordingly, in step 4070, the selection is performed iteratively. The selection of step 4070 begins with comparing the metrics of the first weighted set to the metrics of the second weighted set and selecting the best of these two weighted sets. The selected weighted set is compared to the third weighted set and the best weighted set is selected from this comparison, which is then compared to the fourth weighted set. The best weighted set from this comparison is then selected as the most weighted set for generating the maternal uterine activity signal. As described above, in some embodiments, the most weighted set determined in step 4070 is used as input to step 3905 of method 3900 to generate a weighted average of the data channels.

[0354] Referring now to FIG. 42 , a flowchart of a method 4200 for identifying contractions in a uterine activity signal is shown. In some embodiments, method 4200 is applied during the execution 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. In step 4210, method 4200 receives a current uterine activity channel as input. It will be apparent to those of skill in the art that while method 4200 is described with reference to a single uterine activity channel, method 4200 can also be performed on multiple uterine activity channels, including sequentially and / or simultaneously.

[0355] In step 4220, a smoothed version and an enhanced version of the input signal received in step 4210 are computed. In some embodiments, the smoothed version is computed by convolving the first derivative of the signal with a Hamming window and returning the cumulative sum of the result. In some embodiments, the convolution is performed after padding the signal with its left and right flipped versions at both ends (as will be described in further detail below with reference to step 4230), creating a continuous padded signal to ensure 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 computed by computing the hyperbolic tangent of the z-score normalized smoothed signal. In some embodiments, the enhanced version produced in this way is a smooth, rounded time series in which transient modulations in the peak amplitude of the heartbeat are easily apparent and detectable. FIG. 43 A plot showing an exemplary input uterine activity channel 4310, a smoothed channel 4320, and an enhanced channel 4330 is shown.

[0356] In step 4230, peaks are detected. In some embodiments, after flipping the signal left and right and adding it as padding to both ends, peaks of contractions, their significance, and their width are detected on the enhanced channel 4330. In some embodiments, when a mirrored signal “completes” a half-contraction with its mirrored version, it is possible to detect an incomplete contraction at the recording edge by padding the signal with its mirrored version. In some embodiments, as described above with reference to FIG. 42, the significance of a peak is computed as the ratio of the peak amplitude to the average amplitude of the signal in the vicinity of the peak. In some embodiments, the width of a peak is computed as the width of the peak at half its maximum amplitude. FIG. 13As described at step 1355 of the illustrated method 1300, peaks, significance, and width are calculated. In some embodiments, given the smoothness of the augmented channel 4330 and the above parameters, there is no clutter of detected peaks. In some embodiments, because peaks located near the edges of the signal can be missed (not accounting for padding), a second iteration of peak detection is performed in which the sensitivity of the peak finder is increased. In some embodiments, the sensitivity is increased in the second iteration by lowering the threshold applied to the width of the peaks and the distance between peaks. In some embodiments, in the first iteration, the threshold for the width of the peaks is 30 seconds and the threshold for the distance between peaks is 60 seconds, and in the second iteration, the threshold for the width of the peaks is 20 seconds and the threshold for the distance between peaks is 50 seconds. It will be apparent to those skilled in the art that threshold lowering of different magnitudes is also possible. In some embodiments, if any new peaks are detected in this second iteration, they are only considered if they are located a distance of 160 seconds or less from the edges of the signal. In some embodiments, as a third iteration, short contractions with high significance that can have been missed but constitute physiologically valid contractions are identified. In some embodiments, the sensitivity is further increased in the third iteration by lowering the threshold applied to the width of the peaks and the distance between peaks. In some embodiments, in the third iteration, the threshold for the width of the peaks is 20 seconds and the threshold for the distance between peaks is 40 seconds. In some embodiments, when determining the final contraction peaks, the width of each contraction is calculated using linear interpolation of the left and right points of the signal taken at half the peak significance.

[0357] In step 4240, outlier peaks are identified. In some embodiments, the Euclidean distance between each pair of peak significance is calculated and the error estimate for each peak is calculated as the sum of the distances to other peaks. In some embodiments, outlier error values are detected. In some embodiments, an outlier error value is a value that deviates more than 3 times the scaled median absolute deviation (MAD) from the median. In some embodiments, the scaled MAD is calculated as K*MEDIAN(ABS(A-MEDIAN(A))), where A is the values being evaluated and K is the 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. Further, 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 setting to 0.2 after normalization. In some embodiments, outlier peaks and peaks with a maximum value less than the peak height threshold are discarded as peaks.

[0358] In step 4250, incomplete contractions are detected. In some embodiments, an incomplete contraction is a contraction that is still in progress at the end of the current data segment. In some embodiments, a contraction is detected as incomplete if the contraction peak and contraction offset interval is less than a required minimum time, and if the activity level before and after the contraction differs from a certain threshold. In some embodiments, the required minimum time is one 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 so labeled for further use. For example, in some embodiments, when calculating the overall quality of the trace, contractions that are labeled as incomplete are given less weight (e.g., are weighted by a factor of 0.5 in determining the overall SNR of the trace). In some embodiments, incomplete contractions are completed by considering data from subsequent segments.

[0359] In step 4260, a confidence metric is calculated for each contraction. In some embodiments, three (3) confidence metrics are calculated. In some embodiments, the confidence metrics are calculated as follows:

[0360] 1. Contraction relative energy. This is calculated by first calculating the contraction energy as the sum of data points spanning two-thirds of the contraction width around its peak, and then dividing its energy by the sum of the energy of all contractions.

[0361] 2. Ratio between the upper third average activity (e.g., the average amplitude of the contraction width's third 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 are not located within a contraction and finding the difference between the 5th and 95th percentiles of that distribution.

[0362] 3. Ratio between the range of values during contractions and the range of values between contractions. The range of values between contractions is as calculated in the immediately preceding item 2. The range of values of each contraction is calculated by taking the distribution of data points that constitute the contraction data and finding the difference between the 5th and 95th percentiles of that distribution.

[0363] In step 4270, noisy contractions and small contractions are eliminated based on the confidence metrics determined in step 4260. In some embodiments, the confidence metrics are compared to a predefined threshold and used to eliminate noisy contractions. In some embodiments, a noisy contraction is a contraction whose confidence metric is less than or equal to a predefined threshold. In some embodiments, the predefined threshold is 0.5. Small contractions are also eliminated in some embodiments. In some embodiments, a small contraction is a contraction whose normalized peak is less than 0.2 (e.g., less than 0.2 times the maximum peak contraction).

[0364] In step 4280, a contraction score is determined for each contraction. In some embodiments, the contraction score is used, for example, to determine an initial weight 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 between the average activity level before and after a given contraction normalized by the contraction peak amplitude. In some embodiments, a second contraction score is calculated as a normalized significance calculated as the difference between the contraction peak amplitude and the mean of the activity surrounding the contraction divided by the peak amplitude. FIG. 40

[0365] As described herein, one technical problem in the field of maternal / fetal care is that existing solutions for monitoring uterine activity (e.g., contractions) by using a tocodynamometer and an ultrasound sensor require a pregnant woman to wear an uncomfortable sensor and can produce unreliable data when worn by an obese pregnant woman (e.g., the sensor can not have sufficient sensitivity to produce usable data). As discussed further herein, exemplary embodiments present a technical solution to this technical problem by using various sensors (e.g., biopotential sensors and / or acoustic sensors) integrated into a comfortable wearable device and analyzing data obtainable by such sensors (e.g., electrodes and / or acoustic sensors) to produce signals that can be used to monitor uterine activity. Another technical problem in the field of maternal / fetal care is that existing solutions for analysis based on signals obtainable by sensors (e.g., biopotential sensors and / or acoustic sensors) integrated into a comfortable wearable device are limited to analyzing such signals to extract cardiac data. As discussed herein, exemplary embodiments present a technical solution to this technical problem by analyzing biopotential data and / or acoustic data to produce signals that can monitor uterine activity (e.g., contractions).

[0366] The entire contents of the publications cited herein are incorporated herein by reference. While various aspects of the application have been illustrated and described, it will be understood by those skilled in the art that various changes can be made, and substitutions can be made by equivalents, without departing from the scope of the application as recited in the claims. In addition, many modifications can be made to adapt a particular situation to the teachings of the present disclosure without departing from the central inventive concept described herein. Finally, it is to be understood that the terminology employed herein is used for the purpose of describing particular embodiments only and is not intended to be limiting since the scope of the present application will be defined by the claims as interpreted in light of this disclosure.​

Claims

1. A computer-implemented method comprising: receiving, by at least one computer processor: channel data for a plurality of electrical uterine contraction monitoring signal channels, wherein the channel data for the plurality of electrical uterine contraction monitoring signal channels is different from raw biopotential input; and channel data for a plurality of acoustic uterine contraction monitoring signal channels, wherein the channel data for the plurality of acoustic uterine contraction monitoring signal channels is different from raw acoustic input; computing, by the at least one computer processor, a plurality of channel weights for the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels using a machine learning algorithm, wherein each of the channel weights corresponds to a particular one of the electrical uterine contraction monitoring signal channels or a particular one of the acoustic uterine contraction monitoring signal channels; computing, by the at least one computer processor, a weighted average of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels based at least on the plurality of channel weights; and generating, by the at least one computer processor, a combined uterine contraction monitoring signal channel based on the weighted average.

2. The computer-implemented method of claim 1, wherein the machine learning algorithm comprises a gradient descent optimization process.

3. The computer-implemented method of claim 1, wherein the machine learning algorithm comprises the steps of: defining, by the at least one computer processor, a plurality of channel sets, each of the plurality of channel sets comprising 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, wherein each of the plurality of initial weight sets corresponds 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, wherein each of the plurality of optimized weight sets corresponds to a particular one of the plurality of channel sets; and selecting, by the at least one computer processor, a particular one of the plurality of optimized weight sets as the plurality of channel weights.

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

5. The computer-implemented method of claim 3, wherein the step of selecting a particular one of the plurality of optimized weight sets is performed by a process comprising: generating, by the at least one computer processor, a plurality of temporary uterine activity traces, wherein each of the plurality of temporary uterine activity traces corresponds to a particular one of the plurality of optimized weight sets; computing, by the at least one computer processor, for each of the plurality of optimized weight sets, a signal-to-noise ratio, a cost function, a contraction confidence metric, and a difference index for a particular one of the temporary uterine activity traces that corresponds to the particular one of the plurality of optimized weight sets. ​ computing, by the at least one computer processor, a first mean value for the first temporary uterine activity trace, the first mean value being a mean of a signal-to-noise ratio of the first temporary uterine activity trace, a cost function of the first temporary uterine activity trace, a shrinkage confidence measure of the first temporary uterine activity trace, and a diversity index of the first temporary uterine activity trace; and selecting, by the at least one computer processor, the first one of the plurality of weight sets as the best weight set for the particular one of the plurality of channels based on a determination that the first mean value is greater than the second mean value.

6. The computer-implemented method of claim 3, further comprising: generating, by the at least one computer processor, a first temporary uterine activity trace and a second temporary uterine activity trace corresponding to a particular one of the plurality of channels, wherein the first temporary uterine activity trace corresponds to a first one of the plurality of weight sets for the particular one of the plurality of channels, and wherein the second temporary uterine activity trace corresponds to a second one of the plurality of weight sets for the particular one of the plurality of channels; computing, by the at least one computer processor, a signal-to-noise ratio, a cost function, a shrinkage confidence measure, and a diversity index for the first temporary uterine activity trace for the first one of the plurality of weight sets; computing, by the at least one computer processor, a signal-to-noise ratio, a cost function, a shrinkage confidence measure, and a diversity index for the second temporary uterine activity trace for the second one of the plurality of weight sets; computing, by the at least one computer processor, a first mean value for the first one of the plurality of weight sets, the first mean value being a mean of a signal-to-noise ratio of the first temporary uterine activity trace, a cost function of the first one of the plurality of weight sets, a shrinkage confidence measure of the first one of the plurality of weight sets, and a diversity index of the first one of the plurality of weight sets; computing, by the at least one computer processor, a second mean value for the second one of the plurality of weight sets, the second mean value being a mean of a signal-to-noise ratio of the second temporary uterine activity trace, a cost function of the second one of the plurality of weight sets, a shrinkage confidence measure of the second one of the plurality of weight sets, and a diversity index of the second one of the plurality of weight sets; selecting, by the at least one computer processor, the first one of the plurality of weight sets as the best weight set for the particular one of the plurality of channels based on a determination that the first mean value is greater than the second mean value; and selecting, by the at least one computer processor, the second one of the plurality of weight sets as the best weight set for the particular one of the plurality of channels based on a determination that the second mean value is greater than the first mean value.

7. The computer-implemented method of claim 6, further comprising: enhancing the plurality of channel sets by the at least one computer processor by enhancing data in channels with labeled contractions.

8. The computer-implemented method of claim 3, wherein the step of defining the plurality of channel sets comprises defining a contraction-based channel set, and wherein the contraction-based channel set is determined by a process comprising: 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; and selecting, by the at least one computer processor, one of the plurality of clusters as the contraction-based channel set.

9. The computer-implemented method of claim 8, wherein the step of defining the plurality of channel sets further comprises: improving, by the at least one computer processor, the one of the plurality of clusters by removing channels in the one of the plurality of clusters that decrease intra-cluster consistency between the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels.

10. The computer-implemented method of claim 8, wherein 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 the plurality of acoustic uterine contraction monitoring signal channels that is not included in the one of the plurality of clusters.

11. The computer-implemented method of claim 8, wherein the step of identifying the 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 is performed by a process comprising, for each of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels: generating, by the at least one computer processor, an enhanced version of the one of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels by computing a hyperbolic tangent of a z-score normalized smoothed signal; detecting, by the at least one computer processor, a set of candidate contractions in the enhanced one of the plurality of electrical uterine contraction monitoring signal channels and the plurality of acoustic uterine contraction monitoring signal channels, wherein the set of candidate contractions comprises a plurality of candidate contractions; computing, by the at least one computer processor, a plurality of confidence metrics for each of the candidate contractions; and removing, by the at least one computer processor, at least one of the candidate contractions from the set of candidate contractions based on a confidence metric corresponding to the removed candidate contraction, thereby producing the set of contractions.

12. The computer-implemented method of claim 3, wherein the step of defining, by the at least one computer processor, the plurality of initial weight sets comprises: generating, by the at least one computer processor, a set of channel vote weights and a set of Borne equation weights for each of the set of channels.

13. The computer-implemented method of claim 1, wherein the step of receiving, by the at least one computer processor, the channel data for the plurality of electrical uterine contraction monitoring signal channels comprises generating at least one of the plurality of electrical uterine contraction monitoring signal channels, and wherein the at least one of the plurality of electrical uterine contraction monitoring signal channels is generated by a process comprising: receiving, by the at least one computer processor, the raw biopotential inputs, wherein each of the raw biopotential inputs is received from a corresponding one of a plurality of electrodes, wherein each of the plurality of electrodes is positioned so as to measure a respective one of the raw biopotential inputs of a pregnant human subject; generating, by the at least one computer processor, a plurality of signal channels from the plurality of raw biopotential inputs, wherein the plurality of signal channels from the plurality of raw biopotential inputs comprises at least three signal channels; preprocessing, by the at least one computer processor, respective channel data for each of the signal channels to produce a plurality of preprocessed signal channels, wherein each of the preprocessed signal channels comprises respective preprocessed channel data; extracting, by the at least one computer processor, a respective plurality of R- wave peaks from the preprocessed channel data for each of the preprocessed signal channels to produce a plurality of R-wave peak data sets, wherein each of the R-wave peak data sets comprises 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: at least one signal artifact or at least one outlier data point, wherein the at least one signal artifact is one of an electromyography 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 respective 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, were removed, to produce a plurality of interpolated R-wave peak data sets; generating, by the at least one computer processor, for each respective interpolated R-wave peak data set, a respective R-wave signal data set for a respective R-wave signal channel at a predetermined sampling rate, to produce 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 respective R-wave signal data set of the at least one first specific R-wave signal channel and a respective R-wave signal data set of the at least one second specific R-wave signal channel; and generating, by the at least one computer processor, electrical uterine contraction monitoring data representing an electrical uterine contraction monitoring signal based on at least a respective R-wave signal data set of the first selected R-wave signal channel and a respective R-wave signal data set of the second selected R-wave signal channel, thereby producing the at least one electrical uterine contraction monitoring signal channel.

14. The computer-implemented method of claim 1, wherein the step of receiving, by the at least one computer processor, the signal channel data of the plurality of acoustic uterine contraction monitoring signal channels comprises generating at least one of the plurality of acoustic uterine contraction monitoring signal channels, and wherein the at least one of the plurality of acoustic uterine contraction monitoring signal channels is generated by a process comprising: receiving, by the at least one computer processor, the raw acoustic inputs, wherein each raw acoustic input is received from a corresponding acoustic sensor of a plurality of acoustic sensors, wherein each of the plurality of acoustic sensors is positioned so as 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 plurality of raw acoustic inputs, wherein the plurality of signal channels from the plurality of raw acoustic inputs comprises at least three signal channels; pre-processing, by the at least one computer processor, respective signal channel data of each signal channel to produce a plurality of pre-processed signal channels, wherein each pre-processed signal channel comprises 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 of each of the pre-processed signal channels to produce a plurality of S1-S2 peak data sets, wherein each S1-S2 peak data set comprises 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: at least one signal artifact or 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, 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 respective 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, were removed, to produce 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 a respective S1-S2 signal channel at a predetermined sampling rate based on each respective interpolated S1-S2 peak data set to produce 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 the respective S1-S2 signal data set of the at least one first particular S1-S2 signal channel and the respective S1-S2 signal data set of the at least one second particular S1-S2 signal channel; and generating, by the at least one computer processor, acoustic uterine contraction monitoring data representing an acoustic uterine contraction monitoring signal based on at least the respective S1-S2 signal data set of the first selected S1-S2 signal channel and the respective S1-S2 signal data set of the second selected S1-S2 signal channel to produce the at least one acoustic uterine contraction monitoring signal channel.

Images