A Convolutional Neural Network-based Automatic Identification and Intervention Model for Health Surveillance Data during Postpartum Recovery Periods

1)

Generation and Classification of Pulse Waves

Pulse wave is also known as arterial beat. In the human blood circulation system, when the left ventricle diastole, it will squeeze the blood to shoot into the aorta, the blood vessels will block the blood flow directly into the veins, at this time, because many blood gathered to the proximal aorta makes the aortic pressure increase, which makes its diameter dilation, this phenomenon of pulsation dilation in the human body in some places can be felt, such as the wrist, neck, forefinger, ring finger, etc., which is often referred to as the pulse [37]. The formation of a pulse is based on two conditions: the contractile function of the heart and the elasticity of its walls. When the ventricles of the heart open, the pressure in the arteries increases and when the ventricles of the heart close, the pressure in the arteries decreases, which causes the arterial beat to propagate in the aorta and its various branches in the form of a wave, which is known as a pulse wave.

Figure 1.

The production of PPG

The oxygen saturation detection used in this paper is realized based on the photoelectric volumetric tracing technique. Whether arterial blood contains enough oxygen or not plays a vital role in maintaining human life. The oxygen content of human blood is generally judged by detecting the blood oxygen saturation. Arterial oxygen saturation (SaO2)$$\left( {Sa{O_2}} \right)$$ is the percentage of the actual amount of oxygen bound in the arterial blood to the maximum oxygen capacity of the blood that can bind oxygen. As a result of human respiration, the vast majority of the oxygen in the blood is bound to deoxyhemoglobin (Hb)$$\left( {Hb} \right)$$ and becomes oxyhemoglobin (HbO2)$$\left( {Hb{O_2}} \right)$$, with only about 1-2% of the oxygen being dissolved in other components of the blood. Therefore, the following formula can be used to approximate blood oxygen saturation: (1)SaO2=(CHbO3CHbO2+CHb)×100%$$Sa{O_2} = \left( {\frac{{{C_{Hb{O_3}}}}}{{{C_{Hb{O_2}}} + {C_{Hb}}}}} \right) \times 100\%$$

Where, C_HbO2 is the concentration of oxygenated hemoglobin in human blood tissue and C_Hb is the concentration of deoxyhemoglobin. According to Lambert-Beer law, there is: (2)I=I0×F×e−(E1C1L+E2C2L)=I0′e−(E1C1+E2C2)L$$I = {I_0} \times F \times {e^{ - \left( {{E_1}{C_1}L + {E_2}{C_2}L} \right)}} = {I_0}^{\prime} {e^{ - \left( {{E_1}{C_1} + {E_2}{C_2}} \right)L}}$$

Where I is the light intensity reflected back from the human body, I₀ is the outgoing light intensity, F is the light absorption coefficient of human tissue, E₁ and E₂ represent the absorption coefficients of oxygenated hemoglobin HbO₂ and deoxyhemoglobin Hb in the human blood tissue, C₁ and C₂ represent the concentrations of oxygenated hemoglobin HbO₂ and deoxyhemoglobin Hb in the human blood tissue, and L is the length of the reflective light path. As the heart beats, the human blood tissue fluctuates and changes, the reflected light path changes ΔL and the corresponding reflected light intensity changes (IAC+IDC)$$\left( {{I_{AC}} + {I_{DC}}} \right)$$: (3)IAC+IDC=IDCe−(E1C1+E2C2)ΔL$${I_{AC}} + {I_{DC}} = {I_{DC}}{e^{ - \left( {{E_1}{C_1} + {E_2}{C_2}} \right)\Delta L}}$$

Taking logarithms on both sides of equation (3) gives: (4)ln[(IAC+IDC)/IDC]=−(E1C1+E2C2)ΔL$$\ln \left[ {\left( {{I_{AC}} + {I_{DC}}} \right)/{I_{DC}}} \right] = - \left( {{E_1}{C_1} + {E_2}{C_2}} \right)\Delta L$$

The percentage of AC and DC flow in reflected light is very small, therefore: (5)ln[(IAC+IDC)/IDC]≈IAC/IDC$$\ln \left[ {\left( {{I_{AC}} + {I_{DC}}} \right)/{I_{DC}}} \right] \approx {I_{AC}}/{I_{DC}}$$

Combining equations (4) and (5) gives: (6)IAC/IDC=−(E1C1+E2C2)ΔL$${I_{AC}}/{I_{DC}} = - \left( {{E_1}{C_1} + {E_2}{C_2}} \right)\Delta L$$

The amount of reflected light path change ΔL is an unknown quantity, and choosing two different wavelengths of incidence eliminates the unknown quantity ΔL. Assuming that the two wavelengths of light are λ₁ and λ₂, respectively: (7)IACA1λ1/IDCλ1IACλ2/IDCλ2=E1λ1C1+E2λ1C2E1λ2C1+E2λ2C2$$\frac{{{I_{AC}}A_1^{{\lambda _1}}/{I_{DC}}^{{\lambda _1}}}}{{I_{AC}^{{\lambda _2}}/{I_{DC}}^{{\lambda _2}}}} = \frac{{E_1^{{\lambda _1}}{C_1} + E_2^{{\lambda _1}}{C_2}}}{{E_1^{{\lambda _2}}{C_1} + {E_2}^{{\lambda _2}}{C_2}}}$$

Combining equation (7) with the original defining formula for arterial oxygen saturation gives: (8)SaO2=(E2λ2E1λ1−E2λ1×IACλ1/IDCλ1IAC/IDC−E2λ1E1λ2−E2λ2)×100%$$Sa{O_2} = \left( {\frac{{E_2^{{\lambda _2}}}}{{E_1^{{\lambda _1}} - E_2^{{\lambda _1}}}} \times \frac{{I_{AC}^{{\lambda _1}}/I_{DC}^{{\lambda _1}}}}{{{I_{AC}}/{I_{DC}}}} - \frac{{E_2^{{\lambda _1}}}}{{E_1^{{\lambda _2}} - E_2^{{\lambda _2}}}}} \right) \times 100\%$$

E1λ1$$E_1^{{\lambda _1}}$$, E2λ1$$E_2^{{\lambda _1}}$$, E1λ2$$E_1^{{\lambda _2}}$$, and E2λ2$$E_2^{{\lambda _2}}$$ in equation (8) are constants. Therefore, let: A=E2λ2/(E1λ1−E2λ1)$$A = E_2^{{\lambda _2}}/\left( {E_1^{{\lambda _1}} - E_2^{{\lambda _1}}} \right)$$, B=E2λ1/(E1λ2−E2λ2)$$B = E_2^{{\lambda _1}}/\left( {E_1^{{\lambda _2}} - E_2^{{\lambda _2}}} \right)$$, R=(IACλ1/IDCλ1)/(IACλ2/IDCλ2)$$R = \left( {IA{C^{{\lambda _1}}}/ID{C^{{\lambda _1}}}} \right)/\left( {IA{C^{{\lambda _2}}}/ID{C^{{\lambda _2}}}} \right)$$, carry over to the oxygen saturation formula: (9)SaO2=(A×R−B)×100%$$Sa{O_2} = \left( {A \times R - B} \right) \times 100\%$$

1)

Theoretical basis of blood pressure

Due to the heart’s contraction and contraction, blood flows in the circulatory organs of the human body, and in the process of flow, the pressure generated by the blood on the walls of the blood vessels is the blood pressure. Generally, human blood pressure refers to arterial blood pressure, and human physiological health has a close connection, for the diagnosis of the patient’s disease, analysis and treatment, etc. provides an important basis.

On the basis of the wave conduction velocity formula of the ideal fluid elastic tube, a large number of experiments have been carried out, and a new wave velocity formula has been introduced according to the experimental results: (10)C=KEHρl$$C = K\sqrt {\frac{{EH}}{{\rho l}}}$$

In equation (10): C is the pulse wave velocity. K is called Moons constant and is dimensionless. E is the modulus of elasticity of the ideal fluid elastic tube. H is the thickness of the elastic tube. ρ is the fluid density inside the elastic tube. d is the inner diameter of the ideal fluid elastic tube.

Later, an expression for the relationship between the modulus of elasticity of blood and the pressure at the wall of the blood tube was introduced on its basis: (11)E=E0⋅ep⋅γ$$E = {E_0} \cdot {e^{p \cdot \gamma }}$$

In Eq. (11), E₀ is the modulus of elasticity when the pressure within the blood tube wall reaches equilibrium. γ is a coefficient characterizing the blood tube wall. P is the blood pressure value.

According to the physical velocity-distance relationship expression, there is: (12)C=DT$$C = \frac{D}{T}$$

Where D is the conduction distance of the pulse wave. T is the pulse wave conduction time.

According to Eq. (10) to Eq. (12), it can be introduced: (13)P=1γ[ln(ρlD2K2E0H)−2lnT]$$P = \frac{1}{\gamma }\left[ {\ln \left( {\frac{{\rho l{D^2}}}{{{K^2}{E_0}H}}} \right) - 2\ln T} \right]$$

Under the preconditions of constant vascular elasticity as well as constant pulse wave conduction distance, the first term in equation (13) is a constant; therefore, differentiating both sides of equation (13) yields the relationship between blood pressure change and pulse wave conduction time: (14)ΔP=−2γTΔT$$\Delta P = - \frac{2}{{\gamma T}}\Delta T$$

Where: ΔP is the amount of change in P and ΔT is the amount of change in T. 2)

Methods of acquiring pulse wave conduction time PTT

There are two commonly used methods to obtain pulse wave conduction time (PTT): the first method is to select two points on the same pulse wave conduction pathway, the blood flow drives the propagation of the pulse wave, and by detecting the phase difference between the pulse waveforms of these two points, the conduction time of the pulse wave between these two points can be obtained. Another method is to collect the ECG signal and pulse wave signal at the same time, in a certain cardiac cycle, take the ECG signal as the starting point, and the characteristic points such as the peak or trough of the pulse wave as the end point to calculate the time difference, which can be defined as the pulse wave conduction time PTT.

Each characteristic point of the PPG signal: main wave, descending middle isthmus, heavy beat wave, etc. are located at the extreme value point in the waveform graph. As mentioned above, for the finding of troughs and peaks, the AMPD algorithm is used in this paper, the basic principle of which is shown as follows: let X=[x1,x2,…,xi,…,xN]$$X = \left[ {{x_1},{x_2}, \ldots ,{x_i}, \ldots ,{x_N}} \right]$$ be a univariate uniformly sampled signal (with baseline drift removed) [38].

First the local maximum scale (LMS) needs to be calculated. This is done by using a moving window of different sizes w_k = 2k where k = 1, 2, …, [N/2] − 1, different k are called different scales. For each scale k and i = k + 2, …, N − k + 1, there are: (15)mk,i={ 0, xi−1>xi−k−1∧xi−1>xi+k−1 r+α, otherwise$${m_{k,i}} = \left\{ {\begin{array}{*{20}{c}} {0,}&{{x_{i - 1}} > {x_{i - k - 1}} \wedge {x_{i - 1}} > {x_{i + k - 1}}} \\ {r + \alpha ,}&{otherwise} \end{array}} \right.$$

Where r is a uniformly distributed random number in the range [0, 1] and α is a constant factor (α=1)$$\left( {\alpha = 1} \right)$$.

For i = 1, ⋯, k + 1 and i = N − k + 2, …, N, r + α is assigned to m_k,i. The result is shown in matrix M. (16)M=[ m1,1 m1,2 … m1,N m2,1 m2,2 … m2,N ⋮ ⋮ ⋱ ⋮ mL,1 mL,2 … mL,N]=(mk,j)$$M = \left[ {\begin{array}{*{20}{c}} {{m_{1,1}}}&{{m_{1,2}}}& \ldots &{{m_{1,N}}} \\ {{m_{2,1}}}&{{m_{2,2}}}& \ldots &{{m_{2,N}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{m_{L,1}}}&{{m_{L,2}}}& \ldots &{{m_{L,N}}} \end{array}} \right] = \left( {{m_{k,j}}} \right)$$

where row k contains the value of window length w_k.

Therefore, all elements of matrix M of shape L × N are in the range of [0,1+α]$$\left[ {0,1 + \alpha } \right]$$. This matrix M is called the LMS of the signal X.

This is followed by a row-by-row summation of the LMS matrix M. (17)γk=∑i=1Nmk,i,fork∈{1,2,…,L}$${\gamma _k} = \sum\limits_{i = 1}^N {{m_{k,i}}} ,{\text{ }}for{\text{ }}k \in \left\{ {1,2, \ldots ,L} \right\}$$

Vector γ=[γ1,γ2,…,γi,…,γL]$$\gamma = \left[ {{\gamma _1},{\gamma _2}, \ldots ,{\gamma _i}, \ldots ,{\gamma _L}} \right]$$ contains scale-dependent information about the distribution of zeros (i.e., extreme points). Let λ be the local maximum of the current scale, which is used to reshape the LMS matrix M in such a way that when k > λ, the corresponding m_k,i is removed, resulting in a new λ × N matrix Mr=(mk,i)$${M_r} = \left( {{m_{k,i}}} \right)$$, where i ∈ 1, 2, …, N, k ∈ 1, 2, …, λ.

Finally, the detection of the waveforms is performed in the following way. (18)σi=1λ−1∑k=1λ[(mk,i−1λ∑k=1λmk,i)2]12,fori∈1,2,…,N$${\sigma _i} = \frac{1}{{\lambda - 1}}\sum\limits_{k = 1}^\lambda {{{\left[ {{{\left( {{m_{k,i}} - \frac{1}{\lambda }\sum\limits_{k = 1}^\lambda {{m_{k,i}}} } \right)}^2}} \right]}^{\frac{1}{2}}}} ,{\text{ }}for{\text{ }}i \in 1,2, \ldots ,N$$

By indexing i at record σ_i = 0 and storing the index as a vector P=[p1,p2,…,pq,…,pN]$$P = \left[ {{p_1},{p_2}, \ldots ,{p_q}, \ldots ,{p_N}} \right]$$, then N^$$\hat N$$ is the number of signal X crests and P is the crest index value.

In particular, in order to deal with the elimination of the influence of different scales because of the individual subject’s own differences, the Z-Score normalization is applied to process the signal into data with 0 as the mean and 1 as the variance, in order to reduce the influence of the scale, the characteristics, the difference in distribution, etc., on the feature extraction. The formula is shown below. (19)Zi=xi−μσ$${Z_i} = \frac{{{x_i} - \mu }}{\sigma }$$

Based on the discrete Fourier transform (DFT), its power spectral density (PSD) is further calculated, which can quantify the variations in the PPI sequence, provide basic information on the distribution of energy changes with frequency, and quantitatively assess the modulation of sympathetic and vagal nerves.

The PSD calculation method is derived as follows.

Let f(t) be a PPI sequence, and intercept a segment of |t|≤T2$$\left| t \right| \leq \frac{T}{2}$$ from f(t) to obtain a truncated function f_T(t): (20)fT(t)={ f(t) (t|≤T2) 0 (t|>T2)$${f_T}(t) = \left\{ {\begin{array}{*{20}{c}} {f(t)}&{\left( {t\left| { \leq \frac{T}{2}} \right.} \right)} \\ 0&{\left( {t\left| { > \frac{T}{2}} \right.} \right)} \end{array}} \right.$$

If T is a finite value, then f_T(t) energy is also finite.

Let F[fT(t)]=FT(ω)$$F\left[ {{f_T}(t)} \right] = {F_T}(\omega )$$, the energy E_T of f_T(t) at this point can be expressed as: (21)ET=∫−∞∞fT2(t)dt=12π∫−∞∞|FT(ω)|2d$${E_T} = \int_{ - \infty }^\infty {f_T^2} (t)dt = \frac{1}{{2\pi }}\int_{ - \infty }^\infty {{{\left| {{F_T}(\omega )} \right|}^2}} d$$

By: (22)∫−∞∞fT2(t)dt=∫−T/2T/2f2(t)dt$$\int_{ - \infty }^\infty {f_T^2} (t)dt = \int_{ - T/2}^{T/2} {{f^2}} (t)dt$$

The average power of f(t) is obtained as: (23)P=limT→∞∫−T2T21Tf2(t)dt=12π∫−∞+∞limT→∞1T|FT(ω)|2dω$$P = \mathop {\lim }\limits_{T \to \infty } \int_{ - \frac{T}{2}}^{\frac{T}{2}} {\frac{1}{T}} {f^2}(t)dt = \frac{1}{{2\pi }}\int_{ - \infty }^{ + \infty } {\mathop {\lim }\limits_{T \to \infty } } \frac{1}{T}{\left| {{F_T}(\omega )} \right|^2}d\omega$$

As T is increasing, so is the energy. f_T(t) → f(t) when T → ∞, at which point a limit may exist for 1T|FT(ω)|2$$\frac{1}{T}{\left| {{F_T}(\omega )} \right|^2}$$. Assuming that the limit exists, define it as a function of the power spectral density of f(t), denoted as P(ω). i.e., the power spectrum of f(t) is: (24)P(ω)=liml→∞|FT(ω)|2T(rad2/hz)$$P(\omega ) = \mathop {\lim }\limits_{l \to \infty } \frac{{{{\left| {{F_T}(\omega )} \right|}^2}}}{T}\left( {ra{d^2}/hz} \right)$$

Therefore the average power of f(t) is: (25)P=12π∫−∞∞P(ω)dω$$P = \frac{1}{{2\pi }}\int_{ - \infty }^\infty P (\omega )d\omega$$

The power spectrum can be obtained as a reflection of the variation of signal power with frequency per unit frequency band, that is, the distribution of signal power in the frequency domain. The area of P(ω) is the total power of that signal, known as the bilateral power spectrum, and is thus taken as S(ω) = 2P(ω).

Notice the autocorrelation function of f(t). (26)R(τ)=limT→∞1T∫−T/2T/2f(t)f*(t−τ)dt$$R(\tau ) = \mathop {\lim }\limits_{T \to \infty } \frac{1}{T}\int_{ - T/2}^{T/2} f (t){f^*}(t - \tau )dt$$

According to the relevant theorem, it can be obtained: (27){ R(τ)=12π∫−∞∞ψ(ω)eiωτdω ψ(ω)=∫−∞∞R(τ)e−iωτdτ$$\left\{ {\begin{array}{*{20}{c}} {R(\tau ) = \frac{1}{{2\pi }}\int_{ - \infty }^\infty \psi (\omega ){e^{i\omega \tau }}d\omega } \\ {\psi (\omega ) = \int_{ - \infty }^\infty R (\tau ){e^{ - i\omega \tau }}d\tau } \end{array}} \right.$$

Thus the power spectral function of a power limited signal is a pair of Fourier transform with its autocorrelation function, i.e., the Wiener-Hinchin theorem. In practical experiments, the signals are obtained by sampling and behave in the form of discrete points, so the choice is the DFT, of the form: (28)F(k)=∑n=0N−1f(n)e−j2πnk/N$$F(k) = \sum\limits_{n = 0}^{N - 1} f (n){e^{ - j2\pi nk/N}}$$

The corresponding power spectral density is: (29)P(k)=limN→∞|F(k)|2N(rad2/hz)$$P(k) = \mathop {\lim }\limits_{N \to \infty } \frac{{{{\left| {F(k)} \right|}^2}}}{N}\left( {ra{d^2}/hz} \right)$$

Sample entropy is a method used to measure the complexity of a time series. It can measure the complexity of a time series by measuring the probability of new pattern generation in the signal; the greater the probability of new pattern generation, the greater the complexity of the sequence [39]. The sample entropy can be represented by SampEn(m,r,N)$$SampEn\left( {m,r,N} \right)$$, where m is the number of dimensions, which can be m or m + 1, r is the similarity tolerance, and N is the length of the time series.

In general, the sample entropy is calculated as follows for a time series x(n) = x(1), x(2), x(N) consisting of N data.

Form a sequence of vectors of dimension m by ordinal number, X_m(1), …, X_m(N − m + 1) where X_m(i) = x(i), x(i + 1), …, x(i + m − 1). m consecutive x values starting at point i.

Define the distance d[Xm(i),Xm(j)]$$d\left[ {{X_m}(i),{X_m}(j)} \right]$$ between vectors X_m(i) and X_m(j) as the absolute value of the largest difference in the elements corresponding to both. I.e: (30)d[Xm(i),Xm(j)]=maxk=0,⋯,m−1(x(i+k)−x(j+k))$$d\left[ {{X_m}(i),{X_m}(j)} \right] = \mathop {\max }\limits_{k = 0, \cdots ,m - 1} \left( {x\left( {i + k} \right) - x\left( {j + k} \right)} \right)$$

For a given X_m(i), count the number of j(1≤j≤N−m)$$j\left( {1 \leq j \leq N - m} \right)$$ for which the distance between X_m(i) and X_m(j) is less than or equal to r and denote them as B_i. For 1 ≤ i ≤ N − m, define: (31)Bim(r)=1N−m−1Bi$$B_i^m(r) = \frac{1}{{N - m - 1}}{B_i}$$

Definition Bim(r)$$B_i^m(r)$$ for: (32)B(m)(r)=1N−m∑i=1N−mBim(r)$${B^{(m)}}(r) = \frac{1}{{N - m}}\sum\limits_{i = 1}^{N - m} {B_i^m} (r)$$

Increase the dimension to m + 1 and count the number of distances between X_m+1(i) and Xm+1(j)(1≤j≤N−m,j≠i)$${X_{m + 1}}\left( j \right)\left( {1 \leq j \leq N - m,j \ne i} \right)$$ that are less than or equal to r, denoted A_i. Aim(r)$$A_i^m(r)$$ is defined as: (33)Aim(r)=1N−m−1Ai$$A_i^m(r) = \frac{1}{{N - m - 1}}{A_i}$$

Definition A^m(r) for: (34)Am(r)=1N−m∑i=1N−mAim(r)$${A^m}(r) = \frac{1}{{N - m}}\sum\limits_{i = 1}^{N - m} {A_i^m} (r)$$

Thus, B^m(r) is the probability that two sequences match m points under the similarity tolerance r, while A^m(r) is the probability that two sequences match m + 1 points. The sample entropy is defined as: (35)SampEn(m,r)=limN→∞{−ln[An(r)Bm(r)]}$$SampEn(m,r) = \mathop {\lim }\limits_{N \to \infty } \left\{ { - \ln \left[ {\frac{{{A^n}(r)}}{{{B^m}(r)}}} \right]} \right\}$$

When N is a finite value, it can be estimated using the following equation: (36)SampEn(m,r)=−ln[An(r)Bm(r)]$$SampEn(m,r) = - \ln \left[ {\frac{{{A^n}(r)}}{{{B^m}(r)}}} \right]$$

1)

LSTM algorithm

As an improvement of RNN, LSTM is used to store the state at moment t by adding a new state in the hidden layer. Let X={x(1),x(2),⋯,x(n)}$$X = \left\{ {{x^{(1)}},{x^{(2)}}, \cdots ,{x^{(n)}}} \right\}$$ be a time series of length n, x(t)$${x^{\left( t \right)}}$$ denote the sequence captured at moment t, and x(t)∈Rm$${x^{\left( t \right)}} \in {R^m}$$, be a m-dimensional vector {x1(t),x2(t),⋯,xm(1)}$$\left\{ {x_1^{(t)},x_2^{(t)}, \cdots ,x_m^{(1)}} \right\}$$ meaning the m-dimensional measurements of the tth timestamp in the time series. Where the input of the cell at the tth moment is in addition to the observation x(t)$${x^{\left( t \right)}}$$ at the tth moment, there is also the cell state c_t−1 and the output h_t−1 of the previous moment, i.e., the t − 1th moment, through the phase present the biggest difference between RNN and LSTM is that a special structure that can record and update the state of the cell is added to the LSTM cell c_t, and the function of c_t is to save the state accumulated in the history, and then to participate in the computation of the output. The LSTM does this specifically by setting up three gates: a forgetting gate, an updating gate, and an output gate to fulfill the function. The forgetting gate is responsible for managing what information needs to be forgotten by taking the inputs h_t−1 and x(t)$${x^{\left( t \right)}}$$ using the Sigmoid function to xt(t)$$x_t^{\left( t \right)}$$. Next is the update gate which controls the inputs with new information to update the values. Finally, it is the calculation of the forgotten data and the memorized data to get the new cell state C_t, which is calculated as shown in Eqs. (37) to (40). (37)ft=σ(Wf⋅[ht−1,xt]+bf)$${f_t} = \sigma \left( {{W_f} \cdot \left[ {{h_{t - 1}},{x_t}} \right] + {b_f}} \right)$$ (38)it=σ(Wi⋅[ht−1,xt]+bi)$${i_t} = \sigma \left( {{W_i} \cdot \left[ {{h_{t - 1}},{x_t}} \right] + {b_i}} \right)$$ (39)C˜i=tanh(Wc⋅[ht−1,xt]+bc)$${\widetilde C_i} = \tanh \left( {{W_c} \cdot \left[ {{h_{t - 1}},{x_t}} \right] + {b_c}} \right)$$ (40)Ct=ft×Ct−1+it×C˜t$${C_t} = {f_t} \times {C_{t - 1}} + {i_t} \times {\widetilde C_t}$$

The final step is to calculate the output value h_t at moment t using the output gate by first calculating what needs to be retained by the sigmoid cell o_t, then using the Tanh function to regularize the cell state C_t, to between -1 and 1, and finally multiplying the results of o_t and C_t, as the outputs h_t. Eqs. (41) and (42) are the calculations for h_t and o_t. (41)ot=σ(Wo⋅[ht−1,xt]+bo)$${o_t} = \sigma \left( {{W_o} \cdot \left[ {{h_{t - 1}},{x_t}} \right] + {b_o}} \right)$$ (42)ht=ot×tanh(Ct)$${h_t} = {o_t} \times \tanh \left( {{C_t}} \right)$$ 2)

Attention mechanism for processing time series

Self-attention mechanism as a kind of attention mechanism mainly fuses the relationship between input sequences through queries, keys and values and outputs them.

a={a1,a2,⋯,an}$$a = \left\{ {{a^1},{a^2}, \cdots ,{a^n}} \right\}$$ is an input sequence of length n, and b={b1,b2,⋯,bn}$$b = \left\{ {{b^1},{b^2}, \cdots ,{b^n}} \right\}$$ is an output sequence of length n with weights assigned by the attention mechanism. The following is an example of b¹, describing how the attention mechanism generates the correlation vectors α between the variables based on the individual input vectors. The input vectors are multiplied by different matrices, e.g., a¹ by W^q, a² by W^k, to obtain q and k, and finally q and k are multiplied by a dot product to obtain the correlation vector α, and the correlation vectors of each vector and a¹ are obtained α and then processed by the softmax function to obtain α′, which is computed as shown in Eq. (43). (43)a1,i′=exp(α1,i)∑jexp(α1,j)$${a^{\prime}_{1,i}} = \frac{{\exp \left( {{\alpha _{1,i}}} \right)}}{{\sum\limits_j {\exp } \left( {{\alpha _{1,j}}} \right)}}$$

After obtaining the correlation vector a1,i′$${a'_{1,i}}$$ between the first vector and the ist vector, the important information vⁱ of each input vector is extracted from the input vectors, and then vⁱ multiplied with the corresponding correlation vectors respectively to obtain the output b¹, which is calculated as shown in Eqs. (44) to (47). (44)qi=Wqai$${q^i} = {W^q}{a^i}$$ (45)ki=Wkai$${k^i} = {W^k}{a^i}$$ (46)vi=Wvai$${v^i} = {W^v}{a^i}$$ (47)b′=∑iα1,i′vi$$b^{\prime} = \sum\limits_i {{{\alpha}^{\prime}_{1,i}}} {v^i}$$ (48)A=( α1,1,α2,1,α3,1,α4,1 α1,2,α2,2,α3,2,α4,2 α1,3,α2,3,α3,3,α4,3 α1,4,α2,4,α3,4,α4,4)=( k1 k2 k3 k4)(q1,q2,q3,q4)$$A = \left( {\begin{array}{*{20}{l}} {{\alpha _{1,1}},{\alpha _{2,1}},{\alpha _{3,1}},{\alpha _{4,1}}} \\ {{\alpha _{1,2}},{\alpha _{2,2}},{\alpha _{3,2}},{\alpha _{4,2}}} \\ {{\alpha _{1,3}},{\alpha _{2,3}},{\alpha _{3,3}},{\alpha _{4,3}}} \\ {{\alpha _{1,4}},{\alpha _{2,4}},{\alpha _{3,4}},{\alpha _{4,4}}} \end{array}} \right) = \left( {\begin{array}{*{20}{l}} {{k^1}} \\ {{k^2}} \\ {{k^3}} \\ {{k^4}} \end{array}} \right)\left( {{q^1},{q^2},{q^3},{q^4}} \right)$$

Similarly, all the relation vectors are written in the form of matrix A as shown in equation (48). Normalize A to obtain A′, also known as the attention score Let the output matrix be represented using O and the input vectors be represented using V, then O is computed as shown in Equation (49).The above computational procedure shows that the output corresponding to each input is computed in parallel without causing any increase in time complexity. (49)[ b1 b2 ⋯ bn]=[ v1 v2 ⋯ vn]⋅A′ ⇒O=V⋅A′$$\begin{array}{rcl} \left[ {\begin{array}{*{20}{l}} {{b^1}}&{{b^2}}& \cdots &{{b^n}} \end{array}} \right] & = & \left[ {\begin{array}{*{20}{l}} {{v^1}}&{{v^2}}& \cdots &{{v^n}} \end{array}} \right] \cdot A^{\prime} \\ & \Rightarrow & O = V \cdot A^{\prime} \\ \end{array}$$ 3)

OC-SVM for anomaly detection

OC-SVM uses unsupervised learning to recognize sparsely populated data regions and distinguish data belonging to abnormal regions. Suppose there is a training sample x=(x1,x2,⋯,xt)$$x = \left( {{x_1},{x_2}, \cdots ,{x_t}} \right)$$ of length l. Training sample x belongs to a known class which is generally the normal class. Let ϕ be a kernel function that can project the training data into a high-dimensional space. Then the goal is to solve the objective function of the support vector machine as shown in Eq. (50) such that the data set is separated from the origin. (50)min{12‖w‖2+1vl∑i=1lεi−ρ} s.t:wϕ(xi)≥ρ−εi,i=1,2,3,⋯,l,εi≥0$$\begin{array}{c} \min \left\{ {\frac{1}{2}{{\left\| w \right\|}^2} + \frac{1}{{vl}}\sum\limits_{i = 1}^l {{\varepsilon _i}} - \rho } \right\} \\ s.t:w\phi \left( {{x_i}} \right) \geq \rho - {\varepsilon _i},i = 1,2,3, \cdots ,l,{\varepsilon _i} \geq 0 \\ \end{array}$$

where w is the hyperplane of the decision, ρ is the bias term, ε_i is the non-relaxation variable, and the meta-parameter v ∈ (0, 1), is used to control the number of samples contained in the hypersphere. The decision function is determined using w and ρ and is computed using the following equation: (51)f(x)=wϕ(x)−ρ$$f(x) = w\phi (x) - \rho$$

1)

General Architecture of DA-LSTM-OS Algorithm

This subsection introduces the two-stage attention-based LSTM autoencoder and a class of support vector machines for anomaly detection algorithm (DA-LSTM-OS). This subsection proposes an LSTM algorithm incorporating the attention mechanism implements an encoder and decoder to extract features, and then inputs the extracted features into the OC-SVM for training, thus realizing the detection of abnormal physiological signals.

The hidden state after the attention mechanism is h˜t$${\widetilde h_t}$$, which is computed by multiplying the hidden state h_i of each encoder by its respective attention weight for a weighted sum, thus mapping h_i to a temporal component of the decoder input as shown in Equation (61). (61)h˜t=∑i=1Tβtiht$${\widetilde h_t} = \sum\limits_{i = 1}^T {\beta _t^i} {h_t}$$

The loss function uses the mean square error (MSE) function: (62)MSE(Y,X)=∑i=1T(yi−xi)2T$$MSE\left( {Y,X} \right) = \frac{{\sum\limits_{i = 1}^T {{{\left( {{y_i} - {x_i}} \right)}^2}} }}{T}$$

1)

Experimental design

The dataset in the experiment for measuring heart rate was derived from the Signal processing cup, which consists of 12 sets of signals collected from subjects between the ages of 20-35 years old, with each set of signals including two PPG signals collected using a bracelet integrated with 515 nm green light, one ECG signal, and three-axis acceleration signals. The triaxial acceleration signals were acquired through an acceleration sensor integrated into the bracelet and included motion information in three mutually perpendicular directions, and the ECG signals were electrocardiogram signals that were acquired through electrode pads immediately adjacent to the subject’s chest, while the reference heart rate provided was calculated from the ECG signals. During the data acquisition process, the subject followed a certain pattern of movement on the treadmill in order to achieve the desired state of exercise, which is as follows: a)

Rest (30 seconds) -> 8 km/h (1 minute) -> 15 km/h (1 minute) -> 8 km/h (1 minute) -> 15 km/h (1 minute) -> Rest (30 seconds).

Each set of samples was recorded for 5 minutes of exercise, the overall sample included both forms of exercise described above, and the signal was sampled at 125 Hz, with the corresponding reference heart rate provided every 2 seconds.

Multi-channel parallel adaptive filtering technique and multi-parameter combination of heart rate algorithms were used to measure heart rate. For the experiment, 8s was used as the time window, and each step was 2s, and the algorithm was utilized to calculate the heart rate in each time window.

Figure 6.

The heart rate and the Pearson correlation coefficient of the heart rate

1)

Experimental design

The dataset for the blood pressure experiment consisted of two aspects, on the one hand, the data came from the MIMIC dataset, which included 80,000 physiological waveform records of the patients, and each set of recordings included one or more ECG signals, arterial blood pressure ABP signals as well as PPG signals, and from which we chose 10 sets of recordings that contained both PPG and ABP, and each set of signals contained one hour of consecutive waveform recordings,. The sampling rate of the signals was 125 Hz.

The experimental data set for blood pressure did not include triaxial acceleration signals, and a low-pass filter was used to remove high-frequency noise. For the experiment, the characteristic parameters of PPG were first extracted, and the blood pressure equation was fitted with 3 data sets as the training set and the remaining 7 data sets as the test set. We extracted the reference blood pressure from the ABP signal, with the peak representing diastolic blood pressure and the trough representing systolic blood pressure. The time window used for the experiment was set to 5s.

1)

Experimental design

The dataset for measuring oximetry was constructed using self-acquired data. PPG signals were collected through a physiological indicator monitoring development board, while oximetry was measured as a reference value using a Myriad IPM6 patient monitor. The dataset consists of four sets of data from adult males, each recorded over a 10-minute period, including two PPG and infrared light acquisitions from the development board. The reference oximetry was read every 5 seconds. For the experiment, the left middle finger was pressed on the development board to acquire the PPG signal while the left index finger wore a finger splint to measure the reference oxygen saturation through a patient monitor. The signal’s sampling rate was 100Hz.

The oximetry data does not include the triaxial acceleration signal, so we use low-pass filtering to remove the noise. In the experiment, the R-value is calculated by calculating the AC and DC components of the PPG signal, and then the quadratic equation of the oximetry is fitted. The training set of the experiment comprises the first 5 minutes of the 4 sets of data, and the test set consists of the final 5 minutes of the data. The time window in the experiment was 5 s, which was consistent with the reference oxygen saturation.

Figure 7.

Reference blood oxygen and forecast blood oxygen contrast

During the operation of the whole system, data will be collected at high frequency at any time. At the same time, in order to deal with the huge amount of calculations, it is also necessary to put forward certain requirements on the computing power of the processing chip. STM32F103 has a powerful core of Cortex-M3 and 32-bit ARM microcontroller can fully meet the system’s demand for computing. Its core number and bus width are the best in its class, and it also has a rich set of pins and interfaces required by this design, which can satisfy the system’s various needs for ROM and RAM.

The PPG signal acquisition module of this system selects the SFH7050 photoelectric sensor manufactured by OSRAM, which can support the acquisition of relevant signals such as heart rate, blood oxygen, and blood pressure. Under the working condition, the sensor can emit two channels of light at the same time, i.e., red light with a wavelength of 660nm and infrared light with a wavelength of 940nm. For the system acquisition accuracy needs, this paper at the same time selected the AFE4400 acquisition front-end, and in the original sensor based on the addition of 535nm green light emitting module, greatly improving the accuracy of the PPG signal acquisition.AFE4400 is mostly used in the biological signal acquisition front-end, the collected signals for filtering, amplification, and A/D conversion, etc., this paper makes improvements to its increased I / V conversion circuit and its filtering and conversion circuit. V conversion circuit and modify some parameters of the filtering and amplifying circuit, which can make the pre-processing of PPG signals to achieve more ideal results.

SFH7050 photoelectric sensor collected for the current signal, can not be extracted from the corresponding waveform and preprocessing operations, before the signal needs to be I / V conversion. Due to the size of the collected current signal for the microampere level, so after the I / V conversion is also required to process the signal filtering and amplification, in order to reduce noise interference, to improve the accuracy of the purpose.

The converted voltage value V₀ is calculated from equation (63): (63)V0=i1×Rf$${V_0} = {i_1} \times {R_f}$$

The input current value is i₁, R_f is the feedback resistor in the figure. After the output voltage value V₀ can be input to the subsequent filtering, amplification circuit.

The voltage signal output by the I/V conversion circuit is accompanied by a large number of noise signals, including: baseline drift, industrial frequency interference, ambient light and motion artifacts. Baseline drift is generally caused by breathing, and its frequency is generally lower than 0.5 Hz. Industrial frequency interference is mostly from the operating power supply, the noise signal frequency is generally around 50 Hz. Ambient light interference is due to the photosensitive element does not work in a closed environment and thus cause interference to the PPG signal, through the improvement of the probe structure and other technologies, the impact of ambient light can be improved.

The Butterworth second-order low-pass filter has the best amplitude-frequency characteristic curve when quality factor Q = 0.707 is used. Multi-stage cascading of the filter module is used to improve the filter roll-off characteristics so that it can filter out noise signals other than 5 Hz at one time, which greatly improves the efficiency.

By Eq: (64)Q=13−AUP=0.707$$Q = \frac{1}{{3 - {A_{UP}}}} = 0.707$$ (65)AUP=3−1Q=3−10.707≈1.586$${A_{UP}} = 3 - \frac{1}{Q} = 3 - \frac{1}{{0.707}} \approx 1.586$$

Set the cutoff frequency to 24 Hz, i.e: (66)fH=12πRC=24Hz$${f_H} = \frac{1}{{2\pi RC}} = 24Hz$$

Usually, capacitor C is used with a capacitance of 0.1μF. The amplification circuit in the second half of the module uses the differential amplification circuit of the OP07CDR chip, which amplifies the signal by a factor of: (67)AU=1+R7R6$${A_U} = 1 + \frac{{{R_7}}}{{{R_6}}}$$

Design R₇ = 10 kΩ, R₆ = 5.86 kΩ, the amplification of the filter after cascading is about 1.6 times.

In this paper, the amplification of the inverse amplification circuit is set to 16 times, which can meet the demand of the subsequent module for signal processing. And the amplification multiplier: (68)AV=V0Vi=−R15R17$${A_V} = \frac{{{V_0}}}{{{V_i}}} = - \frac{{{R_{15}}}}{{{R_{17}}}}$$

According to the working principle of the NLMS adaptive filter, in the adaptive filtering of the preprocessed PPG signal, it is necessary to synchronize the input of three-axis acceleration signals from the same body at the same time at the second input. The sensor can meet the three axial acceleration acquisition at the same time, based on the difference between the front and rear two acceleration to realize the acquisition of motion signals. If the difference between the two signals is higher than the threshold, the corresponding pin will be cut off. The four-wire SPI interface can be configured in master-slave mode. In this system, the STM32 master chip can be set as a master mode, and the LIS3DH sensor can be set as a slave mode.

In the process of developing health monitoring system, due to the requirement of its standby time, it is necessary to reduce the power consumption of the equipment under the premise of meeting all the functions of the equipment, so as to achieve the purpose of long standby time, energy saving and emission reduction. Therefore, we should not only put forward the requirements of low power consumption from the system level, but also reduce the power consumption from each module and detail, so as to improve the performance and standby time of the device.

The system software architecture is shown in Fig. 9. The figure mainly introduces two timers and sets the timing period to 10ms and 2ms respectively to interrupt the normal working mode. Before executing the program, initialize it to preset the working mode. Set the PWM wave output of the timer to stabilize the current in the light source circuit at 20mA and set the frequency of the PWM wave to 100KHz. if the master chip does not execute the program related to filtering, it switches to the sleep mode and turns on the timer, which can be woken up by stopping the timer and switching to the fast operation state. When the master chip is interrupted from the 10ms timer, it will start filtering DC components with frequencies below 0.5Hz and return the processed components to the input. In the second half, by setting a timer with a time period of 2ms, the main control chip can be switched to sleep mode, in which synchronization between the I/V conversion circuit and the PPG signal of the signal amplifier circuit will be turned on. After 2ms, the main control chip is awakened by the timer and performs the secondary amplification and filtering process, while simultaneously turning on the acquisition of the three-axis acceleration signal. After the acquisition is completed and the processed PPG signal is divided into two separate inputs to the adaptive filter, the output waveform is decomposed and reconstructed to execute the health extraction algorithm. The final calculated health value and pulse wave are transmitted to the terminal platform of the smart device through Bluetooth module for processing. Once all the processes are completed, the main control chip will switch to low-power mode until the next 10ms timed interrupt occurs.

Figure 9.

System software architecture diagram

The calculation of peak-to-peak intervals and the number of peaks per minute have mostly been used for the calculation of heart rate values from ECG signals or PPG signals at rest. As shown in Eq. (69), some scholars have validated the estimation of heart rate by putting the PPG signal processed by hardware filtering and software algorithms into the frequency domain f_es. (69)fes=60×fHR$${f_{es}} = 60 \times {f_{HR}}$$

where f_HR denotes the true heart rate.

In this paper, the position index P_n of the estimated heart rate in the frequency domain is obtained by spectral peak tracking, and the estimated heart rate can be obtained by substituting into Eq. (70): (70)fes=60×PnNfs$${f_{es}} = 60 \times \frac{{{P_n}}}{N}{f_s}$$

where f_s denotes the sampling frequency and N denotes the number of grids for [0,fs]$$\left[ {0,{f_s}} \right]$$. The grid frequency is set to be about 1 time/second, i.e: (71)60fsN=1$$60\frac{{{f_s}}}{N} = 1$$

where the number of grids N needs to be satisfied N ≫ σf_s, σ is the duration of each window. The FFT transform of the PPG time-domain signal and spectral peak tracking are performed to obtain the position index of the estimated heart rate in the frequency domain P_n.