Rémi Souriau
email: remi.souriau@univ-grenoble-alpes.fr

Julie Fontecave-Jallon

Bertrand Rivet
email: bertrand.rivet@grenoble-inp.fr

Fetal heart rate monitoring by fusion of estimations from two modalities: a Modified Viterbi's algorithm

Keywords:

Objective. In the context of fetal heart rate (fHR) estimation, the article addresses the fusion of two monomodal estimations into a multimodal one. Electrical and mechanical modalities are considered through the use of abdominal electrocardiogram (ECG) and phonocardiogram (PCG). The aim of the fusion is to provide a fHR estimation, robust against noise sources, and especially against the risk of confusion with the mother heart rate (mHR). Approach. A hidden Markov chain is considered to model variations of the two monomodal fHR estimations and the true fHR. Thanks to the Viterbi's algorithm (VA), the two fHR estimations are then merged. However, the classical VA does not ensure that the merged estimation always follows the true fHR: it may alternate between the fHR and the mHR, especially when the mother component is not correctly removed from the abdominal ECG. Therefore, a modified VA is proposed to efficiently avoid confusion with mHR. The aim is to discourage large variations between successive states. Results. Comparisons between classical and modified VA are performed on real pregnant women data. The modified VA reduces the confusion between fHR and mHR in major cases compared to the classical VA. For most recordings, the confusion with mHR decreases under 1%, and for the best case, it is reduced from 59% to 0%. Conclusion. The modified VA succeeds to improve the fHR estimation by reducing confusion with mHR. Significance. Fusion of estimations from two modalities is a promising approach for more robust fHR monitoring.

Introduction

The fetal heart rate (fHR) estimation is a key element in monitoring the health status of the fetus during pregnancy. Study of fHR and its variability may highlight fetal distress that can lead to the need of an emergency delivery. Nowadays, the cardiotocography (CTG) is used as the clinical reference for the fHR monitoring, as recommended by the International Federation of Gynecology and Obstetrics [START_REF] Ayres-De Campos | FIGO consensus guidelines on intrapartum fetal monitoring: Cardiotocography[END_REF]. The CTG is based on Doppler ultrasounds and monitors the fHR as well as the uterine contractions during labor. However, the use of CTG presents some limitations like signal loss or confusion between the fHR and the mother heart rate (mHR) [START_REF] Bhogal | Maternal and fetal heart rate confusion during labour[END_REF]. Some studies show that the use of CTG increases the cesareans rate [START_REF] Buisson | Confusion des enregistrements cardiaques foetal et maternel au cours de l'expulsion : impact sur l'état néonatal et le taux d'extraction instrumentale[END_REF][START_REF] Devane | Cardiotocography versus intermittent auscultation of fetal heart on admission to labour ward for assessment of fetal wellbeing[END_REF]. These limitations justify the interest of research for alternative techniques for fetal monitoring. One can cite methods like the fetal scalp electrodes [START_REF] Feinstein | Fetal Heart Monitoring: Principles and Practices[END_REF][START_REF] Westgate | Suitability of fetal scalp electrodes for monitoring the fetal electrocardiogram during labour[END_REF] or the fetal magnetocardiography [START_REF] Kariniemi | The fetal magnetocardiogram[END_REF][START_REF] Oweis | Non-invasive fetal heart rate monitoring techniques: Review article[END_REF]. These methods can provide higher quality measures but present also different disadvantages (highly invasive, expensive, etc). Non-invasive methods, much more ergonomic for clinical uses, are of high interest for fHR monitoring. Among them, electrocardiography [START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF] allows to measure the electrical activity of the fetal heart by placing electrodes on the mother abdomen. The obtained abdominal electrocardiogram (ECG) signals are corrupted by maternal heart activity and require a preprocessing step to extract the fetal ECG (ECGf). ECGf is widely studied in the literature: many methods have been proposed to remove noises [START_REF] Kahankova | A review of signal processing techniques for non-invasive fetal electrocardiography[END_REF], and there are many databases of synthetic and real ECG recordings available to test algorithms (see [START_REF] Andreotti | An open-source framework for stress-testing non-invasive foetal ECG extraction algorithms[END_REF][START_REF] Matonia | Fetal electrocardiograms, direct and abdominal with reference heartbeat annotations[END_REF] for examples). Another non-invasive approach is the phonocardiography [START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF] which consists in measuring the mechanical activity of the heart using a microphone. The use of the phonocardiogram (PCG) for fetal monitoring regains in interest the past recent years [START_REF] Prashanth Chetlur | Trends in fetal monitoring through phonocardiography: Challenges and future directions[END_REF]. It gives access to fetal PCG (PCGf) and fetal heart sounds, from which fHR can be estimated. In the present study, the ECG and the PCG are used simultaneously to estimate the fHR.

Fig. 1a gives a heart scheme and Fig. 1b is an illustration of ECG and PCG signals. At each cardiac cycle, the passing electric current through the heart (or myocardium) generates P,Q,R,S,T waves during contractions and relaxations of the atria and ventricles. The opening and closing of valves during the systole and the diastole product sounds which can be observed in the PCG signal and are noted 𝑆 1 and 𝑆 2 .

From ECG and PCG signals, there are different methods based on R peaks detection [START_REF] Holsinger | A QRS preprocessor based on digital differentiation[END_REF][START_REF] Mallat | Singularity detection and processing with wavelets[END_REF] or 𝑆 1 and 𝑆 2 sounds detection [START_REF] Ismail | Localization and classification of heart beats in phonocardiography signals-a comprehensive review[END_REF] to estimate the heart rate of an adult. However, for fHR estimation, a first step consists in extracting the ECGf and the PCGf from, respectively, the abdominal ECG (ECGa) and the abdominal PCG (PCGa). This step is not easy due to noise sources in abdominal signals. This is especially true for the ECGa because the mother heart activity is also present and the signal power of the mother is significantly greater than the fetal heart activity in most cases. Different methods to extract the ECGf have been proposed, like template subtraction [START_REF] Vullings | Dynamic segmentation and linear prediction for maternal ECG removal in antenatal abdominal recordings[END_REF], source separation [START_REF] Da | Blind source separation (BSS) of mixed maternal and fetal electrocardiogram (ECG) signal: A comparative study[END_REF][START_REF] Sameni | Multichannel electrocardiogram decomposition using periodic component analysis[END_REF] or neural network [START_REF] Zhong | Fetal electrocardiography extraction with residual convolutional encoder-decoder networks[END_REF]. In order to be robust against noise sources, most of these methods require to add redundancy of information by using more than a single sensor on the woman abdomen (three to four ECGa in [START_REF] Da | Blind source separation (BSS) of mixed maternal and fetal electrocardiogram (ECG) signal: A comparative study[END_REF], five in [START_REF] Sameni | Multichannel electrocardiogram decomposition using periodic component analysis[END_REF] for example). However, such a solution is not ergonomic for the clinician and the mother in real conditions where emergency surgery can be needed.

Therefore, an alternative solution to estimate the fHR from a few numbers of sensors is to consider multimodality. Its interest lies on the complementary and redundant information from modalities to gain in robustness against noise sources. Multimodality has been identified as a promising approach for fetal monitoring [START_REF] Sameni | A review of fetal ECG signal processing; issues and promising directions[END_REF]. For example, the study of [START_REF] Noorzadeh | A multi-modal approach using a non-parametric model to extract fetal ecg[END_REF] used an ECGa signal coupled with a PCGa signal to extract the ECGf. However, only a very few studies focused on fHR estimation from multiple modalities. In [START_REF] Shokouhmand | Multi-modal framework for fetal heart rate estimation: Fusion of low-snr ecg and inertial sensors[END_REF], four ECG signals are combined with one seismocardiographic recording and one gyrocardiographic recording for fHR estimation. The use of PCG for fHR estimation is considered in [START_REF] Khandoker | Validation of beat by beat fetal heart signals acquired from four-channel fetal phonocardiogram with fetal electrocardiogram in healthy late pregnancy[END_REF] and the authors used the ECG modality to validate their fHR.

In a previous study [START_REF] Souriau | Fusion of estimations from two modalities using the Viterbi's algorithm: application to fetal heart rate monitoring[END_REF], we investigated a preliminary approach of fusion of fHR data from ECG and PCG for fHR estimation. By considering the reference fHR as a latent dynamic process and the two estimations as observations processes, a bimodal hidden Markov model (HMM) is used to model these processes. In the HMM, the underlying laws between the variables are probability laws. The reference fHR is assumed to be provided by the CTG (recorded in parallel to the ECG and the PCG). And the probability laws are modelled according to some recording sessions. The Viterbi's algorithm (VA) allows then to decode the true fHR from the two fHR estimations and probability laws.

However, because of the possible presence of the mHR in the fHR estimation inputs, the VA can fail to estimate the fHR. In this current paper, a modified version of the VA is then proposed to improve robustness against this maternal noise source.

Monomodal estimation of the fetal heart rate

The steps of the complete process of fHR estimation are summarized on Fig. 2b. After description of the clinical protocol (Sec. 2.1), this section summarizes the different steps preceding the fusion for multimodal estimation, which will be described in Sec. 3. Before this fusion, the monomodal fHR estimation (Sec. 2.2) using each single modality (ECG or PCG) is necessary and separated into two steps: (1) the preprocessing (Sec. 2.2.1) to extract, on one hand, gives the methodology of fHR estimation from the ECG and the PCG modality. After pre-processing, the monomodal fHR estimation step is performed on the ECG and the PCG separately. This step similar for the two modalities is detailed in [START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF][START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF] and summarized in Sec. 2.2. The fusion step is detailed in Sec. 3.

the ECGf from the ECGa and the thoracic ECG (ECGt) and, on the other hand, to extract the PCGf from the PCGa, and (2) the fHR estimation method (Sec. 2.2.2) which is similar for both modalities. It is worth noting that all the steps from abdominal recordings (ECGa and PCGa) to monomodal fHR estimations (fHR (ecg) and fHR (pcg) ) are recalled for clarity of the paper and do not contain new contribution.

Clinical protocol

Real data considered in this paper comes from six volunteer pregnant women older than 18 years old, between 38(+0) and 38(+4) weeks(+days) of gestation with no maternal or fetal complication, as part of a clinical protocol established at the University Hospital of Grenoble (study No RCB: 2018-A03182-53) [START_REF] Faisant | Non-invasive fetal monitoring: Fetal Heart Rate multimodal estimation from abdominal electrocardiography and phonocardiography[END_REF]. After signature of informed consent, the volunteer was laying on her back, in a comfortable position, in order to minimize movements and electrical interference. As illustrated on Fig. 2a, recording sessions consisted in the acquisition of synchronous signals (Powerlab acquisition system, ADInstruments, sampling frequency 1KHz):

• one ECGt with two thoracic bipolar electrodes on the maternal heart axis and a reference electrode on the maternal wrist (BioAmp, ADInstruments). • one ECGa with two abdominal bipolar electrodes across the axis of the fetal heart (defined after clinical auscultation) and the same reference electrode as for the thoracic ECG (BioAmp, ADInstruments). • one PCGa with a cardio-microphone (MLT201, ADInstruments) put on the skin of the maternal abdomen as close as possible to the fetal heart.

In addition, the fHR was also directly estimated using the reference CTG (Avalon F20/F30, Philips) [START_REF] Ayres-De Campos | FIGO consensus guidelines on intrapartum fetal monitoring: Cardiotocography[END_REF]. The sampling frequency of fHR measure from CTG has been set to 4Hz. The duration of each recording session was about 30 minutes in rest condition for the mother.

Methodology of monomodal fHR estimation

This section describes the two steps of the methodology for monomodal fHR estimation using a single modality (ECG or PCG).

Preprocessing

The preprocessing step consists in attenuating the noise sources of the ECGa and the PCGa in order to obtain an approximation of the ECGf and the PCGf [START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF][START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF]. An approximation of these signals is sufficient in the context of our proposed fHR estimation, since this latter is not based on events detection. For the ECGa, the baseline is filtered with a high-pass filter (finite impulse response(FIR), order 1024, cut-off frequency at 10 Hz) to enhance the fetal (a) Estimation of fHR (ecg) (denoted 𝑌 (1) [0∶𝑇 ] in Sec. 3.1). (b) Estimation of fHR (pcg) (denoted 𝑌 (2) [0∶𝑇 ] in Sec. 3.1). R peaks. Then, the power-line inference is removed using a band-stop filter (Butterworth filter, order 2, stop band [49-51] Hz) and a low pass filter (FIR, order 100, cut-off frequency 80 Hz) is applied to attenuate high frequency noise [START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF]. Finally, the mother heart activity in the ECGa is attenuated using Approximate Linear Dependency Kernel Recursive Least-Squares (ALD-KRLS) [START_REF] Van Vaerenbergh | A comparative study of kernel adaptive filtering algorithms[END_REF] with the ECGt signal as reference (see [9, figure 2] for an example of mother ECG component removal). Considering PCGa, it is first filtered using a band-pass filter (FIR, order 100, pass band Hz) to focus on the frequency range of the heart sounds and get the PCGf signal. The absolute value is then filtered with a low pass filter (FIR, order 100, cut-off frequency 15 Hz) to highlight the envelop of the PCGf [START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF]. The signal output is therefore not exactly the PCGf but a smooth signal for which the periodic aspect is well characterized and from which fHR can be more easily estimated.

Monomodal fHR estimation

After preprocessing steps, fHR estimations are carried out, separately from the ECGf and the PCGf signals. This monomodal fHR estimation step is exactly the same for both modalities, as illustrated on Fig. 3. First, spectrograms of the ECGf and the PCGf are computed with a window (hamming) of 4s and shifted by 250ms (Fig. 3, first line). Choosing a window of 4s allows to highlight the frequential harmonic structure of the ECGf and the PCGf. A shift of 250ms allows, at the end, a fHR estimation at the same sampling frequency than the CTG (i.e. 4 Hz).

Let 𝑍 be the spectrogram matrix from one modality. Using the source-filter model, 𝑍 can be decomposed as 𝑍 = 𝑍 (𝑒) ⊙ 𝑍 (𝜙) where 𝑍 (𝑒) and 𝑍 (𝜙) are, respectively, the excitation and the filter matrices of the spectrogram [START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF][START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF]. The excitation matrix contains the heart rate information and the filter matrix stores the waveform of each heartbeat. ⊙ is the Hadamard product. 𝑍 (𝑒) and 𝑍 (𝜙) are factorized using the NMF algorithm [START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF][START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF]: 𝑊 (𝑒) is a dictionary of Dirac combs whose fundamental frequencies covered the frequency range of the fHR, and the time evolution of the fHR is given by 𝐻 (𝑒) (Fig. 3, second line). The vector of fHR estimations (Fig. 3, third line) is finally extracted from the matrix 𝐻 (𝑒) : following [START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF][START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF], at each time 𝑗, the fHR estimation is given by the frequency associated to the maximal value in column 𝑗 of 𝐻 (𝑒) .

𝑍 = ( 𝑊 (𝑒) 𝐻 (𝑒) ) ⊙ ( 𝑊 (𝜙) 𝐻 (𝜙) ) . ( 1
)
The two monomodal fHR estimations (Fig. 3a from ECG modality and Fig. 3b for PCG modality) are denoted fHR (ecg) and fHR (pcg) , and are used as inputs of the next step, which is the fusion step.

Multimodal estimation of the fetal heart rate

This section gathers the main contribution of this paper and is dedicated to the fusion step where fHR (ecg) and fHR (pcg) are merged into a new estimation denoted fHR (fusion) .

The considered hidden Markov model (HMM) to model the fHR evolution is first presented in (Sec. 3.1). Then, the Viterbi's algorithm (VA) is described in (Sec. 3.2) and its limits in (Sec. 3.3). Finally, a modified version of VA is proposed in (Sec. 3.4) to overcome limits of the classical VA.

Hidden Markov models

The HMM is used in this paper to model the evolution of fHR time series. The simplest case of HMM is given in Fig. 4. where 𝑋 0∶𝑇 = { 𝑋 𝑘 } 𝑘∈ 0,𝑇 is the latent path and 𝑌 0∶𝑇 = { 𝑌 𝑘 } 𝑘∈ 0,𝑇 is the observation path (where 𝑘 is the time index). The time evolution of the latent path is characterized by a probability law called the transition law and the probability rule of the observation state given the latent state is the emission law.

Fig. 5 gives the model used in this study. Let denote:

• 𝑋 0∶𝑇 = { 𝑋 𝑘 } 𝑘∈ 0,𝑇 the true fHR; • 𝑋 (1) 0∶𝑇 = { 𝑋 (1) 𝑘 } 𝑘∈ 0,𝑇 and 𝑋 (2) 0∶𝑇 = { 𝑋 (2) 𝑘 } 𝑘∈ 0,𝑇
the fHR according to, respectively, the first modality and the second modality;

• 𝑌 (1) 0∶𝑇 = { 𝑌 (1) 𝑘 } 𝑘∈ 0,𝑇 and 𝑌 (2) 0∶𝑇 = { 𝑌 (2) 𝑘 } 𝑘∈ 0,𝑇
the estimation of the fHR from both modalities: fHR (ecg) and fHR (pcg) .

All the fHRs and their variations are modelled in the complete HMM (Fig. 5a). 𝑋 0∶𝑇 refers to the true fHR, 𝑋 (1) 0∶𝑇 and 𝑋 (2) 0∶𝑇 are, respectively, the representation of the fHR in the electrical modality and the mechanical modality. And 𝑌 (1) 0∶𝑇 and 𝑌 (2) 0∶𝑇 are, respectively, fHR estimations from these modalities using the ECGf and the PCGf. In practice, the complete HMM cannot be used because 𝑋 (𝑖) 0∶𝑇 (𝑖 ∈ {1, 2}) as well as probability rules Pr

( 𝑋 (𝑖) 𝑘 |𝑋 𝑘
) for all 𝑘 are unknown. Because, for each 𝑘 and 𝑖 ∈ {1, 2}, there is no temporal link between 𝑋 (𝑖) 𝑘+1 and 𝑋 (𝑖) 𝑘 , the complete HMM can be reduced to the HMM given in Fig. 5b where emissions from 𝑋 𝑘 to 𝑋 (𝑖) 𝑘 and 𝑋 (𝑖) 𝑘 to 𝑌 (𝑖) 𝑘 are modelled with an emission matrix for both modalities.

Let consider the HMM given in Fig. 5b. The probability law between two successive latent states (the transition rule) and the probability laws between current latent state and observations (the emission rules) are noted:

∀𝑘 ≥ 1, ∀(𝑖, 𝑗) ∈ 1, 𝑁 2 , ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑎 𝑖𝑗 = Pr ( 𝑋 𝑘 = 𝑗 | | | 𝑋 𝑘-1 = 𝑖 ) 𝑏 (1) 𝑖𝑗 = Pr ( 𝑌 (1) 𝑘 = 𝑗 | | | 𝑋 𝑘 = 𝑖) ) 𝑏 (2) 𝑖𝑗 = Pr ( 𝑌 (2) 𝑘 = 𝑗 | | | 𝑋 𝑘 = 𝑖 ) , ( 2
)
where 𝑁 is the number of states which are directly written with index 𝑖 and 𝑗 for simplification. Coefficients

𝐴 = ( 𝑎 𝑖𝑗 ) 1≤𝑖,𝑗≤𝑁 , 𝐵 (1) = ( 𝑏 (1) 𝑖𝑗 ) 1≤𝑖,𝑗≤𝑁
and

𝐵 (2) = ( 𝑏 (2) 𝑖𝑗 ) 1≤𝑖,𝑗≤𝑁
form, respectively, the transition matrix and the emission matrices. The initial state 𝑋 0 is defined by its probability Pr(𝑋 0 ). Both modalities are assumed to be subjected to independent noise sources, so the joint probability of observations given the state satisfies

∀(𝑖, 𝑗, 𝑚) ∈ 1, 𝑁 3 , Pr ( 𝑌 (1) 𝑘 = 𝑗, 𝑌 (2) 𝑘 = 𝑚 | | | 𝑋 𝑘 = 𝑖 ) = Pr ( 𝑌 (1) 𝑘 = 𝑗 | | | 𝑋 𝑘 = 𝑖 ) Pr ( 𝑌 (2) 𝑘 = 𝑚 | | | 𝑋 𝑘 = 𝑖 ) . ( 3
)
The transition law between two successive states is modelled by the following probability rule:

𝑎 𝑖𝑗 = Pr(𝑋 𝑘 = 𝑗|𝑋 𝑘-1 = 𝑖) = Λ 𝑞 (𝑗 -𝑖) + 𝜖 𝑁 ∑ 𝑗 ′ =1 ( Λ 𝑞 (𝑗 ′ -𝑖) + 𝜖 ) ( 4
)
where 𝜖 is a small coefficient to avoid numerical problems and Λ 𝑞 (.) is the triangular function of width 𝑞:

Λ 𝑞 (𝑥) = 1 𝑞 max (𝑞 -|𝑥|, 0) . ( 5
)
The emission rules for the ECG and the PCG are assumed to be Gaussian mixture centered around the latent state: ) .

(

) 6
The first standard deviation 𝜎 (𝑒𝑐𝑔) 1 (or 𝜎 (𝑝𝑐𝑔)

1

) is small and models the fHR estimation close to the actual fHR. The second standard deviation 𝜎 (𝑒𝑐𝑔)

2

(or 𝜎 (𝑝𝑐𝑔) 2

) is large and models the error of fHR estimation. 𝜋 (𝑒𝑐𝑔) (or 𝜋 (𝑝𝑐𝑔) ) is the ratio between the two Gaussians. Finally, the coefficient 𝛼 (𝑒𝑐𝑔)

𝑖

(or 𝛼 (𝑝𝑐𝑔)

𝑖

) is a normalization coefficient to ensure the sum of all probability is equal to one.

Parameter values. Considering clinical specifications of fHR, the set of all possible states for the hidden state and the observation states are fixed to 30, 250 .

It is worth noting that in the emission laws Eq. 6, only the means of the normal distributions depend on the state value 𝑋 𝑘 . The emission parameters are estimated from fHR estimations of the volunteer pregnant women who participated to the clinical study presented in Sec. In the same way, only the mean of the transition law Eq. 4 depends on the previous state value 𝑋 𝑘-1 . Moreover, 𝑞 is the width of the triangular law which describes the variation of fHR. Considering the sampling frequency of fHR estimations at 4Hz, such a transition law discourages then the fHR variation to be larger than 𝑞 beat per minute (bpm) in 250ms. According to real data from reference CTG, it appears that variations of fHR in 250ms are lower than 3 bpm in 98% of cases, thus the value of 𝑞 is fixed to 3 bpm. Finally, the 𝜖 coefficient allows to avoid some numerical problems, it is arbitrary fixed to 10 -10 .

At last, an initial probability law Pr(𝑋 0 ) is needed to merge the two estimations. According to [START_REF] Von Steinburg | What is the "normal" fetal heart rate?[END_REF], the fHR generally varies between 110-160 bpm. Pr(𝑋 0 ) =  (135, 50) is then proposed to have a law centered around the mean fHR and whose probability Pr(𝑋 0 ∈ [110, 160]) = 40%. Because there are no temporal links between two successive 𝑋 (1) and two successive 𝑋 (2) , the model in (a) can be reduced to the model in (b) which is retained in this study.

Classical Viterbi's algorithm

The VA [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF] is an algorithm which estimates the latent path 𝑋 0∶𝑇 according to the observation paths 𝑌 (1) 0∶𝑇 and 𝑌 (2) 0∶𝑇 . The aim of the VA is to find the latent path X0∶𝑇 which maximizes the conditional probability Pr

( 𝑋 0∶𝑇 | | | 𝑌 (1) 0∶𝑇 , 𝑌 (2)

0∶𝑇

) ,

i.e. :

X0∶𝑇 = argmax 𝑋 0∶𝑇 Pr ( 𝑋 0∶𝑇 | | | 𝑌 (1) 0∶𝑇 , 𝑌 (2) 0∶𝑇 ) . ( 7
)
Because of the conditional probability rule, X0∶𝑇 can also be written as:

X0∶𝑇 = argmax 𝑋 0∶𝑇 Pr ( 𝑋 0∶𝑇 , 𝑌 (1) 0∶𝑇 , 𝑌 (2) 0∶𝑇 ) . ( 8
)
In the HMM, the state 𝑋 𝑘 is assumed to depend on the previous states. This is why, to solve the previous equation, one can use the recursive following term:

𝜇(𝑋 𝑘 ) = max 𝑋 0∶𝑘-1 Pr ( 𝑋 0∶𝑘 , 𝑌 (1) 0∶𝑘 , 𝑌 (2) 0∶𝑘 ) (9)
Because of assumptions made on the HMM, one can write 𝜇(𝑋 𝑘 ) as a temporal sequence which depends on the previous value 𝜇(𝑋 𝑘-1 ) and the new observation(s).

• If 𝑘 > 0, then 𝜇(𝑋 𝑘 ) is given by:

Limits of the classical VA

In the classical VA, the first step consists in computing the vector:

log 𝜇(𝑋 𝑘 ) = ( log 𝜇(𝑋 𝑘 = 1), … , log 𝜇(𝑋 𝑘 = 𝑁) ) 𝑇
for each 𝑘 using Eq. 9 (see Recursion in Algo 1). Then, the path X0∶𝑇 (fHR (fusion) ) is deduced with a back-tracking loop. For fHR estimation, such a selection rule can fail to provide a good estimation. If one of the modality estimations is wrong over a long time, then there is no guarantee that fHR (fusion) will stay on the fHR. This is in particular true for the ECG modality where the mHR can be estimated instead of the fHR if the mother component is not correctly removed from the ECGa signal. One illustration on real data recording is given in Fig. 6. The classical VA has been applied on 30 minutes of fHR estimations and the result of the fusion is displayed on a 50s-window for better visualization. Fig. 6a gives fHR (ecg) , fHR (pcg) and fHR (fusion) according to the VA. Estimations are compared with the reference fHR from the CTG (fHR (ref) ) and the reference mHR estimated from the ECGt (mHR (ref) ) by R-peaks detection. In this example, the algorithm to remove the mother component on ECGa failed and fHR (ecg) modality gives the mHR all the time. On the other hand, fHR (pcg) gives the expected heart rate. With such inputs, fHR (fusion) with the classical VA alternates between the fHR and the mHR.

Fig. 6b gives values of log 𝜇(𝑋 𝑘 ) which allows to estimate fHR (fusion) in Fig. 6a. The Y-axis values are the log 𝜇(𝑋 𝑘 ) and the X-axis corresponds to the states of fHR (fusion) . For each 𝑘, values of log 𝜇(𝑋 𝑘 ) decrease in time: log 𝜇(𝑋 𝑘 = 𝑖) > log 𝜇(𝑋 𝑘+1 = 1), ∀𝑖 ∈ 1, 𝑁 . It is then possible to read the time evolution of log 𝜇(𝑋 𝑘 ) in the figure be reading it from the top to the bottom. Local maxima of log 𝜇(𝑋 𝑘 ) are given by the "*" points for each 𝑘. Two main paths of local maxima can be observed in the figure: the mHRs' path (between 50-100 bpm) and the fHRs' path (between 100-150 bpm). Dashed lines are log 𝜇(𝑋 𝑘 ) synchronized with instant 𝑘 where fHR (fusion) starts and ends to follow the mHRs' path. This example shows that the VA does not ensure that the fHR (fusion) always follows the fHRs' path despite the presence of local maxima on this path.

The VA is therefore modified to overcome this limitation by changing the termination and the back-tracing steps of the classic VA. The main idea of the modified VA consists in using the location of all local maxima of log 𝜇(𝑋 𝑘 ) and the estimated state 𝑋 𝑘-1 to select the new state 𝑋 𝑘 at each 𝑘.

Algorithm 1 Classical Viterbi's algorithm for bimodal HMM Require: 𝑌 (1) 0∶𝑇 , 𝑌 (2) 0∶𝑇 and 𝜋 0 = Pr

( 𝑋 0 ) . 𝑗 1 ← 𝑌 (1) 0 ⊳ Initialization 𝑗 2 ← 𝑌 (2) 0 for 𝑖 ∈ 1, 𝑁 do log ( 𝜇(𝑋 0 = 𝑖) ) ← log 𝜋 0 + log ( 𝑏 (1) 𝑖𝑗 1 ) + log ( 𝑏 (2) 𝑖𝑗 2 ) 𝜓 0 (𝑋 0 = 𝑖) ← 0 end for for 𝑘 ∈ 1, 𝑇 do ⊳ Recursion 𝑗 1 ← 𝑌 (1) 𝑘 𝑗 2 ← 𝑌 (2) 𝑘 for 𝑖 ∈ 1, 𝑁 do log 𝜇(𝑋 𝑘 = 𝑖) ← max 𝑛 { log 𝜇(𝑋 𝑘-1 = 𝑛) + log ( 𝑎 𝑛𝑖 ) + log ( 𝑏 (1) 𝑖𝑗 1 ) + log ( 𝑏 (2) 𝑖𝑗 2 )} 𝜓 𝑘 (𝑋 𝑘 = 𝑖) ← argmax 𝑛 { log 𝜇(𝑋 𝑘-1 = 𝑛) + log ( 𝑎 𝑛𝑖 ) + log ( 𝑏 (1) 𝑖𝑗 1 ) + log ( 𝑏 (2)

Modified Viterbi's algorithm

In the proposed modified VA, the termination step and the path backtracking step (i.e. the selection rule of the state 𝑋 𝑘 given in Algo 1) are modified. The initialization step and the recursion step (see Algo 1) remain the same. The transition rule (i.e. the triangular law) of the state (see Eq. 4) discourages large variations of the state (±𝑞 bpm between two successive states). However, the presence of noise can modify the location of the global maximum of log 𝜇(𝑋 𝑘 ) with a distance to the previous state greater than 𝑞. To ensure the variation of the state remains small, the selection rule of the state 𝑋 𝑘 is modified by choosing the closest local maximum with the previous state 𝑋 𝑘-1 . If none of the local maxima has a distance with the previous state 𝑋 𝑘-1 lower than 𝑞, then the current state is rejected (𝑋 𝑘 = NaN or "not a number") and the next state 𝑋 𝑘+1 is now researched. For successive missing states, the distance to reject the closest local maximum is proportional to the time distance between the current state to evaluate and the latest known state. For example, let ∀𝑖 ∈ 𝑘 -(𝑝 + 1), 𝑘 -1 , 𝑋 𝑖 = NaN and 𝑋 𝑘-𝑝 ≠ NaN. The distance between 𝑋 𝑘 and the closest local maximum at time 𝑋 𝑘-𝑝 has to be lower than 𝑝𝑞 to be not rejected. In the case, where one of the fHR estimation input gives the mHR, the rule to ignore a state at a given time allows not to shift directly to the mHR if the fHR evolution is lost.

Algo 2 gives the modified VA and Fig. 7 is an illustration of the modified VA functioning.

Experiments

Evaluation is carried out on real data (Sec. 2.1) by comparison of fHR estimation performances between the modified VA and the classical VA.

Testing configurations

The classical VA and the modified VA are compared on six patients with the bimodal HMM (Fig. 5b). The HMM's parameters are given in the last paragraph of Sec. 3.1. For the six patients, the initial probability law is centered around the expected fHR (Pr(𝑋 0 ) =  (135, 50)). For real data, the true fHR and the true mHR are unknown. However, fHR estimation from CTG can be considered as the reference fHR and the reference mHR can be estimated thanks to the detection of mother R peaks on the thoracic ECGt signal. The quality of ECGt allows to have an accurate estimation of the mHR. Results of fusion algorithms fHR (fusion) are compared to the reference fHR and the reference mother heart rate (mHR) denoted, respectively, fHR (ref) and mHR (ref) .

Evaluation criteria

The two first evaluation criteria to compare algorithm are based on outliers estimations. The set of outliers is defined as the set of indices of fHR (fusion) values that differ from fHR (ref) by more than 12.5 bpm. The choice of 12.5 bpm is adapted to the clinical criterion for the fHR normal variability (see [START_REF] Martin | Fetal heart rate during labour: definitions and interpretation[END_REF]). The same threshold for the definition of outliers with mHR (ref) is kept to remain comparable, but without any clinical justification.

Let 𝑘 ∈ 0, 𝑇 be the discrete time vector. The fetal non-outlier's ratio 𝑜𝑢𝑡 fHR is defined by: ( * ) : Bold texts give best 𝑜𝑢𝑡 fHR between the two algorithms.

𝑜𝑢𝑡 fHR = 𝑇 ∑ 𝑘=0 𝛿 ( | | | fHR (ref) 𝑘 -fHR (fusion) 𝑘 | | | < 12.5 ) 𝑇 + 1 , ( 14
)
( * * ) : Underlined texts give best 𝑜𝑢𝑡 mHR between the two algorithms.

[1]). The signal loss definition is:

signal loss = 𝑇 ∑ 𝑘=0 𝛿 ( fHR (fusion) 𝑘 = NaN ) 𝑇 + 1 (16)
where NaN (not a number) refers to a missing value.

Results

Test results are summarized in Tab. 1, which is separated in two parts. The first one is dedicated to results on six pregnant women with the initial probability law centered around the expected fHR. And, the second part gives the result for patient n°06 with the initial probability law centered around the expected mHR. For each experiment, 𝑜𝑢𝑡 fHR , 𝑜𝑢𝑡 mHR and the signal loss are given for both the classical VA and the modified VA. Confusion with mHR When the fHR is not correctly estimated, it is interesting to compare the estimated HR with the maternal one to know if there is a confusion with this latter or if the wrong estimation of the fHR comes from another kind of noise. This is evaluated with 𝑜𝑢𝑡 mHR . The use of the modified VA decreases 𝑜𝑢𝑡 mHR for all patients but one. For the classical VA, 𝑜𝑢𝑡 mHR varies from 1% to 59%. For the modified VA, 𝑜𝑢𝑡 mHR is lower or equal to 1% for all patients except patient n°04 where both algorithms gives high 𝑜𝑢𝑡 mHR : 36% and 45% for, respectively the classical VA and the modified VA.

Accuracy of fHR estimation

Patient n°06 is the most spectacular case where fetal heart rate accuracy, 𝑜𝑢𝑡 fHR , increases from 41% with the classical VA to 96% with the modified VA. The maternal confusion, 𝑜𝑢𝑡 mHR , decreases from 59% with the classical VA to 0% with the modified VA. Additionally, to test the influence of the initial probability law, this latter is chosen as Pr(𝑋 0 ) =  (60, 50), so that the mean of this probability is now centered about the classical mHR. As one can see in the last line of Tab. 1, the modified VA now provides an heart estimation related to the mHR for 81% of the time while 𝑜𝑢𝑡 mHR = 0% with Pr(𝑋 0 ) =  (135, 50). This highlights the importance of initial state, but a rough value of fHR is sufficient for a good quality estimation, as shown previously.

Signal loss

The signal loss induced by the modified VA is between 4% and 20%. This result seems, at the first sight, disappointing compared to the classical VA for which the signal loss is always 0%. It has to be noted that a 0% signal loss does not mean necessary that the fetal heart rate estimation is accurate: it has to be analyzed with regard to the 2 other criteria 𝑜𝑢𝑡 mHR and 𝑜𝑢𝑡 fHR .

Fig. 8 illustrates results of fHR estimations from the two algorithms, for 3 subjects (Patients n°01, 05 and 06). In each figure, the two plots on the left are the input estimations from the ECG modality (top) and the PCG modality (bottom). And the two right plots are the output estimations from the classical VA (top) and the modified VA (bottom). Results are given for the initial probability rules centered around the expected fHR, Pr(𝑋 0 ) =  (135, 50). In Fig. 8a, for patient n°01, the modified VA does not increase 𝑜𝑢𝑡 fHR and does not decrease 𝑜𝑢𝑡 mHR . The signal loss (with the modified VA) is equal to 7%. This is shown with less isolated points in the plot of fHR estimation with the modified VA. Fig. 8b are results of patient n°05 where the modified VA decreases the confusion with the maternal heart rate, 𝑜𝑢𝑡 mHR , from 5% to 0%. This can be shown around 210-240s and about 520s. Finally, results of patient n°06 are given in Fig. 8c.

In this experiment, one can see that the fetal heart rate estimation from the ECG modality follows the mHR. The fetal heart rate estimation with the classical VA switches between the two heart rates (maternal and fetal ones) while the fetal heart rate estimation with the modified VA succeeds to avoid confusion with the maternal heart rate.

Discussion

The multimodal analysis, by merging ECG and PCG, is a promising way to have a robust fHR estimation from noisy observations. Our previous work on fusion of estimations using multimodal HMM [START_REF] Souriau | Fusion of estimations from two modalities using the Viterbi's algorithm: application to fetal heart rate monitoring[END_REF] showed that merging estimations from different modalities helps to get better estimation in many cases. However, if one modality is very noisy, the estimation can fail. In particular, the presence of the mother component in ECGa can make more difficult the estimation of the fHR: indeed, if the mother component is not correctly removed, the mHR might be estimated instead of the fHR. This issue can lead to a misestimation after the fusion of monomodal fHR estimations if mHR is estimated from the ECGa for several seconds. In this paper, the classical VA is modified in order to be adapted to the fHR estimation application and in particular, to avoid the confusion with the mHR. The proposition is to discourage large variations between successive estimated values.

Advantages of the proposed modified VA

The modified VA succeeds to reduce the confusion with the mHR in almost all cases with the mother non-outlier's ratio 𝑜𝑢𝑡 mHR lower or equal to 1%. This improvement is particularly important in the clinical context [START_REF] Bhogal | Maternal and fetal heart rate confusion during labour[END_REF][START_REF] Buisson | Confusion des enregistrements cardiaques foetal et maternel au cours de l'expulsion : impact sur l'état néonatal et le taux d'extraction instrumentale[END_REF] where confusion between fetal and maternal heart rates can lead to emergency cesarean sections. In addition, the accuracy of the fHR estimation, evaluated with the fetal non-outliers' ratio 𝑜𝑢𝑡 fHR , is improved for 50% of patients, with the modified VA in comparison to the classical VA. Moreover, the accuracy of the fHR estimation never decreases by more than 5%. The design of the modified VA allows to ignore some estimations when all possible estimations at a given time are aberrant. The signal loss induced by the new algorithm remains low and always acceptable for clinical uses (≤ 20% for all patients [START_REF] Ayres-De Campos | FIGO consensus guidelines on intrapartum fetal monitoring: Cardiotocography[END_REF]). On the contrary, the classical VA is never subject to signal loss: this means that an estimation is provided at each time even if the estimation is obviously wrong. The philosophy behind the design of the modified VA is to limit outliers to help the interpretation of the estimated fetal heart rate by the clinicians: it is better to provide no estimation than giving a wrong one. Indeed, according to clinicians, having some missing isolated fHR values over fHR estimation has a very few impact when reading the fHR signal, while aberrant values can make the reading more difficult.

During delivery, clinicians control the fetal heart rate in sequences of time duration of a few minutes [START_REF] Ayres-De Campos | FIGO consensus guidelines on intrapartum fetal monitoring: Cardiotocography[END_REF]: e.g. the baseline (the mean level of the most constant fHR) is estimated on segments of 10 minutes while the variability (oscillations in the fHR) is estimated on 1-minute segments. Both the classical VA and the modified VA merge 10 minutes of fHR estimations in less than 1 second on a processor Intel(R) Core(TM) i7-10610U CPU @1.80GHz-2.30GHz with Matlab 2022a. Such results show that these algorithms may be adapted for clinical uses in real time.

Limits

The fHR estimation after fusion depends on the input estimations. If input estimations are too noisy (as illustrated with patient n°04), both algorithms fail to provide a good estimation.

The classical VA provides better results than the modified VA when the initial probability rule is centered around the expected mHR. In that case, the classical VA allows the estimation to switch back to the expected heart rate, while the modified VA prevents the estimation from large variations (Fig. 7), so that the heart rate estimation can stay stack around the maternal heart rate. However, this situation can be simply avoid by choosing an initial rough probability.

In the configuration where one monomodal estimation is centered around the fHR and the other estimation is centered around mHR, if the merged estimation at previous time is centered around mHR, then the modified VA algorithm will force the next estimation to stay around the mHR.

Perspectives

From a methodological point of view, there are many perspectives to improve the fusion algorithm. Modifying the transition law and/or emission laws to be closer from reality is one of them, as well as adding a control of the prominence of each local maximum to reject local maximum likely associated with noise. Some changes have been tested but not presented in this paper, since most changes only slightly improve some criteria and degrade the other with small percentages.

As previously mentioned, the quality of input estimations has an impact on the fusion output. Changes on the preprocessing step and the monomodal fHR estimation should not be neglected and are always under development [START_REF] Lafaye De Micheaux | Residual convolutional autoencoder combined with a non-negative matrix factorization to estimate fetal heart rate[END_REF][START_REF] Reséndiz | Hidden markov model in nonnegative matrix factorization for fetal heart rate estimation using physiological priors[END_REF].

From a clinical point of view, this paper presents a proof of concept to estimate the fHR from two modalities. Tests have been performed on a small number of pregnant women and without considering the health status. Additional signal acquisitions are needed on more patients for a clinical validation of the proposed methodology.

Conclusion

In this study, two fetal heart rate estimations, obtained separately from an ECGf and a PCGf, are merged into a new estimation thanks to a modified Viterbi's algorithm. Because of the presence of the maternal heart rate in the fetal heart rate estimation, even the fusion of estimations using the classical Viterbi's algorithm can lead to confusion between the maternal and fetal heart rates. Therefore a modified Viterbi's algorithm has been proposed to limit large variations in the fetal heart rate estimation.

Results on real pregnant women data showed that the proposed modified Viterbi's algorithm succeeds to improve the fetal heart rate estimation by reducing both the presence of outliers estimation and the confusion with the maternal heart rate. Thus, fusion of estimations from two modalities seems a promising approach for more robust fetal heart rate monitoring in clinical situations.

  (a) Scheme of the heart. (b) Synchronized ECG and PCG signals.

Figure 1 :

 1 Figure 1: Heart and biosignals. (a) is a scheme which summarizes the different part of the heart and (b) gives an illustration of the theoretical ECG (top) and PCG (bottom) signals. The location of P,Q,R,S,T waves are given for the ECG signal. S 1 and S 2 in the PCG signal give the location of, respectively, systoles and diastoles.

  (a) Signals recording protocol. (b) Methodology of fHR estimation.

Figure 2 :

 2 Figure 2: Recording and fHR estimation. (a) is a scheme of the recording protocol of the studied database. The CTG (not represented here) has also been recorded for evaluation. (b)gives the methodology of fHR estimation from the ECG and the PCG modality. After pre-processing, the monomodal fHR estimation step is performed on the ECG and the PCG separately. This step similar for the two modalities is detailed in[START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF][START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF] and summarized in Sec. 2.2. The fusion step is detailed in Sec. 3.

Figure 3 :

 3 Figure 3: Monomodal fHR estimation steps. (a) and (b) give, respectively, fHR estimation for the ECG modality and the PCG modality. From top to bottom, spectrograms of the ECGf and the PCGf (Fs = 1kHz, window = 4s, shift = 250ms, zero padding ratio = 4) are first given. Then, excitation matrices 𝐻(𝑒) obtained from spectrograms using the NMF algorithm[START_REF] Dia | Fetal heart rate estimation from a single phonocardiogram signal using non-negative matrix factorization[END_REF][START_REF] Dia | Fetal heart rate estimation by non-invasive single abdominal electrocardiography in real clinical conditions[END_REF] are provided. fHR estimations deduced from excitation matrices are shown in the last row.

Figure 4 :

 4 Figure 4: A classic HMM with one latent variable 𝑋 𝑘 and one observation 𝑌 𝑘 . Links represent the conditional dependencies between variables. These dependencies are modelled with a transition matrix 𝐴 and an emission matrix 𝐵.

{ 1 ) 2 ) + ( 1 -= 1 ) 2 ) + ( 1 -

 121121 Pr(𝑌 (𝑒𝑐𝑔) 𝑘 = 𝑗|𝑋 𝑘 = 𝑖) = 𝛼 (𝑒𝑐𝑔) 𝑖 ( 𝜋 (𝑒𝑐𝑔)  (𝑗; 𝑖, (𝜎 (𝑒𝑐𝑔) 𝜋 (𝑒𝑐𝑔) ) (𝑗; 𝑖, (𝜎 (𝑗|𝑋 𝑘 = 𝑖) = 𝛼 (𝑝𝑐𝑔) 𝑖 ( 𝜋 (𝑝𝑐𝑔)  (𝑗; 𝑖, (𝜎 (𝑝𝑐𝑔) 𝜋 (𝑝𝑐𝑔) ) (𝑗; 𝑖, (𝜎 (

2 . 1 . 1 = 2 . 4 , 2 = 1 = 1 . 6 , 2 =

 2112421162 Parameter values of the emission matrix for the ECG modality are 𝜎 (𝑒𝑐𝑔) 𝜎 (𝑒𝑐𝑔) 56.2 and 𝜋 (𝑒𝑐𝑔) = 0.6. Parameter values of the emission matrix for the PCG modality are 𝜎 (𝑝𝑐𝑔) 𝜎 (𝑝𝑐𝑔) 63.1 and 𝜋 (𝑝𝑐𝑔) = 0.5.

  (a) Complete model. (b) Model for this study.

Figure 5 :

 5 Figure 5: HMMs to model fHR's evolution. The model given in (a) is closer to the reality where the true fHR, the fHRgiven by the different modalities and the estimations are all represented. Because there are no temporal links between two successive 𝑋(1) and two successive 𝑋(2) , the model in (a) can be reduced to the model in (b) which is retained in this study.

  X𝑇 = argmax 𝑖 log 𝜇(𝑋 𝑇 = 𝑖). ⊳ Termination for 𝑘 ∈ 0, 𝑇 -1 do ⊳ Path backtracking (𝑘 from 𝑇 -1 to 0) X𝑘 = 𝜓 𝑘+1 ( X𝑘+1 ) end for Return: X0∶𝑇 .

  (a) Heart rate estimation, real example. (b) Vector values (log 𝜇(𝑋 𝑘 )) evolution.

Figure 6 :

 6 Figure 6: Limit of fHR estimation using the classical VA illustrated on 50s. (a) gives fHR (ecg) on the top, fHR (pcg) at the middle and fHR (fusion) at the bottom. fHR (ref) and mHR (ref) are, respectively, the reference fHR given by the CTG and the reference mHR estimated from the ECGt. (b) gives the vector values log 𝜇(𝑋 𝑘 ) which allows to get fHR (fusion) in (a). Each line corresponds to log 𝜇(𝑋 𝑘 ) values for all states at a given date 𝑘. Because the values of log 𝜇(𝑋 𝑘 ) decrease in time, one can read the evolution of log 𝜇(𝑋 𝑘 ) in time by shifting from the top to the down of the figure. For each instant 𝑘, local maxima are displayed on the figure. Dashed lines in (a) and (b) are synchronized and give log 𝜇(𝑋 𝑘 ) with 𝑘 corresponding to the location when fHR (fusion) starts and ends to give the mHR and not the fHR.

  Modified Viterbi's algorithm for fHR monitoring (a) Patient n°01. (b) Patient n°05. (c) Patient n°06.

Figure 8 :

 8 Figure 8: Result examples on patients n°01, n°05 and n°06. On each figure, the left plots are the input estimations from the modalities separately. The top right plot gives the result of fusion estimation using the classical VA and the bottom right figure gives the result with the modified VA. For each example, the initial probability rule Pr(𝑋 0 ) =  (135, 50) is centered around the expected fHR.

Table 1

 1 Non-outliers ratio using classical VA and modified VA.

	Pr(𝑋 0 )	Patient number	𝑜𝑢𝑡 fHR	Classical VA 𝑜𝑢𝑡 mHR signal loss 𝑜𝑢𝑡 fHR	Modified VA 𝑜𝑢𝑡 mHR signal loss
		01	92% ( * )	1%	0%	87%	1%	7%
		02	75%	8%	0%	84%	1% ( * * )	13%
	 (135, 50)	03 04	90% 38%	2% 36%	0% 0%	93% 23%	1% 45%	6% 20%
		05	92%	5%	0%	90%	0%	6%
		06	41%	59%	0%	96%	0%	4%
	 (60, 50)	06	41%	59%	0%	0%	81%	6%

  It is provided by the criterion 𝑜𝑢𝑡

	𝑜𝑢𝑡 fHR are larger for the modified VA in 50% of cases, resulting to a better estimation of the fHR with the modified VA
	than with the classical VA. The use of the modified VA never decreases 𝑜𝑢𝑡	fHR by more than 5%. The only exception is
	with patient n°04 where both algorithms fail to provide a high 𝑜𝑢𝑡 fHR : 38% and 23% for, respectively the classical VA
	and the modified VA.

fHR

. In tests where Pr(𝑋 0 ) =  (135, 50), values of

Acknowledgement

This work is supported by the French National Research Agency, as part of the SurFAO project (ANR-17-CE19-0012).

The authors thank the Clinical Investigation Centre for Innovative Technology of Grenoble, for their support to data acquisition.

𝜇(𝑋

Pr

( 𝑋 0∶𝑘-1 , 𝑋 𝑘 , 𝑌 (1) 0∶𝑘-1 , 𝑌 (1) 𝑘 , 𝑌 (2) 0∶𝑘-1 , 𝑌 (2) 𝑘 )

By using conditional probability rule and assumptions of the HMM, one can write from Eq. 10:

{ Pr ( 𝑋 0∶𝑘-1 , 𝑌 (1) 0∶𝑘-1 , 𝑌 (2) 0∶𝑘-1

)

The passage from the first line to the second one in the above equation is performed using the following property: If ∀𝑎, 𝑓 (𝑎) ≥ 0 and ∀(𝑎, 𝑏), 𝑔(𝑎, 𝑏) ≥ 0, then:

• If 𝑘 = 0, then 𝜇(𝑋 0 ) is given by:

In practice, logarithms of Eqs. 11 & 12 are used to avoid numerical issues in algorithms.

Once 𝜇(𝑋 𝑘 ) is known, then the last state 𝑋 𝑇 can be estimated by finding the state which maximizes 𝜇(𝑋 𝑇 ):

And, all the previous states can be deduced using back tracking and maximum index of each 𝜇(𝑋 𝑘 ). Algo 1 summarizes the classical VA.

Algorithm 2 Modified Viterbi's algorithm for bimodal HMM Require: 𝑞, 𝑌 (1) 0∶𝑇 , 𝑌 (2) 0∶𝑇 and 𝜋 0 = Pr where * 𝑘 refers to the value of * at time 𝑘 and 𝛿(𝑥) = 1 if 𝑥 is true, 0 else. And, the mother non-outlier's ratio 𝑜𝑢𝑡 mHR is equal to:

It is worth noting that these criteria values are between 0 and 1: the higher 𝑜𝑢𝑡 fHR , the better, while the lower 𝑜𝑢𝑡 mHR , the better. In addition to the outlier's ratio, the signal loss ratio is also computed. For the classical VA, there is no removed data, but for the modified VA an estimation at a given time can be rejected (if Algo 2). The signal loss evaluates the quantity of removed data. For clinical uses, the signal loss should be the lower as possible (≤ 20% to be acceptable