A new single-channel fetal heart detection algorithm

Abstract

The application provides a new type of single-channel-based fetal heart detection algorithm, first uses a single-channel body surface bioelectricity collection device to collect pregnant woman abdominal bioelectricity data, wherein the data contains fetal electrocardiosignal and other mother and baby physiological signals; then uses wavelet transform and singular value decomposition algorithm to reconstruct the pregnant woman reference electrocardiosignal, uses an adaptive filtering algorithm to eliminate most of the interference physiological signals; then uses wavelet transform, singular value decomposition and digital band-pass filter to reconstruct the fetal electrocardiosignal, obtains the fetal electrocardiogram curve, and simultaneously uses time-frequency domain transform algorithm, non-negative matrix decomposition algorithm and power spectral density to obtain accurate fetal heart rate. The algorithm is a single-channel non-invasive passive fetal electrocardiosignal analysis algorithm, is convenient to operate, is reliable in performance, can accurately obtain fetal electrocardio information, and provides a basis for non-invasive fetal heart monitor product design based on body surface bioelectricity signal collection and analysis.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of biomedical signal processing, and specifically relates to a novel fetal heart detection algorithm, which realizes single-channel non-invasive passive fetal electrocardiogram signal detection and analysis. BACKGROUND

[0002] Fetal electrocardiogram signal is a source signal generated by the heart of a fetus in the development stage, reflects the development of the fetus in the uterus, can help doctors judge whether the fetus has defects in development, is an important means for effective monitoring of the health of the fetus during pregnancy, and has very important significance in the medical field.

[0003] At present, the mainstream fetal electrocardiogram detection scheme is active invasive detection, and representative methods include Doppler ultrasound method, scalp electrode method and the like. These methods have high accuracy, but in the process of signal acquisition, they may cause some potential harm to the pregnant woman or the growth and development of the fetus. At the same time, these detection instruments are usually large in size and need to be operated by professional personnel, which is not conducive to real-time monitoring of fetal heart, and there is still a problem that diseases cannot be found as early as possible. Although there are small fetal heart instrument products on the market that use ultrasonic detection principles, they can meet the needs of part of the family for daily monitoring, but there are a series of problems such as poor sensitivity and poor precision.

[0004] The maternal abdominal wall electrode method is an important development direction of fetal electrocardiogram signal detection in the future, which places one or more electrodes at different positions on the abdomen of a pregnant woman, collects the body surface bioelectric signal, and obtains the fetal electrocardiogram signal by using a separation algorithm. The operation process of this method is similar to adult electrocardiogram detection, and since it is non-invasive passive detection, it will not have any impact on the development of the fetus, and has very high application value. However, the development of this technology currently mainly faces the following problems:

[0005] Firstly, since the fetal electrocardiogram signal is usually weaker than other interfering physiological signals, the mixed signal obtained has rich noise types and low signal-to-noise ratio, which increases the difficulty of overall recognition, detection and extraction, and the accuracy of the obtained results is low;

[0006] Secondly, most researches or patents need to use abdominal multi-channel signals (also called "multi-abdominal lead") to complete heart rate analysis, and there are problems of complicated operation and easy to make mistakes in the process of product miniaturization and marketization;

[0007] Thirdly, most researches and patents use neural networks and other means to ensure the accuracy of detection, which has high requirements for the design of product hardware circuit, and is often accompanied by the problem of high price cost, which is not conducive to the marketization of related products.

[0008] Therefore, it has practical significance and good application prospect to design a fetal heart detection scheme which is simple to operate, has high accuracy and low manufacturing cost. SUMMARY

[0009] The present application aims to overcome the problem that the existing fetal heart monitor on the market is difficult to balance the operation simplicity, high accuracy and low manufacturing cost, and provides a single-channel non-invasive passive fetal electrocardio signal detection and extraction method, which can better balance the three, and provides a high-value scheme for ordinary families to monitor fetal electrocardio in daily life, and helps to find and treat the disease as soon as possible.

[0010] In order to achieve the above-mentioned purpose, the specific disclosed technical scheme of the present application mainly includes the following steps:

[0011] S1, a single-channel sampling device collects the abdominal skin bioelectric signal A of a pregnant woman (also called "single abdominal lead" signal).

[0012] S2, the signal A is sent into a wavelet transform-singular value decomposition reconstructor to obtain a reference electrocardio signal B of the pregnant woman. The symN series wavelet basis function is used in the wavelet transform link, and the purpose is to remove baseline noise, electromyographic noise and power frequency interference. For the use of wavelet basis function, the traditional discrete wavelet transform (such as Mallet algorithm) can also remove the three kinds of noise, but in this application scenario, the Mallet algorithm will introduce a relatively serious pseudo-Gibbs phenomenon, which will greatly damage the details of the electrocardio signal, so the traditional discrete wavelet transform method is replaced by a stationary wavelet transform method, then N-layer symN wavelet decomposition is used, and the N-layer approximation coefficients are set to zero after reconstruction, so that the signal without baseline noise can be obtained. The signals obtained by wavelet decomposition reflect the characteristics of the source signal in a certain frequency range, and by filtering the coefficients of a certain layer, the characteristics of the signal in that frequency range can be changed, which is the theoretical basis of wavelet denoising. In this scenario, the signal without baseline noise is subjected to 5-layer symN stationary wavelet decomposition, and the hard threshold function and the VisuShrink improved threshold combination are used for denoising of the detail coefficients of each layer, and finally all the approximation coefficients and detail coefficients are used for reconstruction to obtain the target signal, and the filtering operation is completed.

[0013] S3, the reference electrocardio signal B of the pregnant woman and the abdominal skin bioelectric signal A of the pregnant woman are sent into an adaptive filter, and most of the electrocardio components of the pregnant woman are eliminated by using the recursive least square method to obtain a noisy fetal electrocardio signal C with good signal-to-noise ratio.

[0014] S4, a wavelet transform-singular value decomposition reconstructor is used to filter the residual interference signals in signal C and compensate for the lost detail signals to obtain signal D.

[0015] S5, a digital band-pass filter is used to further filter and shape signal D to obtain a fetal electrocardio curve; at the same time, Choi-Williams time-frequency distribution transformation, non-negative matrix decomposition and power spectral density algorithm are used to calculate the accurate heart rate of the fetus.

[0016] Compared with the prior art, the present application has the advantages that the present application overcomes the problem of low accuracy in the single-rib guide case, simplifies the user operation process, and can well maintain the detection accuracy in the case of low signal-to-noise ratio. Meanwhile, the traditional signal processing algorithm is used, and the embedded chip can be programmed and some common digital signal processing devices can be used in the product manufacturing process, without the need to re-produce a special circuit, thereby reducing the development cost. BRIEF DESCRIPTION OF DRAWINGS

[0017] Figure 1 It is an overall structure diagram of the algorithm of the present application.

[0018] Figure 2 It is a measured pregnant woman abdominal lead signal waveform diagram provided by the DaISy database (the DaISy database is established by a Belgian scholar Lathauwer, and the example in the figure uses the No. 3 abdominal guide clinical data);

[0019] Figure 3 It is a fetal electrocardiosignal waveform diagram finally extracted by the present application;

[0020] Figure 4 It is a fetal electrocardiosignal local amplification waveform diagram finally extracted by the present application;

[0021] Figure 5 It is an iteration update algorithm flowchart used in the non-negative matrix factorization algorithm of the present application;

[0022] Figure 6 It is a fetal electrocardiosignal power spectral density diagram calculated by the present application. DETAILED DESCRIPTION

[0023] The implementation of the present application will be further described in detail in combination with the drawings and specific implementation methods, and the following examples or drawings are used to illustrate the present application, but not to limit the scope of the present application.

[0024] Please refer to Figure 1The embodiment provides a complete algorithm flow structure of the application. The Chinese meanings and common abbreviations of the key English words are as follows: amplitude, frequency, time, electrocardiogram, fetal electrocardiogram, wavelet transform, singular value decomposition, adaptive filter, band pass filter, Choi-Williams distribution, non-negative matrix factorization, and power spectral density. If no special description is given, the English words and abbreviations appearing in the following text refer to the Chinese meanings herein.

[0025] The detection flow of the fetal electrocardiogram signal of the application is as follows:

[0026] S1, signal collection, specifically: using a single-channel body surface bioelectricity collection device to collect the body surface electrical signal of a pregnant woman, and recording it as signal A, and the collection result is shown in Figure 2 .

[0027] S2, reference signal reconstruction, separating the electrocardiogram signal B of the pregnant woman from the body surface bioelectricity signal A, specifically:

[0028] S2-1, first, the signal A is sent into a wavelet transform module to remove most of the baseline drift, myoelectricity signal and power frequency interference noise and the like in the signal. Experiments show that N=8 or 10 in the symN wavelet base function can achieve good results, and in the embodiment, the N value is set to 8 by default, that is, the sym8 wavelet base function is used, and the following will not be described. In the module, first, the input signal is subjected to 8-layer sym8 stationary wavelet decomposition, and the 8th layer approximation coefficient is subjected to zero processing, and then the intermediate signal A1 is obtained by gradually reconstructing. Then, the signal A1 is subjected to 5-layer sym8 stationary wavelet decomposition, and the VisuShrink improved threshold is used for hard threshold filtering of the detail coefficients of each layer, and then the target signal A2 is obtained by gradually reconstructing the approximation coefficients and the detail coefficients. The definition formula of the VisuShrink improved threshold is as follows:

[0029]

[0030] where δ represents the standard variance of noise, N represents the length of wavelet coefficient acted by threshold, and j represents the decomposition scale of wavelet coefficient acted by threshold (also referred to as decomposition layer, same below). The improvement of the threshold is that the traditional VisuShrink threshold is proposed under the Gaussian model for the joint distribution of multi-dimensional independent normal variables, which is a global processing method. The improved threshold considers the characteristics that noise decreases with the increase of wavelet decomposition scale, adjusts the threshold used for each layer, prevents the damage of inappropriate threshold to useful components in the signal, is a local processing method, and is more suitable for non-stationary signal processing represented by an electrocardiogram signal.

[0031] The hard threshold function calculation formula is:

[0032]

[0033] where λ is the set threshold. The meaning of the formula is that when the absolute value of the signal X is greater than or equal to the threshold, it can be passed without loss, and when it is less than the threshold, it will be zeroed, that is, not allowed to pass. In this way, the wavelet filtering operation is completed.

[0034] S2-2, input the signal A2 into the source signal channel of the SVD reconstruction module, and input the signal A into the reference signal channel of the SVD reconstruction module. The working principle of the module is the SVD singular value spectrum theory: large singular values correspond to signal components with large energy, and small singular values correspond to signal components with small energy. As long as the singular values of the corresponding components are zeroed, the components can be removed. Therefore, the signal processing method of the SVD reconstruction module is: first, estimate the maternal heart rate value in the signal A using the power spectral density, from which the number of points in a period can be obtained; then, the one-dimensional signal A2 is split into a signal matrix according to the period, and SVD operation is performed thereon, and the calculation formula is:

[0035] A2=U∑V T

[0036] where U and V are orthogonal matrices, ∑ is a singular value matrix, and is arranged in descending order. In the obtained singular value matrix, the smallest non-zero singular value is zeroed, then inverse SVD operation is performed, and the matrix is restored to a one-dimensional target signal A3.

[0037] S2-3, input the signal A3 into the wavelet transform module, and the operation of this step is the same as that of S2-2, to obtain the reconstructed maternal reference electrocardiogram signal B.

[0038] S3, adaptive cancellation filter, remove most of the pregnant ECG components from signal A, specifically: the pregnant reference ECG signal B is input to the reference signal channel of the adaptive filter, the collected signal A is input to the source signal channel of the adaptive filter, then the recursive least squares method is used to realize the adaptive filtering process, eliminate the pregnant ECG components in the collected signal A, obtain the noisy fetal ECG signal C.

[0039] The update formula of the adaptive filter weight used in this step is:

[0040]

[0041] Wherein, The adaptive filter weight at the n th moment is a p x 1 column vector, wherein p is the length of the filter weight and p = 5; The adaptive filter weight at the n-1 th moment is a p x 1 column vector; K[n] is a p x 1 column vector gain factor, the greater the gain factor, the greater the correction of e[n] is the error, and the calculation formula is:

[0042]

[0043] Wherein, h[n] is a reference data vector, which is expressed as:

[0044] h[n] = [x R [n] x R [n-1] … x R [n-p+1]] T

[0045] Wherein, x R [n] is a reference signal, which is the reference ECG signal B of the pregnant woman in this patent. The update formula of the gain factor K[n] is:

[0046]

[0047] Wherein, ∑[n-1] is the covariance matrix of λ is a forgetting factor, and λ = 0.99, which reduces the influence of previous samples on the next moment estimation. The update formula of ∑[n] is:

[0048] ∑[n] = (I-K[n]h T [n])∑[n-1]

[0049] Finally, the filter zero input state weight and covariance matrix are:

[0050]

[0051] ∑[-1] = 105 I

[0052] where ∑[-1]=10 5 The selection of I is to reduce the influence of the weight and the covariance matrix in zero-input state on the subsequent estimation.

[0053] S4, signal filtering and compensation, further improve the signal-to-noise ratio of the signal, specifically:

[0054] S4-1, input the signal C into the wavelet transform module, the operation of this step is the same as S2-2, and the intermediate signal C1 which removes the baseline noise, electromyographic noise and power frequency interference is obtained.

[0055] S4-2, input the signal C1 into the source signal channel and the reference signal channel of the SVD reconstruction module, and obtain the intermediate signal C2 which removes the maternal electrocardiogram signal. According to the singular value spectrum theory, it is necessary to set the larger singular value to zero to remove the noise. Therefore, the processing of the SVD reconstruction module to the signal is: calculate the rough fetal heart rate value in the signal C1 using the power spectral density, and divide and recombine the one-dimensional signal C1 according to the period to form a signal matrix; then perform SVD operation on the matrix, and set the largest three singular values in the singular value matrix to zero; finally, perform inverse SVD operation, and restore the matrix to one-dimensional signal C2.

[0056] S4-3, perform wavelet transform on the obtained signal C2 according to the steps of S2-2, and the final target signal D can be obtained. This is because in the previous processing, some detail components of the signal will be distorted and additional high-frequency noise will be introduced, so the wavelet transform in this step is needed to perform certain waveform compensation.

[0057] S5, extract the fetal electrocardiogram curve from the signal D obtained from S4, specifically:

[0058] Use the digital elliptical band-pass filter to filter out the burr in the signal, and the fetal electrocardiogram curve can be obtained. The design of the digital elliptical band-pass filter adopts the commonly used bilinear transformation method, and the integral step T=2 is set, and the corresponding index in the analog frequency domain is: the lower stop band edge frequency f s1 =3Hz, the passband start frequency f c1 =5Hz, the passband cutoff frequency f c2 =82Hz, the upper stop band start frequency f s2 =100Hz, the lower stop band ripple fluctuation α s1 =60dB, the upper stop band ripple fluctuation α s2 =80dB, the passband ripple fluctuation α p =1dB. The specific extraction effect of the fetal electrocardiogram curve can be referred to Figure 3 , and the local detail amplification effect of the fetal electrocardiogram signal extraction result can be referred to Figure 4 .

[0059] S6, the signal D obtained by S4 is sequentially subjected to CWD transformation, NMF operation and PSD operation, so that the fetal heart rate can be accurately calculated in an application scenario with low signal-to-noise ratio, and specifically:

[0060] S6-1, using CWD transformation, the one-dimensional signal D is changed into a high-dimensional intermediate signal D1 to fully expose the instantaneous frequency information in the signal. Time-frequency analysis is a common means of analyzing non-stationary signals, which can clearly depict the relationship between frequency and time, and can detect small signals. In the field of biological signal detection, it has very important application in strong noise scene signal extraction and analysis. CWD transformation is a common time-frequency analysis transformation, which can weaken the cross terms generated in the traditional Wigner-Wille distribution, reduce the interference of artifacts on the real signal, and also improve the resolution of the self term and improve the detection accuracy of weak signals. The discrete form of CWD is as follows:

[0061]

[0062] Where τ and ω represent time delay and angular frequency respectively, σ is a controllable factor, τ is time delay, and the corresponding kernel function is:

[0063]

[0064] S6-2, using NMF algorithm to decompose the high-dimensional intermediate signal D1, a low-dimensional one-dimensional signal D2 is obtained. NMF algorithm is a multivariate data analysis method, which can approximately decompose a large non-negative matrix into two smaller non-negative matrices through iterative algorithm, so as to achieve the purpose of dimensionality reduction and feature extraction. The definition of NMF is as follows:

[0065] V m×n =W m×r H r×n +E m×n

[0066] Where V m×n represents the source signal matrix, all elements are non-negative; W m×r represents the decomposition base matrix; H r×n represents the extracted source signal feature matrix; E m×n is the error matrix, which reflects the deviation between the decomposition result and the source signal; the coefficient r should satisfy The problem of finding the NMF decomposition expression can be equivalent to an optimization problem:

[0067]

[0068] The meaning of the expression is to find a pair of W and H matrices, under the condition that both W and H matrices remain non-negative, so that the L2 norm of the error matrix E produced by the decomposition is minimum, that is, the decomposition error is minimum. In actual operation, the gradient descent method is often used on the basis of the Euclidean distance model to find the optimal solution, and this idea is also called the multiplication iteration algorithm, and its iteration expression is as follows:

[0069]

[0070]

[0071] After the iterative decomposition is completed, from the perspective of characteristic domain approximation equivalence, the embodiment selects the H r×n second row row vector as the target low-dimensional one-dimensional signal D2, and thus the NMF operation ends. The flowchart of the NMF multiplication iteration update algorithm used in the embodiment can be referred to in Figure 5 .

[0072] S6-3, calculate the fetal heart rate using the PSD algorithm. Calculate the PSD of the signal D2 to obtain the power spectrum density graph; by positioning the highest peak in the graph except near the zero frequency, the fetal heartbeat frequency corresponding to the peak value can be obtained; multiplying the value by 60 can calculate the fetal heart rate. Thus, the entire detection process is completed. The PSD calculation result of the embodiment is 131 times / min, compared with the true result 132 times / min, the deviation is within the normal error range, indicating the effectiveness and feasibility of the algorithm. The PSD extraction result can be referred to in Figure 6 .

[0073] Those skilled in the art will realize that the embodiments described herein are for the purpose of helping the reader to understand the principles of the present application, and should be understood as not limiting the protection scope of the present application to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations according to the technical inspirations disclosed in the present application without departing from the essence of the present application, and these modifications and combinations are still within the protection scope of the present application.

Claims

1. A novel fetal heart rate detection algorithm, comprising, in sequence, a wavelet-singular value-wavelet reconstruction module, an adaptive filtering module, a wavelet-singular value-wavelet compensation module, and a feature extraction and calculation module, characterized in that: First, the wavelet-singular value-wavelet reconstruction module is used as follows: a single-channel surface bioelectrical signal acquisition device is used to acquire the abdominal bioelectrical signal of the pregnant woman; the acquired raw data includes fetal electrocardiogram signal and other interfering physiological signals, including but not limited to the pregnant woman's electrocardiogram signal, electromyography signal, etc.; a relatively good reference electrocardiogram signal of the pregnant woman is reconstructed by wavelet transform-singular value decomposition-wavelet transform in sequence. Secondly, the adaptive filtering module is used as follows: the original sampled signal and the pregnant woman's reference ECG signal are input into the adaptive filter, and most of the pregnant woman's ECG components are eliminated by the recursive least squares method to obtain a noisy fetal ECG signal with a good signal-to-noise ratio. Next, the wavelet-singular value-wavelet compensation module is used as follows: the wavelet transform-singular value decomposition-wavelet transform method is used sequentially to filter and compensate the output signal of the adaptive filtering module; Finally, the feature extraction and calculation module is used as follows: a digital bandpass filter is used to shape the output signal of the wavelet-singular value-wavelet compensation module to obtain the fetal electrocardiogram curve; at the same time, the accurate fetal heart rate is calculated by sequentially using the Choi-Williams time-frequency distribution transform, non-negative matrix factorization and power spectral density algorithm.

2. The novel fetal heart rate detection algorithm for pregnant women as described in claim 1, characterized in that, The wavelet transform part is performed according to the following steps: First, the signal is decomposed into N-level symN stationary wavelet decompositions, and the approximation coefficients of the Nth level are set to zero; the new signal is then reconstructed step by step. Finally, the signal is subjected to 5-level symN stationary wavelet decomposition, and hard thresholding is performed on the detail coefficients of each level using the VisuShrink improved threshold. The target signal is then reconstructed step by step using the approximation coefficients and detail coefficients at each level. The VisuShrink improved threshold is defined as follows: Where δ represents the standard deviation of the noise, N represents the length of the wavelet coefficients affected by the threshold, and j is the decomposition scale of the wavelet coefficients affected by the threshold.

Images