CN107411741A - Multichannel myoelectricity Coupling Characteristics method based on coherence-Non-negative Matrix Factorization - Google Patents

Abstract

The invention discloses a multi-channel myoelectric coupling characteristic analysis method based on coherence-non-negative matrix decomposition, which combines traditional coherence analysis and non-negative matrix decomposition, first synchronously collects multi-channel myoelectric signals and preprocesses them, Then calculate the coherence of multi-channel EMG signals, and finally obtain the functional connection strength between muscles in each frequency band through non-negative matrix decomposition. The invention quantitatively analyzes the multi-channel myoelectric coupling characteristics and the intermuscular functional connection strength of each frequency band, provides an effective observation means for in-depth exploration of the movement control mechanism of the central nervous system, and has important application value in the field of rehabilitation medicine.

Description

Analysis method of multi-channel myoelectric coupling characteristics based on coherence-nonnegative matrix factorization

technical field

The invention relates to the research field of nervous system motion control mechanism, in particular to a multi-channel myoelectric coupling characteristic analysis method based on coherence-nonnegative matrix decomposition.

Background technique

Intermuscular coupling is the interrelationship and mutual coordination between different muscles during the movement of the limbs. By studying the coupling characteristics of each characteristic frequency band between multi-channel surface electromyography (sEMG), the functional connection between multi-channel muscles and the mechanism of execution and coordination of limb movements controlled by the central nervous system can be obtained. In recent years, based on the traditional coherence analysis method, research on the coupling characteristics between muscles in the process of limb movement has been carried out one after another. Some scholars use the consistency analysis method to calculate the normalization of the cross-spectral density of two EMG signals to the signal's self-spectral density function to reflect the coupling relationship of EMG signals in the frequency domain. However, the traditional intermuscular coherence analysis method can only reflect the coherence in the frequency domain, and cannot extract the characteristic information of EMG signals at different time-frequency scales, and cannot reflect the functional connection strength between muscles in each frequency band. In addition, during limb movement, multiple muscles move simultaneously, resulting in a single channel or two channel EMG signals that cannot fully reflect the functional coupling relationship between muscles during limb movement.

Contents of the invention

The purpose of the present invention is to provide a multi-channel myoelectric coupling characteristic analysis method based on coherence-nonnegative matrix decomposition, which can obtain the coupling characteristics between multi-channel electromyographic signals and can also reflect the functional connection strength between muscles in each frequency band.

In order to achieve the above object, the following technical solutions are adopted: the method of the present invention combines coherence analysis and non-negative matrix decomposition, and the specific steps are as follows:

Step 1, synchronously collect multi-channel EMG signals and preprocess them;

Step 2, calculating the coherence of the multi-channel EMG signal;

Step 3, obtain the intermuscular functional connection strength of each frequency band through non-negative matrix decomposition.

Further, in step 1, when collecting multi-channel EMG signals, the Trigno ^TM Wireless EMG acquisition device from Delsys Company of the United States was used, the resolution was set to 16bit, and the sampling rate was 2000Hz; The whole body is fixed on the bracket with a bandage to ensure that the posture remains unchanged during the experiment. Adjust the bracket so that the forearm is parallel to the ground, and the angle between the forearm and the upper arm is about 90 ^° . At the same time, the EMG signals of multiple muscles under the forearm pronation movement are collected .

Further, in step 1, when the electromyographic signal is preprocessed, an adaptive 50Hz power frequency notch filter is used to process the electromyographic signal to remove power frequency interference; select Butterworth third-order bandpass FIR filter to The electromyographic signal is processed so that the electromyographic signal is mainly concentrated between 5-200Hz.

Further, in step 3, the specific content of the analysis method for the inter-muscular coupling characteristics of each frequency band is as follows:

First, the coherence analysis of the multi-channel EMG signal is carried out, and then the coherence results between the multi-channel EMG signals are decomposed by using the non-negative matrix decomposition method to obtain the coherence of each frequency band, and then quantitatively analyze the coupling characteristics between the multi-channel muscles;

The coherence performance reflects the degree of correlation between two time series in the frequency domain. Let X and Y be two sets of time series. The formula for calculating the coherence of the two signals is as follows:

where S _XY (f) is the cross-spectral density of X and Y at frequency f, S _XX (f) and S _YY (f) are the autospectral densities of X and Y respectively; C _XY is the coherence of X and Y Its value range is 0-1; if C _XY (f) = 1, it means that X and Y are completely linearly related at frequency f; if C _XY (f) = 0, it means that X and Y are completely linearly related at frequency f Independent; if the value of C _XY (f) is between 0 and 1, it means that X and Y are partially linearly related at frequency f, and there may be a nonlinear relationship;

The basic idea of the non-negative matrix factorization method is: for any given non-negative matrix V _i×j , the non-negative matrix factorization method can find a non-negative matrix W _i×p and a non-negative matrix H _p×j such that

V≈WH (2)

or

In the formula, the matrix V _i×j is the connection matrix, the matrix W _i×p is the base matrix, and the matrix H _p×j is the coefficient matrix;

In non-negative matrix factorization, the objective function is used to measure the approximation of the decomposition result;

The rule for the iteration of the objective function of the Euclidean distance

The objective function ||V-WH|| ² is monotonous, but not an increasing function, and the condition for ||V-WH|| ² to remain unchanged is that the matrices W and H are fixed;

Rules for Iteration of Objective Functions for K-L Divergence

The objective function D(V||WH) is monotonous, but not an increasing function, and the condition for D(V||WH) to remain unchanged is that the matrices W and H are fixed;

The value of the cooperative number p, the number of EMG signal channels i, and the time series length j of the signal satisfy

(i+j)×p<i×j (8)

Since the matrix before and after decomposition only contains non-negative elements, the column vectors of the original matrix V can be interpreted as the weighted sum of all column vectors in the base matrix W, and the weight coefficients are the elements in the corresponding column vectors in H.

Compared with the prior art, the present invention has the following advantages: a new method for the research of multi-channel EMG coupling characteristics is proposed, which can quantitatively describe the functional coupling characteristics of EMG signals on different time-frequency scales, and then quantitatively analyze the multi-channel EMG coupling characteristics. The coupling characteristics between channel EMG signals and the functional connection strength between muscles in each frequency band are effective methods for analyzing the coupling characteristics between muscles in characteristic frequency bands. Research methods of multi-channel EMG signals.

Description of drawings

Fig. 1 is the flowchart of the method of the present invention.

Fig. 2 is a diagram of the myoelectric signal acquisition position of the method of the present invention.

Fig. 3 is a comparison chart before and after the electromyographic signal preprocessing.

Figure 4 is a diagram of the coherence of the multi-channel EMG signals of the subjects - non-negative matrix decomposition results.

detailed description

The present invention will be further described below in conjunction with accompanying drawing:

The EMG signal of the human body changes with the movement state of the limbs, and multiple muscles are involved in each movement at the same time. Therefore, it is particularly important to analyze the coupling characteristics of multi-channel EMG signals for the study of the motor control mechanism of the central nervous system. The present invention combines the traditional coherence analysis method with non-negative matrix decomposition, as shown in Figure 1, specifically including synchronous acquisition of multi-channel electromyographic signals, preprocessing of electromyographic signals, calculation of coherence of multi-channel electromyographic signals, and non-negative matrix decomposition , Analysis of intermuscular coupling characteristics of each frequency band, and evaluation of functional coupling characteristics. The specific method includes the electromyographic signal acquisition part and the signal processing part:

The electromyographic signal was collected by using Trigno ^TM Wireless EMG acquisition equipment from Delsys Company of the United States, the resolution was set to 16bit, and the sampling rate was 2000Hz. Before collecting the signal, the subject sat quietly with the upper arm drooping naturally, and the elbow was fixed on the bracket with a bandage to ensure that the posture remained unchanged during the experiment. The bracket was adjusted so that the forearm was parallel to the ground, and the angle between the forearm and the upper arm was about 90°. Wipe the surface of the skin under test with 75% alcohol before the electrode to remove the oil and dander on the skin surface, and collect the EMG signals of multiple muscles under the forearm pronation action synchronously, as shown in Figure 2. Flexor carpi digitorum superficiale (FDS), extensor carpi ulnaris (ECU), extensor digitorum (ED), extensor carpi radialis (ECR), flexor carpi radialis (FCR), palmaris longus (PL), brachialis biceps (BB) and brachioradialis (B).

The signal processing part includes preprocessing and analysis of the intermuscular coupling characteristics of each frequency band

Signal preprocessing: Myoelectric signal is a weak signal, which is easily disturbed by noise. The collected multi-channel original EMG signal needs to be preprocessed, and the EMG signal is processed by using an adaptive 50Hz power frequency notch filter. Remove power frequency interference; use Butterworth third-order band-pass FIR filter to process the electromyographic signal, so that the electromyographic signal is mainly concentrated between 5-200Hz. The comparison before and after signal preprocessing is shown in Figure 3. It can be seen from Figure 3 that the preprocessing effectively filters out the 50Hz power frequency interference and its double frequency interference in the original EMG signal.

Analysis of Muscle Coupling Characteristics in Each Frequency Band

The present invention firstly analyzes the coherence of the electromyographic signals of each channel, and then uses a non-negative matrix decomposition method to decompose the coherence results between muscles to obtain the coherence of each frequency band, and then quantitatively analyzes the coupling characteristics between muscles.

Coherence can reflect the degree of correlation between two time series in the frequency domain. Let X and Y be two sets of time series. The formula for calculating the coherence of the two signals is as follows:

In the formula, S _XY (f) is the cross-spectral density of X and Y at frequency f, and S _XX (f) and S _YY (f) are the auto-spectral densities of X and Y, respectively. C _XY is the coherence between X and Y, and its value range is 0-1; if C _XY (f)=1, it means that X and Y are completely linearly correlated at frequency f; if C _XY (f)=0, it means X and Y are completely independent at frequency f; if the value of C _XY (f) is between 0 and 1, it means that X and Y are partially linearly related at frequency f, and there may be a nonlinear relationship.

After the coherence analysis, non-negative matrix decomposition is performed on the coherence results, and then the intermuscular coupling characteristics of each frequency band of the multi-channel are analyzed.

The basic idea of NMF can be simply described as: For any given non-negative matrix V _i×j , NMF can find a non-negative matrix W _i×p and a non-negative matrix H _p×j such that

V _i×j ≈W _i×p H _p×j (2)

or

In the formula, the matrix V _i×j is the connection matrix, the matrix W _i×p is the base matrix, and the matrix H _p×j is the coefficient matrix.

In non-negative matrix factorization, the choice of objective function is very important, and it is used to measure the approximation of the decomposition result. There are many choices for the objective function of the non-negative matrix factorization algorithm, among which the objective function based on K-L divergence (Kullback-Leibler divergence) and Euclid distance (Eulidean distance) is more commonly used.

The objective function based on Euclidean distance is:

In the formula, the condition for the objective function ||V-WH|| ² to obtain the minimum value is V=WH, and the minimum value is 0.

The objective function based on K-L divergence is:

In the formula, the condition for the objective function D(V||WH) to obtain the minimum value is V=WH, and the minimum value is 0.

In the optimization problem, the matrices W and H are both variables. No matter which objective function is selected, the matrices W and H are not convex functions, so it is difficult to obtain the optimal solution. Therefore, the following iterative rules are adopted, which can not only ensure the operation speed, but also facilitate the operation.

Iteration rule for the objective function of Euclidean distance

The objective function ||V-WH|| ² is monotonous, but not an increasing function, and the condition for ||V-WH|| ² to remain unchanged is that the matrices W and H are fixed.

Iteration Rules for Objective Functions of K-L Divergence

The objective function D(V||WH) is monotonous, but not an increasing function, and the condition for D(V||WH) to remain unchanged is that the matrices W and H are fixed.

The value of the synergy number p is generally strictly selected by a suitable method, and it meets the number i of EMG signal channels and the time sequence length j of the signal.

(i+j)×p<i×j (10)

Since the matrix before and after decomposition only contains non-negative elements, the column vectors of the original matrix V can be interpreted as the weighted sum of all column vectors in the base matrix W, and the weight coefficients are the elements in the corresponding column vectors in H. This representation based on the combination of base vectors has a very intuitive semantic interpretation, which reflects the concept of "parts form the whole" in human thinking.

For verifying a kind of multi-channel myoelectric coupling characteristic analysis method based on coherence-non-negative matrix decomposition of the present invention, gather the upper limb sEMG of 6 healthy subjects (age is (25 ± 3) years old), test related The information is shown in Table 1.

Table 1 Relevant information of the subjects

The subjects were required to have no muscle fatigue before the experiment, be in a good mental state, and be familiar with the experimental procedure. According to the multi-channel myoelectric signal acquisition and processing process of the present invention, the upper limb surface myoelectric signals of healthy subjects are collected, and the coupling characteristics of the multi-channel myoelectricity are analyzed, and then the movement control mechanism of the central nervous system is studied.

Figure 4 is a diagram of the coherence-nonnegative matrix decomposition results of the multi-channel EMG signals of the subjects (different colors represent different coupling strengths). In Figure 4, the abscissa on the left is the frequency of the EMG signal, the ordinate is the signal spectrum, and the abscissa and ordinate on the right are the 8 muscles of the upper limbs of the subject; it can be seen directly from Figure 4 (a) to (d) It was found that the sEMG signals of the upper limbs of the subjects were decomposed into 4 frequency bands, and the different colors of the squares represented the coupling strength between different muscles, reflecting the difference in the coupling strength between muscles. Therefore, through coherence-nonnegative matrix decomposition of multi-channel EMG signals, the strength of intermuscular coupling characteristics corresponding to EMG signals in different frequency bands can be obtained, which provides effective observation data for in-depth exploration of the motor control mechanism of the central nervous system.

The above-mentioned embodiments are only descriptions of preferred implementations of the present invention, and are not intended to limit the scope of the present invention. All such modifications and improvements should fall within the scope of protection defined by the claims of the present invention.

Claims (4)

1. one kind is based on the multichannel myoelectricity Coupling Characteristics method of coherence-Non-negative Matrix Factorization, it is characterised in that institute State method to be combined coherent analysis and Non-negative Matrix Factorization, comprise the following steps that：

Step 1, synchronous acquisition multichannel electromyographic signal and it is pre-processed；

Step 2, multichannel electromyographic signal coherence is calculated；

Step 3, function connects intensity between each frequency range flesh is obtained by Non-negative Matrix Factorization.

2. the multichannel myoelectricity Coupling Characteristics side according to claim 1 based on coherence-Non-negative Matrix Factorization Method, it is characterised in that：In step 1, when gathering multichannel electromyographic signal, U.S. Delsys company's Ts rigno is utilized^TMWireless EMG collecting devices, resolution ratio are set to 16bit, sample rate 2000Hz；Before gathering signal, subject is sat quietly and naturally droops upper arm, Ancon is fixed on support with bandage, ensure experimentation in posture it is constant, adjusting bracket makes forearm parallel to the ground, forearm with The angle of upper arm is about 90 °, while gathers the electromyographic signal of the polylith muscle under the preceding action of forearm rotation.

3. the multichannel myoelectricity Coupling Characteristics side according to claim 1 based on coherence-Non-negative Matrix Factorization Method, it is characterised in that：In step 1, when being pre-processed to electromyographic signal, using adaptive 50Hz notch filters wave filter to flesh Electric signal is handled, and removes Hz noise；Electromyographic signal is handled from the rank band logical FIR filter of Butterworth three, Electromyographic signal is set to be concentrated mainly between 5-200Hz.

4. the multichannel myoelectricity Coupling Characteristics side according to claim 1 based on coherence-Non-negative Matrix Factorization Method, it is characterised in that in step 3, the particular content of Coupling Characteristics method is as follows between each frequency range flesh：

Multichannel electromyographic signal is subjected to coherent analysis first, then believed multichannel myoelectricity using non-negative matrix factorization method Coherence's result between number is decomposed, and obtains the coherence of each frequency range, and then is coupled between quantitative analysis multichannel muscle special Property；

Coherence property embodies degree of correlation of two time serieses on frequency domain, if X and Y is two groups of time serieses, two signal phases Dryness calculation formula is as follows：

<mrow> <msubsup> <mi>C</mi> <mrow> <mi>X</mi> <mi>Y</mi> </mrow> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>S</mi> <mrow> <mi>X</mi> <mi>Y</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> <mrow> <msub> <mi>S</mi> <mrow> <mi>X</mi> <mi>X</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> <mo>*</mo> <msub> <mi>S</mi> <mrow> <mi>Y</mi> <mi>Y</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

In formula, S_XY(f) it is cross-spectral densities of the X and Y on frequency f, S_XX(f)、S_YY(f) be respectively X and Y auto spectral density；C_XYFor X and Y coherence, its span are 0-1；If C_XY(f)=1, illustrate that X and Y are fairly linear related on frequency f；If C_XY(f)=0, illustrate that X and Y are completely independent on frequency f；If C_XY(f) value illustrates X and Y on frequency f between 0 to 1 Partial linear is related, it is understood that there may be non-linear relation；

The basic thought of non-negative matrix factorization method is：For any given nonnegative matrix V_i×j, non-negative matrix factorization method energy Enough find a nonnegative matrix W_i×pWith a nonnegative matrix H_p×jSo that meet

V≈WH (2)

Or

<mrow> <msub> <mi>V</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>&amp;ap;</mo> <msub> <mrow> <mo>(</mo> <mi>W</mi> <mi>H</mi> <mo>)</mo> </mrow> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>r</mi> </munderover> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>k</mi> </mrow> </msub> <msub> <mi>H</mi> <mrow> <mi>k</mi> <mi>j</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>

In formula, matrix V_i×jFor connection matrix, matrix W_i×pFor basic matrix, matrix H_p×jFor coefficient matrix；

In Non-negative Matrix Factorization, object function is used for weighing the approximation ratio of decomposition result；

For the rule of the iteration of the object function of Euclidean distance

<mrow> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>p</mi> </mrow> </msub> <mo>&amp;LeftArrow;</mo> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>p</mi> </mrow> </msub> <mfrac> <msub> <mrow> <mo>(</mo> <msup> <mi>VH</mi> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mrow> <mi>i</mi> <mi>p</mi> </mrow> </msub> <msub> <mrow> <mo>(</mo> <msup> <mi>WHH</mi> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mrow> <mi>i</mi> <mi>p</mi> </mrow> </msub> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow> 1

<mrow> <msub> <mi>H</mi> <mrow> <mi>p</mi> <mi>j</mi> </mrow> </msub> <mo>&amp;LeftArrow;</mo> <msub> <mi>H</mi> <mrow> <mi>p</mi> <mi>j</mi> </mrow> </msub> <mfrac> <msub> <mrow> <mo>(</mo> <msup> <mi>W</mi> <mi>T</mi> </msup> <mi>V</mi> <mo>)</mo> </mrow> <mrow> <mi>p</mi> <mi>j</mi> </mrow> </msub> <msub> <mrow> <mo>(</mo> <msup> <mi>W</mi> <mi>T</mi> </msup> <mi>W</mi> <mi>H</mi> <mo>)</mo> </mrow> <mrow> <mi>p</mi> <mi>j</mi> </mrow> </msub> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>

Object function | | V-WH | |²It is dull, but is not increasing function, and | | V-WH | |²It is matrix W and H to keep constant condition It is fixed；

For the rule of the iteration of the object function of K-L divergences

<mrow> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>p</mi> </mrow> </msub> <mo>&amp;LeftArrow;</mo> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>p</mi> </mrow> </msub> <mfrac> <mrow> <munder> <mo>&amp;Sigma;</mo> <mi>j</mi> </munder> <msub> <mi>H</mi> <mrow> <mi>p</mi> <mi>j</mi> </mrow> </msub> <msub> <mi>V</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>/</mo> <msub> <mrow> <mo>(</mo> <mi>W</mi> <mi>H</mi> <mo>)</mo> </mrow> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> </mrow> <mrow> <munder> <mo>&amp;Sigma;</mo> <mi>j</mi> </munder> <msub> <mi>H</mi> <mrow> <mi>p</mi> <mi>j</mi> </mrow> </msub> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mi>H</mi> <mrow> <mi>p</mi> <mi>j</mi> </mrow> </msub> <mo>&amp;LeftArrow;</mo> <msub> <mi>H</mi> <mrow> <mi>p</mi> <mi>j</mi> </mrow> </msub> <mfrac> <mrow> <munder> <mo>&amp;Sigma;</mo> <mi>i</mi> </munder> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>p</mi> </mrow> </msub> <msub> <mi>V</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>/</mo> <msub> <mrow> <mo>(</mo> <mi>W</mi> <mi>H</mi> <mo>)</mo> </mrow> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> </mrow> <mrow> <munder> <mo>&amp;Sigma;</mo> <mi>i</mi> </munder> <msub> <mi>W</mi> <mrow> <mi>i</mi> <mi>p</mi> </mrow> </msub> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>

Object function D (V | | WH) is dull, but is not increasing function, and it is matrix W and H that D (V | | WH), which keeps constant condition, It is fixed；

The value and electromyographic signal channel number i, the length of time series j of signal for cooperateing with number p meet

(i+j) × p ＜ i × j (8)

Non- negative element is only included due to decomposing in front and rear matrix, therefore, original matrix V column vector can be construed to basic matrix The weighted sum of all column vectors in W, and weight coefficient is the element in H in corresponding column vector.