CN105997064A - Method for identifying human lower limb surface EMG signals (electromyographic signals) - Google Patents

Abstract

The invention relates to a method for identifying human lower limb surface EMG signals (electromyographic signals). The method comprises the following main steps: acquiring surface EMG signals of corresponding muscle blocks under the stimulation of movement actions in real time; carrying out pretreatment on the collected EMG signals, thus obtaining the EMG signals with artifact signals eliminated; carrying out decomposition on the obtained EMG signals by adopting a discrete wavelet transform method, thus obtaining a low-frequency coefficient vector and a high-frequency coefficient row; carrying out singular value decomposition on obtained wavelet components by adopting a filtering method combining time domain with frequency domain, and constituting a characteristic matrix by adopting the singular values obtained by decomposition; carrying out training on a characteristic sample by adopting a support vector machine, and generating a support vector machine classifier for carrying out classification and identification on a blind sample. The method has the beneficial effects that compared with the prior art, the novel lower limb EMG signal preprocessing method is provided.

Description

A kind of discrimination method for human body lower limbs surface electromyogram signal

Technical field

Surface electromyogram signal (the Electromyographic that the present invention is designed in the bioelectrical signals of human body lower limbs Signal, EMG) artefact eliminate, feature extraction and identification technique field.

Background technology

Data show according to statistics, and China is the most formally marching toward aging society at the beginning of 21 century, and aging process exceedes it He is national, it is contemplated that the year two thousand twenty China aging population can reach 2.48 hundred million, and the year two thousand fifty is up to 400,000,000.Old people and extremity disabled persons Population structure expands the most rapidly, the distinguishing feature of above-mentioned crowd be its daily behavior activity be both needed to provide assist help. Paralysis is one of most common reason causing above-mentioned crowd's LOM, and especially lower part of the body paralysis, it relates to limbs, body Dry part or afunction completely.At present clinical expert it is believed that limb motion rehabilitation be considered as one effectively Solution, it requires that affected lower limb actively assists in actively exercise.But, obstacle is existed for limb function Crowd, be generally difficult to smoothly complete such as stand, squat down, the lower limb motion such as walking.Therefore, grinding by EMG signal Study carefully the biofeedback mechanism that can help to explore nervus motorius with muscular tissue, it was predicted that with perception limb motion situation, assessment is old Year people, people with disability and the musculation ability of sub-health population, be applicable to the health of the elderly and physically disabled for development Purgation again limb exercise aid device provides theoretical foundation and application foundation.

Surface EMG signal is a kind of physiology letter analyzing the human body lower limb athletic performance relevant to activities of daily living Breath source, this signal be by electrode guide, amplify and from the neuromuscular system activity that muscle surface is recorded time non- Steadily One-dimension Time Series bioelectrical signals, it can reflect the motion feature of muscle strength and people.The motion intention of people is usual The contraction and the diastole that are stimulated muscle cell by neural excitation are realized, owing in different limb motions, the pattern of muscle contraction is not With, cause the feature of corresponding surface electromyogram signal also to have difference, at the surface EMG signal energy being intended under control produce of people Reflected well limb motion or motion characteristic, assess human body motion intention.Limb is mainly passed through in the realization that human motion is intended to Body completes, and owing to lower limb EMG signal is increasingly complex relative to upper limb, simultaneously by bigger noise jamming, therefore transports human body Dynamic research is concentrated mainly on upper limb EMG signal, and the research to lower limb EMG signal and identification thereof needs the most perfect.Mesh Front in the identification of EMG signal, conventional method is all based on traditional classification or clustering algorithm, such as support vector machine, nerve Network (Neural Network Algorithm is called for short " NNA "), linear discriminant analysis (Linear Discriminant Analysis, is called for short " LDA ") etc..Wherein, LDA algorithm i.e. can identify certain single action, it is also possible to multiple actions is added Upper label is identified as a special class.In the research of lower limb EMG signal, 2016, John A.Spanias et al. Use LDA algorithm, have studied and only carry out classifying and by EMG signal and the other types of instrument sensor return with EMG signal The data method that carries out together classifying；2014, AJ Young et al. used the method pair of Sensor Time History EMG signal is classified, but the method only account for whole during the time span of signal.

Summary of the invention

The purpose of the present invention, it is simply that for the problems referred to above, it is provided that a kind of identification for human body lower limbs surface electromyogram signal Method.

By experimental analysis, find that lower limb EMG signal yardstick is the faintest, simultaneously because hardware limitations and limbs The reasons such as movement, are the most extremely easily subject to power frequency, baseline drift and white Gaussian noise interference, therefore test the lower limb collected EMG signal artefact is serious, directly carries out feature extraction based on initial data and Classification and Identification is unpractical.Tradition artefact eliminates Method is primary signal to carry out bandpass filtering treatment (the main component integrated distribution of lower limb EMG signal is in 20～500Hz frequency ranges On).The present invention carries out notch filter and low-pass filtering to Hz noise, baseline drift respectively, the work proposed according to the present invention Frequently lower limb EMG signal is estimated by the interference noise factor and baseline drift noise factor respectively, when the energy of noise is beyond threshold During value, use trap and low-pass filtering, otherwise use filtering front signal.And for white Gaussian noise, first calculate lower limb EMG letter Number zero passage count, counted by zero passage carries out interval division to lower limb EMG time-domain signal, then by the white Gaussian noise factor by The secondary noise to each interval is estimated, and when the energy of noise is beyond threshold value, uses filter result, before otherwise using filtering Signal.On the premise of surface EMG signal has non-stationary, non-linear behavior, the present invention is by wavelet transform (Discrete Wavelet Transform, DWT) is applied in the feature extraction of lower limb EMG signal, and it is by time domain, frequency domain Analyze and combine the time to surface EMG signal and information that frequency is comprised is analyzed.And traditional artefact eliminates (such as band Bandpass filter), feature extraction (such as Fourier transform, time and frequency domain analysis) method is the most independently in time domain or frequency domain Analytical data, and EMG signal is considered as steadily or short-term stationarity signal processes, the most traditional method can not be accurate Ground obtains the EMG physiological reaction feature of human body lower limbs limb motion action.Singular value decomposition (Singular Value Decomposition, SVD) it is a kind of effective algebraic characteristic extracting method, owing to singular value features is describing signal numerically More stable, and there is the critical natures such as transposition invariance, invariable rotary shape, shift invariant, therefore singular value features is permissible A kind of effective algebraic characteristic as signal describes.Finally, the present invention is by wavelet transform (DWT) and singular value decomposition (SVD) combine, based on new lower limb EMG signal preprocess method, it is further proposed that time, the filtering method that combines of frequency domain, root According to eigenmatrix obtained above, use the method for support vector machine (Support Vector Machine, SVM) that signal is entered Row Classification and Identification.

The technical scheme is that a kind of discrimination method for human body lower limbs surface electromyogram signal, it is characterised in that Comprise the following steps:

A. putting viscous for disposable electromyographic electrode to lower hind limb musculature epidermis, Real-time Collection correspondence muscle masses move in activity Make the surface EMG signal under stimulating；

B. the EMG signal gathered in step a is carried out pretreatment, it is thus achieved that eliminate the EMG signal after artefact signal；Described pre- Processing method includes Hz noise filtering, baseline drift filtering and white Gaussian noise filtering；In this step, Hz noise filtering, The order of baseline drift filtering and white Gaussian noise filtering can carry out arbitrary arrangement；

C. the EMG signal using discrete small wave converting method to obtain step b decomposes, and obtains low frequency coefficient after decomposition Vector cA1 and high frequency coefficient vector cD1；Low frequency coefficient vector cA1 is decomposed by the method using wavelet transform, it is thus achieved that Low frequency coefficient vector cA2 and high frequency coefficient cD2；Continue the method reusing wavelet transform low frequency coefficient vector is carried out Decompose, until till obtaining low frequency coefficient vector cA5 and 5 high frequency coefficient row cD1, cD2, cD3, cD4, cD5；

When d. using, frequency domain combine filtering method, by step c obtain Wavelet Component carry out singular value decomposition, and The singular value constitutive characteristic matrix that decomposition is obtained；

E. with the eigenmatrix of acquisition in step d as sample, use support vector machine that feature samples is trained, and raw Become support vector machine classifier for blind sample is carried out Classification and Identification.

2, a kind of discrimination method for human body lower limbs surface electromyogram signal according to claim 1, its feature exists In, the concrete grammar of pretreatment described in step b includes:

B1. Hz noise filtering；Method particularly includes: using EMG signal c (t) collected in step a as being originally inputted Signal C_PLIT () carries out notch filter, obtaining filter result is a (t), defines Hz noise factor ε_PLI, then ε_PLICan be by as follows Formula 1 calculates:

Wherein, var is signal variance operator, is used for calculating seasonal effect in time series variance, by Hz noise factor ε_PLITo filter Ripple result is modified shown in equation below 2:

Wherein, s_PLIT () is the final result filtering Hz noise noise, formula 2 shows if the energy of Hz noise noise Amount accounting exceedes the 10% of original energy, then use the filter result of notch filter；

B2. baseline drift filtering；Method particularly includes: the signal s that step b1 is obtained_PLIT () is as original input signal c_BW T () carries out low-pass filtering, filter result is d (t), and the definition baseline drift factor is ε_BW, then ε_BWCan be calculated by equation below 3:

According to baseline drift factor ε_BWBaseline drift noise d (t) obtaining filtering is modified obtaining b (t), and it is expressed Formula is as shown in Equation 4:

Finally, the signal s after removing baseline drift noise is obtained_BWShown in (t) equation below 5:

s_BW(t)=c_BW(t)-b (t) (formula 5)；

B3. white Gaussian noise filtering；Method particularly includes: the signal s that step b2 is obtained_BWT () is as original input signal c_WGNT () is filtered；To time-domain signal c_WGNT () carries out interval division based on zero crossing, then by thresholding method to letter Number c_WGNT () carries out white Gaussian noise filtering；Particularly as follows: comprise j zero crossing z in hypothesis signal_i(i=1,2 ..., j), then For interval z_i＜ t_j＜ z_i+1Interior signal has:

Wherein, T is the threshold value of signal c (t), t_mFor the extreme point in interval；Threshold value T is true by equation below 7 and formula 8 Fixed:

σ=median (| c_wGN(t) |: t=1,2 ..., L)/0.6745 (formula 8)

Herein, σ is input signal c_WGNThe noise level of (t)；L is signal c_WGNT the length of (), in particular for discrete letter Number L is discrete point and counts；Median is signal median operator, is used for obtaining seasonal effect in time series median.

Beneficial effects of the present invention is, relative to conventional art, the present invention proposes a kind of new lower limb EMG signal and locates in advance Reason method, is simultaneously based on the lower limb EMG signal feature after artefact eliminates, and first proposes to combine DWT and SVD, transports lower limb typical case The EMG signal of dynamic action is such as walked, and supports, walking are walked, squatting down and stand etc. carries out feature extraction, finally by SVM to being carried The eigenmatrix gone out carries out identification of classifying.

Accompanying drawing explanation

Fig. 1 is the flow chart of the present invention；

Fig. 2 is to utilize wavelet function that signal is carried out the logic diagram of 5 layers of wavelet decomposition process.

Detailed description of the invention

Below in conjunction with the accompanying drawings, technical scheme is described in detail:

As it is shown in figure 1, the present invention is for the discrimination method of human body lower limbs surface electromyogram signal, mainly comprise the steps that

Step 1) acquisition of lower limb movement raw EMG signal:

Putting viscous for disposable electromyographic electrode to lower hind limb musculature epidermis, Real-time Collection correspondence muscle masses are in active actions Surface EMG signal c (t) under Ci Jiing；

Step 2) primary signal pretreatment (artefact removing method):

The raw EMG signal collected is carried out artefact elimination, it include to Hz noise noise, white Gaussian noise, Filtering of baseline drift noise.The filtering of three kinds of noises in this step does not has fixing sequencing, presses Hz noise here The filtering of noise, the filtering of white Gaussian noise, this order of filtering of baseline drift are described；Secondly, in narration herein In " denoising " and " filtering to noise " implication identical.

A) denoising of Hz noise: the surface EMG signal that will collect, c (t), original as this denoising process Input signal, c_PLIT (), carries out notch filter, obtain filter result, a (t), simultaneously the definition Hz noise factor, ε_PLI, such as public affairs Shown in formula (1):

${\mathrm{\ϵ}}_{PLI}=\frac{\mathrm{var}\[c_{PLI}\left(t\right)-a\left(t\right)\]}{\mathrm{var}\[c_{PLI}\left(t\right)\]}---\left(1\right)$

Wherein, var is signal variance operator, is used for calculating seasonal effect in time series variance.By Hz noise factor ε_PLITo filter Ripple result is modified as follows:

$s_{PLI}\left(t\right)=\left\{\begin{array}{cc}{c}_{PLI}\left(t\right),& {\mathrm{\&epsiv;}}_{PLI}<0.1\\ a\left(t\right),& {\mathrm{\&epsiv;}}_{PLI}\&GreaterEqual;0.1\end{array}\right.---\left(2\right)$

Wherein, s_PLIT () is the final result filtering Hz noise noise.Formula (2) shows if Hz noise noise Energy accounting exceedes the 10% of original energy, then use the filter result of notch filter.

B) denoising of baseline drift: the signal obtained after Hz noise filtering being carried out, s_PLIT (), as this denoising The original input signal of journey, c_BWT (), carries out low-pass filtering, filter result is d (t).Similarly, define the baseline drift factor, ε_BW, as shown in formula (3):

${\mathrm{\&epsiv;}}_{BW}=\frac{\mathrm{var}\&lsqb;d\left(t\right)\&rsqb;}{\mathrm{var}\&lsqb;c_{BW}\left(t\right)\&rsqb;}---\left(3\right)$

According to the baseline drift factor, ε_BW, to filtering the baseline drift noise obtained, d (t), it is modified obtaining, b (t), Shown in its expression formula such as formula (4):

$b\left(t\right)=\left\{\begin{array}{cc}0,& {\mathrm{\&epsiv;}}_{BW}<0.05\\ d\left(t\right),& {\mathrm{\&epsiv;}}_{BW}\&Element;\&lsqb;0.05,0.5\&rsqb;\\ c_{BW}\left(t\right),& {\mathrm{\&epsiv;}}_{BW}>0.5\end{array}\right.---\left(4\right)$

Finally, the signal s after removing baseline drift noise is obtained_BWT () is as follows:

s_BW(t)=c_BW(t)-b(t) (5)

C) denoising of white Gaussian noise: the signal s obtained after baseline drift noise filtering being carried out_BWT () goes as this The original input signal c of process of making an uproar_WGNT () is filtered.To time-domain signal c_WGNT () carries out interval division based on zero crossing, then By thresholding method to signal c_WGNT () carries out white Gaussian noise filtering.It is assumed that comprise j zero crossing z in signal_i(i= 1,2 ..., j), then for interval z_i＜ t_j＜ z_i+1Interior signal has:

$s_{WGN}\left(t_j\right)=\left\{\begin{array}{cc}{c}_{WGN}\left(t_j\right)\frac{\left|c_{WGN}\left(t_m\right)\right|-T}{\left|c_{WGN}\left(t_m\right)\right|},& \left|c_{WGN}\left(t_m\right)\right|>T\\ 0,& \left|c_{WGN}\left(t_m\right)\right|\&le;T\end{array}\right.---\left(6\right)$

Wherein, T is the threshold value of signal c (t), t_mFor the extreme point in interval.Threshold value T is determined by formula (7) and (8):

$T=\mathrm{\&sigma;}\sqrt{2\ln L}---\left(7\right)$

σ=median (| c_WGN(t) |: t=1,2 ..., L)/0.6745 (8)

Herein, σ is input signal c_WGNThe noise level of (t)；L is signal c_WGNT the length of (), in particular for discrete letter Number L is discrete point and counts；Median is signal median operator, is used for obtaining seasonal effect in time series median.

Step 3) wavelet transform (DWT):

Wavelet transformation is a kind of new analysis method formed that time-domain analysis and frequency-domain analysis combined, and reflection is table The change that the facial muscle signal of telecommunication was presented on time and two dimensions of frequency, therefore the method is for both the above method Should have certain advantage in theory, it is possible to make full use of the information that surface electromyogram signal is comprised.Wavelet analysis is a kind of window The Time-Frequency Localization signal analysis method that size is fixed, shape is variable of mouth, i.e. has higher frequency discrimination in low frequency part Rate and relatively low temporal resolution, have higher temporal resolution and relatively low frequency resolution at HFS.

The signal S of a given a length of N, discrete wavelet transformation (DWT) at most can become log signal decomposition₂N number of frequency Rate level.The first step disassembles and starts from signal S, and after decomposition, decomposition coefficient is made up of two parts: low frequency coefficient vector cA1 and high frequency system Number vector cD1, the former also referred to as approximation (Approximation) composition, the latter is also referred to as details (Detail) composition.To Amount cA1 is obtained through convolution algorithm with low pass resolution filter by signal S, and vector cD1 is to be decomposed filter by signal S with high pass Ripple device obtains through convolution algorithm.In next step decomposes, by same method, low frequency coefficient cA1 is divided into two parts, i.e. The cA1 of signal S above is replaced, after decomposition, returns the low frequency coefficient cA2 and high frequency coefficient cD2 of yardstick 2, simultaneously yardstick 1 High frequency coefficient cD1 keeps constant, and the rest may be inferred continues to decompose.Finally obtain a low frequency coefficient row cA5 and five high frequency coefficients Row cD1, cD2, cD3, cD4, cD5, as shown in Figure 2.

The component that decomposition based on the EMG signal to lower extremity movement obtains, we are carrying out bandpass filtering treatment to it, are making CD1, the frequency band range of cD2, cD3, cD4, cA5 be respectively 256-512Hz, 128-256Hz, 64-128Hz, 32-64Hz and 16-32Hz。

Step 4) singular value decomposition (SVD):

According to the method in step 3, the EMG signal of lower extremity movement can be decomposed by DWT method, from step 3 The low frequency coefficient row cA5 obtained and five high frequency coefficient row cD1, in cD2, cD3, cD4, cD5, chooses 5 data and carrys out structure Become matrix expression, as selected: cD2, cD3, cD4, cD5 and cA5 constitute matrix, and the expression matrix obtaining 5 × I is as follows, M_ij= [d_2i；d_3i；d_4i；d_5i；a_5i], i=1,2 ..., I, wherein I is lower extremity movement sample action number, d2i, d3i, d4i, d5i and A5i correspondence DWT decomposes the coefficient obtained, and in order to obtain the eigenvalue that can be able to use in svm classifier, will use SVD here To M_ijCarry out decomposing and obtain singular value.

For M, N ∈ C^m×n, make equation U if there is m rank unitary matrice U and n rank unitary matrice V^TMV=N sets up, then it is assumed that Matrix M and N Unitary Equivalent.IfThen matrix M^TThe eigenvalue of M should meet following relation:

λ₁≥λ₂≥…≥λ_r≥λ_r+1=...=λ_n=0 (9)

Then can obtain, the singular value of matrix M is

IfThen certainly exist m rank unitary matrice U and n rank unitary matrice V so that matrix M, matrix U and Matrix V meets formula (10):

$U^{T}MV=N=\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]---\left(10\right)$

Σ=diag (σ₁, σ₂... σ_r), wherein σ_i(i=1,2 ..., r is the non-zero singular value of matrix M；Formula (10) It is equivalent to formula (11):

$M={\mathrm{UNV}}^T=U\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]V^T---\left(11\right)$

Here,

U=MM^T (12)

V=M^TM (13)

Wherein, formula (11) is also referred to as the singular value decomposition (SVD) of matrix M.

Step 5) svm classifier:

In order to lower extremity movement is classified, decomposed these matrixes M obtained by DWF_ijSingular value, be used as lower limb The time and frequency domain characteristics matrix of EMG signal.For chosen data segment, matrix M_ijDetermined by below equation:

Wherein, i is the singular value number of the data of experiment (Trial) every time, and j is EMG number of samples.

Vapnik et al. creates support vector machine (SVM) method.The basic thought of SVM is to map the data into higher-dimension sky In between, and find out the hyperplane with maximum decision-making edge.The principle of minimization risk that SVM formula is characterized has been demonstrated to be better than Other traditional experimental principles of minimization risk.In order to obtain a good real-time grading device of robustness, in the present invention, Eigenvalue needed for svm classifier will obtain in terms of time domain or frequency domain two, or the analysis method that use time-frequency domain combines.

Assume to exist the training vector x of l sample_iWith corresponding tag along sort y_i, (x₁, y₁) ..., (x₁, y₁)∈R^D ×-1,1}, then SVM is substantially to solve for quadratic programming problem:

$\min \frac{1}{2}\left|\right|w|{|}^2+C\underset{i=1}{\overset{1}{\mathrm{\&Sigma;}}}{\mathrm{\&zeta;}}_i$

y_i(w·x_i+b)≥1-ζ_i, ζ_i＞ 0 (i=1,2 ..., 1) (15)

Wherein, C is constant, ζ_iBeing slack variable, w is weight vectors, and b is deviation value.For example x, its discriminant function As follows:

$f\left(x\right)={\mathrm{\&Sigma;}}_{i=1}^{N_s}{y}_i{\mathrm{\&alpha;}}_{i}G(x,x_i)+b---\left(16\right)$

Wherein, N_sIt is the number supporting vector, α_iIt is positive Lagrange multiplier, G (x, x_i) it is kernel function.In SVM the heaviest The parameter of regularity wanted also is determined by the ten of training set times of cross-validation process.

This patent uses linear kernel function K processing feature matrix M_ij, classification function f (x) of support vector machine can describe For:

$f\left(x\right)=\mathrm{sgn}\left({\mathrm{\&Sigma;}}_{j=1}^n{y}_j{\mathrm{\&alpha;}}_{j}K\right(\stackrel{\&RightArrow;}{x_j},\stackrel{\&RightArrow;}{x})+b)---\left(17\right)$

N is number of samples, x_jIt is jth sample, y_jBe jth sample SVM output, K be for data conversion linear Kernel function, α_jIt it is the Lagrange multiplier of primal-dual optimization problem.

Sample data first passes through normalized (zero-mean+unit norm deviation processing), then by ten times of intersections Verification technique assessment classification identification effect.The present invention can make the recognition accuracy of human body lower limbs limb motion action obtain effectively Improve.

Claims (2)

1. the discrimination method for human body lower limbs surface electromyogram signal, it is characterised in that comprise the following steps:

A. putting viscous for disposable electromyographic electrode to lower hind limb musculature epidermis, Real-time Collection correspondence muscle masses sting in active actions Surface EMG signal under Jiing；

B. the EMG signal gathered in step a is carried out pretreatment, it is thus achieved that eliminate the EMG signal after artefact signal；Described pretreatment Method includes Hz noise filtering, baseline drift filtering and white Gaussian noise filtering；

C. the EMG signal using discrete small wave converting method to obtain step b decomposes, and obtains low frequency coefficient vector after decomposition CA1 and high frequency coefficient vector cD1；Low frequency coefficient vector cA1 is decomposed by the method using wavelet transform, it is thus achieved that low frequency Coefficient vector cA2 and high frequency coefficient cD2；Continue the method reusing wavelet transform low frequency coefficient vector to be carried out point Solve, until till obtaining low frequency coefficient vector cA5 and 5 high frequency coefficient row cD1, cD2, cD3, cD4, cD5；

When d. using, the filtering method that combines of frequency domain, the Wavelet Component obtained in step c is carried out singular value decomposition, and will point The singular value constitutive characteristic matrix that solution obtains；

E. with the eigenmatrix of acquisition in step d as sample, use support vector machine that feature samples is trained, and generation Hold vector machine classifier for blind sample is carried out Classification and Identification.

A kind of discrimination method for human body lower limbs surface electromyogram signal the most according to claim 1, it is characterised in that step Described in rapid b, the concrete grammar of pretreatment includes:

B1. Hz noise filtering；Method particularly includes: using EMG signal c (t) that collects in step a as original input signal C_PLIT () carries out notch filter, obtaining filter result is a (t), defines Hz noise factor ε_PLI, then ε_PLIEquation below can be passed through 1 calculates:

Wherein, var is signal variance operator, is used for calculating seasonal effect in time series variance, by Hz noise factor ε_PLITo filtering knot Fruit is modified shown in equation below 2:

Wherein, s_PLIT () is the final result filtering Hz noise noise, formula 2 shows if the energy of Hz noise noise accounts for Ratio exceedes the 10% of original energy, then use the filter result of notch filter；

B2. baseline drift filtering；Method particularly includes: the signal s that step b1 is obtained_PLIT () is as original input signal c_BW(t) Carrying out low-pass filtering, filter result is d (t), and the definition baseline drift factor is ε_BW, then ε_BWCan be calculated by equation below 3:

According to baseline drift factor ε_BWBaseline drift noise d (t) obtaining filtering is modified obtaining b (t), and its expression formula is such as Shown in formula 4:

Finally, the signal s after removing baseline drift noise is obtained_BWShown in (t) equation below 5:

s_BW(t)=c_BW(t)-b (t) (formula 5)；

B3. white Gaussian noise filtering；Method particularly includes: the signal s that step b2 is obtained_BWT () is as original input signal c_WGN T () is filtered；To time-domain signal c_WGNT () carries out interval division based on zero crossing, then by thresholding method to signal c_WGNT () carries out white Gaussian noise filtering；Particularly as follows: comprise j zero crossing z in hypothesis signal_i(i=1,2 ..., j), the most right In interval z_i＜ t_j＜ z_i+1Interior signal has:

Wherein, T is the threshold value of signal c (t), t_mFor the extreme point in interval；Threshold value T is determined by equation below 7 and formula 8:

σ=median (| c_WGN(t) |: t=1,2 ..., L)/0.6745 (formula 8)

Herein, σ is input signal c_WGNThe noise level of (t)；L is signal c_WGNT the length of (), in particular for discrete signal L It is discrete point to count；Median is signal median operator, is used for obtaining seasonal effect in time series median.