CN106295212A - A kind of electromyographic signal synchronization processing method based on polynary empirical mode decomposition - Google Patents

Abstract

The invention discloses a kind of electromyographic signal synchronization processing method based on polynary empirical mode decomposition, the method synchronous acquisition polylith difference muscle group nervus motorius when limbs move encourages bioelectrical signals produced by muscle group, obtain a multichannel electromyographic signal (MEMG), then use polynary empirical mode decomposition (MEMD), the MEMG signal of collection is carried out synchronization process.The MEMG data collected are mapped in multi-C vector space by MEMD method by the Hammersley sequence of low difference, and then resolve into several intrinsic mode functions (IMFs) based on empirical mode decomposition is adaptive under the conditions of multivariate.At present, in MEMG process, presently disclosed method can effectively utilize body diverse location muscle group when being shunk by motion-activated produced bioelectrical signals time, spatial information.Compared with traditional method, the present invention can ensure that MEMG decompose after the concordance of interchannel IMFs exponent number, after the most also realizing decomposing, the corresponding each IMF information of interchannel has stronger dependency.

Description

A kind of electromyographic signal synchronization processing method based on polynary empirical mode decomposition

Technical field

The invention belongs to the technical field of multichannel (or multivariate) non-stationary nonlinear properties Synchronization Analysis, relate to one Multichannel electromyographic signal synchronization processing method based on polynary empirical mode decomposition.

Technical background

Electromyographic signal (Electromyography, EMG) is the physiology of human nerve musculature in muscle contraction The movable sign produced at body skin surface, is the nonlinear and nonstationary electro-physiological signals of a kind of complexity.Owing to SEMG is machine The direct forms of characterization of somatic nerves musculature, the physiological mechanism therefore excavating SEMG signal can disclose body to a certain extent The functional characteristic of neuromuscular system, can be that the physiological feedback mechanism of research human nerve musculature provides reason in a physiologically Opinion foundation；Clinical medical is alternatively Clinics and Practices neuromuscular system functional disease (such as motor neuron disease and The muscle diseases etc. such as amyotrophy) effectively help is provided；In rehabilitation engineering, the aspect such as including human motion rehabilitation, flesh energy evaluation Also directiveness suggestion it is provided that.

Existing SEMG Digital Signal Analysis and Processing method is primarily directed to signal under single channel and is analyzed, and mainly includes Five kinds of methods below.(1) wavelet analysis, the method is by wavelet basis function, by continuous differential so that be decomposed component of signal Time-frequency characteristics have preferable independence, process in real time, on SEMG signal, there is advantage, but owing to wavelet decomposition has intersection , can completely does not describes multicomponent, the natural quality of non-stationary SEMG signal, more can not reasonable dismissal motor neuron spy Levy rule.(2) high-order statistic, SEMG signal is modeled by the method according to linear time invariant model, can preferably assess flesh Meat constant force shrink (static) SEMG signal non-linear, non-stationary and non-Gaussian system degree, but shortcoming be the method can not be right SEMG characteristic under dynamic load or muscular force effect is analyzed.(3) artificial neural network, although the method can preferably be retouched State SEMG signal, but the training process that there is also certain limitation, such as artificial neural network generally requires substantial amounts of priori Data.(4) independent component analysis, the method is assumed that noise information and nervus motorius metamessage belong to different source signals, is passed through Signal mixed model is counter solves computing for structure, extracts signal mode, can process the signal of non-Gaussian system from sophisticated signal, but It is that the one-tenth merotype obtained after the method processes exists loss, often cannot Reverse reconstruction primary signal.(5) standard empirical pattern Decomposing and set empirical mode decomposition, above two method is all based on Hilbert-Huang transform, is different from traditional Fourier Linear and stable state frequency spectrum analysis method based on leaf transformation, these methods are special based on signal sequence time scale local The decomposition method of property, it is possible to direct, posterior and adaptive process SEMG signal.At present, standard empirical Mode Decomposition or Set ensemble empirical mode decomposition method is widely recognized in Analysis of nonlinear signals.

But, said method, owing to being analyzed merely with SEMG signal characteristic on time dimension, does not consider The spatial information of multi-channel surface myoelectric (Multi-channel Electromyography, MEMG) signal, is processing MEMG The information between the muscle group of body diverse location is have ignored during signal.At present, multiple physiological parameter monitoring and analysis and complexity The modeling of biological feedback system more can be relatively reasonable disclose body movement physiological reaction mechanism problem, the most progressively become forward position Study hotspot.Therefore, traditional SEMG signal analysis based on single channel is owing to lacking physiologic information on space scale Analyze, seriously constrain the development of above-mentioned technology.For this problem, this patent discloses one and divide based on polynary empirical mode The multichannel electromyographic signal synchronization processing method solved, the MEMG signal Hammersley sequence by low difference that the method will gather Map the data in multi-C vector space, and if then resolving into based on empirical mode decomposition is adaptive under the conditions of multivariate Dry intrinsic mode function (Intrinsic Mode Functions, IMF), effectively make use of comprised MEMG time space Information.The present invention can ensure that MEMG decompose after the concordance of interchannel IMFs exponent number, interchannel pair after the most also realizing decomposing Each IMF information should have stronger dependency.

Summary of the invention

Present invention is primarily aimed at and overcome the deficiencies in the prior art, it is provided that be a kind of many based on polynary empirical mode decomposition Passage electromyographic signal synchronization processing method, utilizes the time space information that MEMG signal is comprised, and is passed through by the MEMG signal of collection The Hammersley sequence of low difference maps the data in multi-C vector space, and then based on empirical mode under the conditions of multivariate Decompose and adaptive resolve into several natural mode of vibration components, effectively ensure that each passage IMFs exponent number concordance, secondly with Standard empirical Mode Decomposition or set ensemble empirical mode decomposition method are compared, the interchannel that polynary empirical mode decomposition (MEMD) obtains IMFs information has preferable dependency.

In order to realize object above, the technical scheme is that a kind of electromyographic signal based on polynary empirical mode decomposition Synchronization processing method, mainly comprises the steps that

Step 1: by p, p >=3 passage electromyographic signal collection instrument synchronizes pickup human body at special exercise operation condition Electromyographic signal at lower typical muscle group, obtains electromyographic signal sample data, constitutes one p and ties up primary signal s (t), and s (t)= s₁(t),s₂(t),…,s_p(t) }, wherein t is moment value；

Step 2: by the Hammersley sequence of low difference, is mapped to the s (t) acquired in step 1 p-1 and ties up spheroid The data point set u that upper formation is new_p={ u₁,u₂,…,u_p}；

Step 3: point set step 2 obtained, along direction vector collectionWherein, p-1 ties up spheroid Upper angle point set θ_v={ θ_v1,θ_v2,…,θ_p-1|θ_v∈R^p, direction vector v=1,2 ..., V carries out projecting and obtain one group respectively Set of projections

Step 4: the extreme value of the set of projections that calculation procedure 3 obtainsAnd the moment of correspondence

Step 5: the extreme value matrix that step 4 is obtainedInterpolation, obtains envelope curve collection

Step 6: envelope curve collection step 5 obtained is averaging, and obtains average envelope curve,

Step 7: by formula (5) extract component d (t):

$d\left(t\right)=s\left(t\right)-\stackrel{\‾}{a}\left(t\right)---\left(5\right)$

Now, if d (t) meets screening conditions, then d (t) is defined as natural mode of vibration component c (t), simultaneously by s (t)-d T () is considered as s (t), repeat above step 1 to 7；If d (t) is unsatisfactory for screening conditions, then d (t) is considered as s (t), repeats above step Rapid 1 to 7 obtains residual components r (t), stops cycle calculations；

After repeating screening by above step, primary signal s (t) will be broken down into N number of natural mode of vibration component, { c_N(t)}, N=1,2 ..., N and residual components r (t) sum, as shown in formula (6):

$s\left(t\right)={\mathrm{\Σ}}_1^N{c}_N\left(t\right)+r\left(t\right)---\left(6\right)$

Wherein, screening conditions areWherein a (t) is the amplitude of envelope averaged curve, and σ is regulation ginseng Number, its value determines and σ ∈ (0,1) according to practical situation.

Further, the Hammersley sequence obtaining low difference of described step 2 method particularly includes:

First the multichannel electromyographic signal data construct Halton sequence that will obtain, i.e. assumes that whole prime number collection is ascending After arrangement, front n prime number is { x₁,x₂,…,x_n, then the i-th sample of the Halton sequence on certain dimension is:

$r_i^x=\frac{a_0}{x}+\frac{a_1}{x^2}+\frac{a_2}{x^3}+...+\frac{a_s}{x^{s+1}}---\left(1\right)$

Wherein, radix x is integer, a_j(j=1,2 ..., s) it is the number on the j position after sample index i is launched by x system, s Representing that sample index i presses the total number after x system launches, therefore sample index i is as follows by the method for expressing of radix x:

I=a₀+a₁×x+a₂×x²+…+a_s×x^s (2)

Prime number is brought into formula (2) respectively, and the Halton sequence obtaining n sample is as follows:

$(r_i^{x_1},r_i^{x_2},...,r_i^{x_n})---\left(3\right)$

The Hammersley sequence being obtained n sample by Halton sequence is as follows:

$(\frac{i}{n},r_i^{x_1},r_i^{x_2},...,r_i^{x_{n-1}})---\left(4\right).$

Further, described step 3 method particularly includes:

New data point set u is formed on (p-1) dimension spheroid Hammersley sequential sampling obtained_p={ u₁,u₂,…,u_p} All direction vectors along (p-1) dimension spheroidProject to respectively on spheroid, obtain one group of set of projectionsNow, a length of unit length of set of projections, not there is physical significance, set of projections direction vector value represents u_pIn Each element.

Further, the method for described step 4 is:

By setting retrieval data length, and calculate zero crossing number, draw the extreme value that set of projections is correspondingAnd the moment

Further, described electromyographic signal of the present invention includes by getting involved and the myoelectricity telecommunications of non-intervention two ways pickup Number.

A kind of electromyographic signal synchronization processing method based on polynary empirical mode decomposition of the present invention, according to above-mentioned steps each The natural mode of vibration component that channel decomposition goes out contains the time space information of MEMG signal；Manifold based on polynary empirical mode decomposition Road electromyographic signal synchronization processing method, can ensure that MEMG decompose after the concordance of interchannel IMFs exponent number, secondly with standard empirical Mode Decomposition or set ensemble empirical mode decomposition method are compared, and the method is decomposed the corresponding each IMF information of the interchannel obtained and had relatively Strong dependency.

Accompanying drawing explanation

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

Fig. 2 lower limb leg experimental protocol；In Fig. 2, (a) is that sitting is stretched one's legs schematic diagram, and (b) is flexing schematic diagram of standing, (c) For walking schematic diagram；

Fig. 3 surface myoelectric electrode placement positions；In Fig. 3, (a) is human body front schematic view, and (b) is human body schematic rear view；

The sitting of Fig. 4 secondary sources normal group tested 1 stretch one's legs experiment MEMD result in vastus medialis part, Fig. 4 In (a) be 1-4 rank IMF component time-domain diagrams, (b) is 5-8 rank IMF component time-domain diagrams, when (c) is 9-12 rank IMF component Territory figure, (d) is 13-16 rank IMF component time-domain diagrams.

Detailed description of the invention

Describing in detail below in conjunction with the accompanying drawings, the present invention (compares at present when processing for multichannel electromyographic signal (MEMG) Set empirical mode decomposition (EEMD) method widely used in electromyographic signal analysis and process field) specific implementation method with And results of performance analysis.

Typical muscle group includes and is not limited to trapezius muscle, latissimus dorsi m., levator scapulae, rhomboideus, erector spinae, pectoralis major, abdomen Outer oblique, rectus abdominis m., pectoralis minor, intercostales externi, obliquus internus abdominis m., levator scapulae, triangular muscle, biceps brachii m., coracobrachialis, triangle Flesh, extensor pollicis brevis, extensor carpi radialis brevis, extensor carpi radialis longus, extensor digiti minimi, extensor digitorum, brachioradialis, palmaris longus, pronator ters, Flexor carpi ulnaris m., flexor carpi radialis, gluteus maximus, biceps femoris, gastrocnemius, extensor digitorum longus and musculus soleus；

The multichannel electromyographic signal synchronization processing method based on polynary empirical mode decomposition proposed by checking this patent Performance, first implement experiment and acquire 11 tested lower extremity movement SEMG data of normal male.In experimentation, Every subjects has been required three kinds of leg exercise actions, and including the leg extension under sitting state, (abbreviation " stretch by sitting Lower limb "), leg flexing (being called for short " flexing of standing ") under standing state and walking, as shown in Figure 2.In order to eliminate individual variation Impact, every tested at least perform 4 experiment actions.

The sample rate of multichannel electromyographic signal collection equipment is 1000Hz, sampling precision 14, unit mV.Equipment has 8 Individual digital channel and 4 analog channels.Experimentation have employed wherein 5 passages, including 4 digital channels for acquisition tables Face SEMG signal and 1 analog channel are for gathering the angle angle of thigh and shank internode.Initial data is stored directly in and sets In standby embedded memory card, and returned to the software equipment of data collecting system in real time by bluetooth.The SEMG of the 1st to 4 passage Signal is corresponding in turn to following muscle: rectus femoris (RF), biceps femoris (BF), vastus medialis (VM) and semitendinosus m. (ST).Myoelectricity electricity Pole is distributed as shown in Figure 3.

The four-way electromyographic signal data that will collect, carry out polynary empirical mode decomposition according to step 2-7 successively, point IMFs is obtained after solution.Sitting by tested 1 is stretched one's legs as a example by experiment, obtains 16 rank after the EMG Signal Decomposition Based at its vastus medialis IMFs, concrete outcome is as shown in Figure 4.Additionally, rectus femoris, biceps femoris and the EMG signal of three passages of semitendinosus m., the most same to time-division Solution obtains 16 rank IMFs.

In order to verify that this patent disclosed method is good and bad in the performance processing multichannel electromyographic signal, we have employed two Performance Evaluation index, i.e. intrinsic mode function (IMF) exponent number concordance and pattern extent queue, believe simultaneously with at present at myoelectricity Number analyze compared with set empirical mode decomposition (EEMD) method widely used with process field, above-mentioned decomposition result is entered Row assessment.Intrinsic mode function (IMF) exponent number concordance is to compare the diversity of the IMFs exponent number quantity that each channel decomposition goes out. And pattern extent queue can assess the degree of relevancy between the same order IMF component of different passage, for multichannel EMG signal This dependency the most each the highest passage correspondence IMF concordance is the best, and the degree of pattern queue is the highest.For above-mentioned IMFs, it is put down The computational methods of related property coefficient are as follows.First, respectively each IMF component is standardized and normalized IMF Power spectrum；Secondly, calculate the n-th 1 further (n1=1,2 ..., N) passage and the n-th 2 (n2=1,2 ..., N) m (m=of passage 1,2 ..., M) rank IMF component power spectrum between relative coefficient c (n1, n2), finally give a N rank relative coefficient square Battle array C_mC_m:

Work as n₁≠n₂n₁≠n₂Time calculate Matrix C_mC_mMeansigma methods ρ of element_mρ_m, this value is the m of this tested corresponding experiment The average correlation coefficient of rank IMF.Calculate the average correlation coefficient of each rank IMF successivelyCalculate the flat of this system number AverageThen can obtain this tested average correlation coefficient ρ of IMFs under corresponding experiment action.

Performance indications 1: intrinsic mode function (IMF) exponent number concordance

Table 1 tested four-way SMEG signal processes total exponent number of gained IMFs through EEMD and MEMD

As shown in Table 1, for the total exponent number of IMFs of different passages, EEMD is the most equal, and MEMD method gained IMFs Total exponent number of its different passages is equal.

Performance indications 2: pattern extent queue

According to decompose IMFs power spectrumanalysis, calculate sitting stretch one's legs, flexing of standing, walking three kinds of experiment actions average Relative coefficient, in order to weigh the pattern queue of IMFs.In order to contrast further, all tested average to normal group here Relative coefficient is added up.For one group of IMFs of certain decomposition of movement gained of single-subject, first calculate different interchannel The average correlation coefficient of same order IMFs, then reject the 80% energy respective frequencies average correlation system less than the IMFs of 20Hz Number, calculates the meansigma methods between residual coefficient.Finally ask for tested identical experiment action of normal group under identical decomposition method The average of coefficient of correspondence and standard deviation, result is as shown in table 2.

Table 2 normal group EEMD and the average correlation coefficient table of MEMD gained IMFs

As shown in Table 2, the average correlation coefficient similar of MEMD, for EMG signal tested in three kinds of experiment actions, Outside the average correlation coefficient of MEMD all more than 0.78, and in three kinds of experiment actions, the average correlation coefficient of EEMD is 0.69 Below.From tested angle analysis, in tested, flexing of the standing experiment of 81.82% of stretching one's legs in experiment for sitting 90.91% Tested and walking experiment in 90.91% tested, the average correlation coefficient of MEMD be higher than EEMD；Sitting is stretched one's legs in experiment In tested, the flexing of standing of 90.91% in the tested and walking experiment of 90.91% 63.64% tested, MEMD's is average relevant Property coefficient is higher than EEMD.In sum, except the average correlation coefficient similar of MEMD in walking experiment, MEMD's is average generally Relative coefficient is higher than EEMD.

Claims (5)

1. an electromyographic signal synchronization processing method based on polynary empirical mode decomposition, mainly comprises the steps that

Step 1: by p, p >=3 passage electromyographic signal collection instrument synchronizes pickup human body allusion quotation under special exercise operation condition Electromyographic signal at type muscle group, obtains electromyographic signal sample data, constitutes p dimension primary signal s (t), s (t)={ s₁ (t),s₂(t),…,s_p(t) }, wherein t is moment value；

Step 2: by the Hammersley sequence of low difference, is mapped to shape on a p-1 dimension spheroid by the s (t) acquired in step 1 The data point set u of Cheng Xin_p={ u₁,u₂,…,u_p}；

Step 3: point set step 2 obtained, along direction vector collectionWherein, angle on p-1 dimension spheroid Degree point set θ_v={ θ_v1,θ_v2,…,θ_p-1|θ_v∈R^p, direction vector v=1,2 ..., V carries out projecting and obtain one group of projection respectively Collection

Step 4: the extreme value of the set of projections that calculation procedure 3 obtainsAnd the moment of correspondence

Step 5: the extreme value matrix that step 4 is obtainedInterpolation, obtains envelope curve collection

Step 6: envelope curve collection step 5 obtained is averaging, and obtains average envelope curve,

Step 7: by formula (5) extract component d (t):

$d\left(t\right)=s\left(t\right)-\stackrel{\‾}{a}\left(t\right)---\left(5\right)$

Now, if d (t) meets stop condition, then d (t) is defined as natural mode of vibration component c (t), simultaneously by s (t)-d (t) It is considered as s (t), repeats above step 1 to 7；If d (t) is unsatisfactory for stop condition, then d (t) is considered as s (t), repeats above step 1 Obtain residual components r (t) to 7, stop cycle calculations；

After repeating screening by above step, primary signal s (t) will be broken down into N number of natural mode of vibration component, { c_N(t) }, N=1, 2 ..., N and residual components r (t) sum, as shown in formula (6):

$s\left(t\right)={\mathrm{\Σ}}_1^N{c}_N\left(t\right)+r\left(t\right)---\left(6\right)$

Wherein, screening stop condition isWhereinBeing the changing value of average envelope curve, a (t) is that envelope is average The amplitude of curve, and σ is regulation parameter, its value determines and σ ∈ (0,1) according to practical situation.

A kind of electromyographic signal synchronization processing method based on polynary empirical mode decomposition, its feature It is the Hammersley sequence obtaining low difference of described step 2 method particularly includes:

First the multichannel electromyographic signal data construct Halton sequence that will obtain, i.e. assumes the whole ascending arrangement of prime number collection After, front n prime number is { x₁,x₂,…,x_n, then the i-th sample of the Halton sequence on certain dimension is:

$r_i^x=\frac{a_0}{x}+\frac{a_1}{x^2}+\frac{a_2}{x^3}+...+\frac{a_s}{x^{s+1}}---\left(1\right)$

Wherein, radix x is integer, a_j(j=1,2 ..., it is s) number on the j position after sample index i is launched by x system, s represents Sample index i presses the total number after x system launches, and therefore sample index i is as follows by the method for expressing of radix x:

I=a₀+a₁×x+a₂×x²+…+a_s×x^s (2)

Prime number is brought into formula (2) respectively, and the Halton sequence obtaining n sample is as follows:

$(r_i^{x_1},r_i^{x_2},...,r_i^{x_n})---\left(3\right)$

The Hammersley sequence being obtained n sample by Halton sequence is as follows:

$(\frac{i}{n},r_i^{x_1},r_i^{x_2},...,r_i^{x_{n-1}})---\left(4\right).$

A kind of electromyographic signal synchronization processing method based on polynary empirical mode decomposition, its feature It is described step 3 method particularly includes:

New data point set u is formed on (p-1) dimension spheroid Hammersley sequential sampling obtained_p={ u₁,u₂,…,u_pAlong (p-1) all direction vectors of spheroid are tieed upProject to respectively on spheroid, obtain one group of set of projectionsNow, a length of unit length of set of projections, not there is physical significance, set of projections direction vector value represents u_p In each element.

A kind of electromyographic signal synchronization processing method based on polynary empirical mode decomposition, its feature The method being described step 4 is:

By setting retrieval data length, and calculate zero crossing number, draw the extreme value that set of projections is correspondingAnd the moment

A kind of electromyographic signal synchronization processing method based on polynary empirical mode decomposition, its feature It is that described electromyographic signal of the present invention includes by getting involved and the myoelectricity signal of telecommunication of non-intervention two ways pickup.