Skip to main content

Detection method of flexion relaxation phenomenon based on wavelets for patients with low back pain

Abstract

The flexion relaxation phenomenon (FRP) can be defined as a reduction or silence of myoelectric activity of the lumbar erector spinae muscle during full trunk flexion. It is typically absent in patients with chronic low back pain (LBP). Before any broad clinical utilization of this neuromuscular response can be made, effective, standardized, and accurate methods of identifying FRP limits are needed. However, this phenomenon is clearly more difficult to detect for LBP patients than for healthy patients. The main goal of this study is to develop an automated method based on wavelet transformation that would improve time point limits detection of surface electromyography signals of the FRP in case of LBP patients. Conventional visual identification and proposed automated methods of time point limits detection of relaxation phase were compared on experimental data using criteria of accuracy and repeatability based on physiological properties. The evaluation demonstrates that the use of wavelet transform (WT) yields better results than methods without wavelet decomposition. Furthermore, methods based on wavelet per packet transform are more effective than algorithms employing discrete WT. Compared to visual detection, in addition to demonstrating an obvious saving of time, the use of wavelet per packet transform improves the accuracy and repeatability in the detection of the FRP limits. These results clearly highlight the value of the proposed technique in identifying onset and offset of the flexion relaxation response in LBP subjects.

Introduction

Low back pain (LBP) is the most prevalent musculoskeletal condition; its economic burden is very large and appears to be growing [1]. Persons seeking care for LBP constitute a substantial proportion of patients seen in primary care offices [2]. Although it is difficult to predict who among those with first episodes of LBP will develop recurrent or chronic symptoms, there has been a recent surge in the study of both clinical and biological markers of “potential chronicity”. Several authors have tried to characterize typical neuromuscular and kinematic patterns that distinguish the LBP patients from healthy participants. Changes in feedforward EMG (Electromyogram) activation [3, 4], changes in erector spinae muscle fatigue characteristics [5], changes in the control of force production and movement production as well as changes in lumbopelvic kinematics have all been described [6] as potential indicators of physiological adaptation to chronic LBP. Amongst the various physiological indicators of LBP, the flexion relaxation phenomenon (FRP) has been one of the most studied surface electromyographic (sEMG) responses.

The FRP is defined by a reduction of silence of myoelectric activity of the lumbar erector spinae (LES) muscle during full trunk flexion [7]. This neuromuscular response is likely triggered by growing mechanical load in the ligaments and disks of the lumbar spine which are highly innervated by mechanoreceptors and nociceptors monitoring proprioceptive and noxious stimuli [8]. The onset and offset of the flexion relaxation response, as well as the amplitude and kinematics of the FRP, can be influenced by several factors such as loading of the trunk [9], lumbopelvic posture [9, 10], angular trunk velocity [11], task repetition [12], and muscular fatigue [13, 14]. Most importantly, persistent activation of the LES muscles during a full trunk flexion is typically observed in individuals with LBP as a protective “splinting” response to increase lumbar stabilization in response to pain or tissue injury [15, 16]. The absence of the FRP in LBP patients has also been associated with muscle spasm [17], reduced range of motion [18], and exaggerated stretch reflexes [18].

The sEMG FRP response is characterized, in healthy individuals, by three distinct phases that can clearly be delimited from one another; these limits are far less distinguishable in patients with LBP. FRP parameters are not only useful in fundamental research to study the underlying mechanisms of spinal stability, but their use in the clinical realm could also lead to significant development in the evaluation and clinical management of LBP. Effective, standardized, and accurate methods of identifying FRP limits are, however, needed before any broad clinical utilization of this neuromuscular response can be contemplated.

Determining the limits of the beginning and end of the FRP is performed manually. This manual process causes errors especially in the presence of strong non-stationary signals and presents a strong consumption of time. Therefore, automated algorithm methods of identifying these physiological events directly from sEMG signals are clearly needed before any widely used clinical applications of the FRP can be developed. Because muscular activities appear during the FRP relaxation phase in case of LBP subjects, the determination of the relaxation phase limits automatically is a non-trial problem for this kind of patients compared to healthy subjects. Indeed, in case of LBP subjects, the difference of amplitude between FRP phases is small and strong local impulses and artifacts appear during the relaxation phase. The contribution of this article consists in the development of an automated method based on wavelet transformation which would improve time limit detection of sEMG FRP in case of LBP people. Conventional Discrete Wavelet Transform (DWT) and Wavelet Packet Transform (WPT) [16] are compared with visual methods results obtained manually. Experimental system with real worlds data are used to validate our proposed approach based on observation and quantization of two physiological proprieties. The wavelet transforms (WTs) DWT and WPT permit to decompose a signal in several sub-signals of different frequencies ranges while maintaining a strong accuracy in the time domain [16]. Depending to the obtained wavelet decompositions which correspond to different groups of frequency, difference between FRP phases and artifacts can decrease. The idea of the proposed method of FRP limits detection consists of taking advantage of time domain sub-signals obtained with wavelets transforms from the collected sEMG signals. Instead of searching the limits of a single signal, we searched limits from all these sub-signals and realized an accurate detection based on all obtained limits.

Due to its ability to accurately analyze the time-frequency domain of non-stationary signals, WT appear in many works based on some bio-signals like sEMG signals [1625] and in non-destructive evaluation to determine the time of fly information in the ultrasonic signals [26]. Furthermore, the use of wavelets for the sEMG signals analysis is also legitimized by the fact that there is a similarity between the formation of a sEMG signal and the WT. The sEMG signal can be seen as a sum of little pulse waves named motor unit action potentials [17], and the WT is based on the fact that all signals are a sum of waves of different dilatations and translations [16]. WT is a powerful tool to cancel noises, as presented in [17] on sEMG signals. WT is also successfully used on sEMG signal in the following works: in [20] to automatically differentiate hand motions, in [21] to identify human emotions from face muscular activities, in [22] to separate deep muscles contribution during calf contraction, and in [23] to study muscles contraction. Furthermore, the study of time and frequency domains of sEMG signals is essential: the time domain gives information about the strength, and frequency domain analysis returns information on fatigue and on used muscle fibers [24, 27]. Moreover, in [25], the authors compare WT and Fourier transform for the estimation of fatigue. Authors of [18, 19] address LBP and use WT in their data processing. However, even if WT is successfully used, in both cases the analyzed experimental tasks are different from the FPR task and these articles do not address the search for the limits of the relaxation phase. In [18], wavelet decomposition is successfully applied on torso muscles sEMG signals to discriminate LBP and healthy subjects during a seating trunk rotation. In [19], a weight is attached to LBP patients and is suddenly dropped. WT helps to detect the drop on the spinal erector’s sEMG signal recorded during this task.

The following section describes the experimental set-up used to generate FRP sEMG signals. Background on WT and proposed methods are described in Section “Proposed method based on WT”. Experimental results of the proposed FRP analysis based on wavelet methods are shown in Section “Experimental results and analysis”. Finally, the conclusions are proposed in the last section.

Experimental system and sEMG signal

LBP participants were asked to perform a trunk flexion–extension task illustrated in Figure 1a. From an upright standing position with arms crossed over their chest, subjects were instructed to bend forward as far as possible during a 5-s movement period (flexion phase). Then, they were required to hold the fullyflexed position for a 3-s period and to return to the initial upright position. The extension phase lasted 5 s. Speed and duration of all the movement phases were standardized with a metronome. Each participant completed five flexion–extension cycles or trials, noted U = 5.

Figure 1
figure 1

Body motion of FRP task (a), FRP sEMG signal (b), and FRP kinematic angles (c) with limits t 1 and t 2 .

sEMG data were collected using bar bipolar active surface electrodes applied bilaterally over the LES at the level of lumbar vertebrae 5 (L5) approximately 3 cm from the mid-line (electrodes were applied in-line with muscle fiber direction). The two electrodes, placed on the left and right side of L5, simultaneously recorded muscle activities. The appearance of a sEMG FRP signal for the right sensor is shown in Figure 1b for a healthy subject. In this figure, all three phases of the back muscle activity can clearly be observed. Recall that the goal of the proposed algorithm is to automatically and accurately determine interval relaxation limits of time point t1 and t2 in Figure 1b. We noted error areas near each muscular activity limit in grey in Figure 1b. Those areas depict possible errors during limits detection. In our case, detection methods are applied to LBP subject and as shown in Figure 2 the limit determination is much difficult to realize. Indeed, in Figure 2, we can observe the high non-stationarity of FRP signal for an LBP subject compared to a healthy subject. However, it is essential to determine LBP FRP limits with good accuracy and good repeatability, in order to make inter-individual and intra-individual studies of FRP more reliable.

Figure 2
figure 2

FRP sEMG signal obtained for a healthy patient (a) and anLBP patient (b).

Electrode material was 99.9% Ag and interelectrode distance was fixed at 10 mm. A ground electrode was placed on the left anterior superior iliac spine (ASIS) of each participant. Skin impedance was reduced by (i) shaving body hair, (ii) gently abrading the skin with fine-grade sandpaper (Red Dot Trace Prep, 3 M, St. Paul, MN, USA), and (iii) wiping the skin with alcohol swabs. sEMG activity was recorded using a single differential Delsys Surface sEMG sensor with a common mode rejection ratio of 92 dB at 60 Hz, a noise level of 1.2 μV, a gain of 10 V/V ±1%, an input impedance of 1015 Ω (Model DE2.1, Delsys Inc., Boston, MA, USA) and sampled at F S  = 1000 Hz with a 12-bit A/D converter (PCI 6024E, National Instruments, Austin, TX, USA). The data were collected by LabView (National Instruments, Austin, TX, USA).

Additionally, kinematic data were collected by a motion analysis system (OptotrakCertus, Northern Digital, Waterloo, ON, Canada). Light-emitting diodes were positioned on the right side of each participant on eight landmarks: lateral malleolus, lateral condyle of the femur, greater trochanter, ASIS, S2, L1, T12, and C7. The data were sampled at 100 Hz and low-pass filtered by a dual-pass, fourth-order Butterworth filter with a cutoff frequency at 5 Hz. An example of resulting kinematic angles data is shown in Figure 1c synchronized with the recorded sEMG signal in Figure 1b.

Due to the high non-stationarity of signal in case of LBP subjects, there are no obvious limits and the evaluation of their detection methods is a non-trivial problem. By observing Figure 2, we can easily understand that visual detections could not be reliable enough to become a reference to evaluate automatics methods. However, these kinematic data will allow us to evaluate our limits detection methods. Indeed, two interesting physiological properties between muscles activities relaxation interval limits t1 and t2 and kinematics allow us to determine two assessment criteria. First, as shown in Figure 1c, limits t1 and t2 obtained from this sEMG signal can be reported on the corresponding kinematic signal and convert into angles φ1 and φ2. The physiological properties are as follow:

  1. i.

    If t 1 and t 2 are well determined, angles φ 1 and φ 2 for left and right sensors of a given trial have to be equivalent. This property represents our « accuracy criterion » described in expression (17) of Section “Proposed evaluation method”;

  2. ii.

    If t 1 and t 2 are well determined, angles φ 1 and φ 2 for a given subject k and a given sensor have to be equivalent for all trials. This property constitutes our « repeatability criterion», defined in Section “Proposed evaluation method” by the expression (18).

We evaluated the quality of our proposed methods using these properties in comparison with the conventional approach based on visual detection of relaxation interval limits. From Figures 1 and 2, we observe the similitude to voice activity detection method developed to determine the limit of voice and non-voice intervals in the signal. Some previous work is done using WPT to determine these limits [28].

Proposed method based on WT

General processing process

The general principle of the proposed limits detection method of the FRP is described in Figure 3. This diagram shows the following steps of the proposed process:

  1. 1.

    The first step called « sEMG signal acquisition » sampled at 1000 Hz and filtered the sEMG signal with a band-pass filter between 10 and 450 Hz in addition to dismissing the power supply harmonics of 60 and 100 Hz ones corresponding to the Certus kinematic sensor.

  2. 2.

    The second step consists in decomposing the filtered signal into signals using the DWT and the wavelet per packet transform (WPT) [16] described in Section “WT and wavelet per packet transform background” and applied on sEMG signal in Section “WTs application to sEMG signal”. So, instead of relying solely on the filtered signal, we base our limits detection on several time domain sub-signals of the same signal at different frequency levels. In this way, the accuracy and repeatability of the search limits is improved.

  3. 3.

    Before searching for limits, in the third step, each FRP sub-sEMG signal passes through a shaping algorithm presented in Section “Proposed shaping algorithm”. This algorithm has specifically been developed to determine the shape of the sEMG signal while eliminating residual forms mainly due to non-stationarity of the muscle activities.

  4. 4.

    The next step, depicted in Section “Limits detection of interval relaxation phase”, automatically determines the limits of main muscle activity on the left and right bound of relaxation phase with a simple threshold applied to each shape of each sub-signal from the sEMG signal. This way, we obtain as much estimated limits as the number sub-signals and a decision is made about searched point time limits, t 1 and t 2 in Figure 1.

  5. 5.

    We consider here one last step that evaluates the limits estimated based on movements captured by the body kinematic. The evaluation of obtained limits is not a trivial problem. Criteria of accuracy and repeatability have been defined and presented in Section “Proposed evaluation method”.

Figure 3
figure 3

General process of the proposed limits detection method with n = 1, 2, …, N, j = 1, 2, …, J, λ(j) = 1, 2, …, 2j, u = 1, 2, …, U, andk = 1, 2, …, K.

WT and wavelet per packet transform background

Used in many applications, WT and wavelet per packet transform (WPT) are signal processing methods which provide a set of decomposed signals in independent frequency bands. It returns accurate time-frequency information effective in case of non-stationary signals such as sEMG time series. Mathematically, the continuous wavelet transform (CWT) is defined as a convolution of a time series signal x(t) with a wavelet ψ(α,β)(t) of various scales α and translations β:

W α , β = x t ψ α , β * t d t ,
(1)

with

ψ α , β t = 1 α ψ t β α .
(2)

Expression (1) represents the continuous wavelet analysis equation. ψ(α,β)(t) is named « child wavelet » and is obtained from a basic function referred to the « mother wavelet» ψ(t) as shown in (2). This basic wavelet can take several forms with different properties [16]. The time-frequency analysis result W(α,β) is named the wavelet coefficients (WCs) at scale α and time β. These coefficients of the continuous wavelet decomposition accurately identify important time scale in a time series signal. However, this method does not directly reconstruct the signal at specific timescales; for this reason, the process proposed in this article is based, not on WCs, but on their corresponding time series representation. The DWT produces these reconstructions from any WCs from discrete signal x n, for n = 1, 2, …, N. This is accomplished by determining the WC in the same manner as described in (1) and (2), but with discrete dyadic version of α and β. The analysis equation of DWT is written

W j , p = n = 0 N x n ψ j , p * n p ,
(3)

with

ψ ( j , p ) n = 1 2 j ψ n 2 j .
(4)

Compared to the CWT α = 2j and β = 2jp. j = 1, 2, …, J are called the discrete levels of timescales and p = 0, 1, …, P are the discrete locations P j = N 2 j 2 j where · is the function of round down to the next integer. The DWT converts the real time series signal x[n] sampled at discrete time instants n into a real two dimensional space W(j,p), at each discrete level j and discrete location p.

Each of these WCs of the DWT can be thought of as « detail WCs » and can be used in the reconstruction of the time series signal at the specified level. It results the following synthesis equation named the « detail » time series signal of level j:

d j n = p = 0 P j W j , p ψ j , p n p
(5)

In order to complete the signal representation, it is necessary to add the low frequencies Q(j,p) corresponding to the scales larger than 2j. A scaling function χ(j,p)[n] is defined in (7) based on the basic function χ(t) called « father wavelet » and used to computed « approximation » coefficients (6).

Q j , p = n = 0 N x n χ j , p * n p
(6)

with

χ ( j , p ) n = 1 2 j χ n 2 j .
(7)

Similar to « detail » time series in (5), the « approximation » time series signal of the j th level can be obtained with the synthesis equation (8).

a j n = p = 0 P j Q j , p ϕ j , p n p
(8)

As generalized by Mallat [16], the DWT can be seen as filters bank structure described in Figure 4. Each stage of the structure represents a wavelet level j composed by two quadrature mirror filters h n and g n. h n is a low-pass filter and g n is a pass-band filter. These two filters are, respectively, obtained by inner product of the scaling function χ(t) and the wavelet function ψ(t) [16]. Considering the wavelet level j, a(j)n is computed by passing the previous level approximation components a(j−1)n in the filter h n and d(j)n by passing the same previous level approximation components a(j−1)n in the filter g n. a(j)n is called low-frequency approximation signal at the resolution j and d(j)n high-frequency detail signal. Each level detail signal d(j)n of the wavelet transformation corresponds to a representation of the signal at a given frequencies range [16]. Figure 5 gives an example of FRP sEMG signal decomposed with four levels of DWT. It is interesting to observe that for all decomposed levels time point t1 and t2 are preserved.

Figure 4
figure 4

Tree structure of DWT for J= 4 levels.

Figure 5
figure 5

sEMG signal decomposition with J= 4 levels of DWT.

The WPT method is an expansion of the described DWT decomposition, which presents more possibilities for signal processing information extractions. In Figure 6, at each level j, instead of only considering the decomposition of the (j – 1)th approximation signal a(j−1)[n], the (j – 1)th detail signal d(j−1)[n] is also decomposed. As depicted in Figure 6, and contrary to the case of the DWT, the wavelet per packet transform can be seen as a complete tree structure offering more resolutions and information about the input signal x[n]. To simplify the writing in Figure 6, synthesis signals of WPT are noted δ ( j , λ ( j ) ) [ n ] where λ(j) is the index of sub-signal of the level j where λ(j) = 1, 2, …, 2j. The wavelet per packet transform returns 2j representations of the input signal at each level. The DWT returns one detail signal at each level j and one detail signal and one approximation signal at the last level J. In this article, our goal is to take advantage of this large number of decompositions to determine the angles φ1 and φ2 in Figure 1.

Figure 6
figure 6

Tree structure of WPT for J= 3 levels.

WTs application to sEMG signal

During the FRP experimentation, two sEMG signals of back muscle activities are recording: one sEMG signals from a sensor positioned to left of L5 and a second one from a sensor positioned to right of L5. We define s(k,u)L[n] and s(k,u)R[n], the n th samples of the filtered sEMG signal coming respectively from the sensor left (index · L) and the sensor right (index · R) at the u th trial of the subject k. We note k= 1, 2, …, K, u = 1, 2, …, U and n = 1, 2, …, N, with K the number of tested subjects, U the number of repetition trials of the experimentation, and N the number of samples of the sEMG signal sampled at the frequency F S .

The proposed method globally depicted in Section “General processing process” is identically performed on both sEMG signals, s(k,u)L[n] and s(k,u)R[n]. Estimated t1 and t2 limits for these antagonist two FRP signals are used to evaluate all limits detection methods. In order to simplify the writing, the processing process is described for the FRP sEMG response of a given left-right sensor, noted s(k,u)[n].

Based on the explanation of the DWT and the WPT, the sEMG signal s(k,u)[n] is decomposed through J wavelet levels. Considering the WPT, we note s ( k , u , j , λ ( j ) ) [ n ] the λ(j)th time-series representation of the signal s(k,u)[n] at level j; with λ(j) = 1, 2, …, 2j and j = 1, 2, …, J. Considering J levels of the WPT method, the total number of obtained sub-signals N S is defined by N S = j = 0 J 2 j .

WPT is a generalization of the DWT, which is a particular path (left path) in the structure tree of Figure 6. Therefore, based on synthesis signals s ( k , u , j , λ ( j ) ) [ n ] of WPT, detail signals of the DWT are simply obtained by imposing λ(j) =2 for all j.

Proposed shaping algorithm

By nature, sEMG signal is a sum of electrical potential generated by muscle cells [17]. As a result, this signal and its decompositions are composed by random peaks of different amplitudes. Furthermore, some artifacts can appear during the acquisition. For these reasons, the detection of muscles activities in a relaxation phase, defined in time interval t1 and t2 (in Figure 1) could not perform well visually directly from s(k,u)n or s ( k , u , j , λ ( j ) ) [ n ] . It is necessary to develop an accurate shaping method in order to make the relaxation interval main muscular activities appear clearly. In order to efficiently realize this sEMG signal shaping, we proposed a method based on moving multi-size windows, which results are shown in Figure 7. Based on several moving windows of different sizes, the shaping method smooth the EMG signal by decreasing the contribution of local artefacts while preserving the accuracy of the smoothing. The proposed shaping process is an iterative method done from an initial point of the relaxation, noted ρinit (see Figure 7). ρinit is defined between the interval t1t2, preferably around the middle of the relaxation phase. Our multi-windows method is performed identically and independently of both sides of this initial starting point: (i) for time points ρ greater than ρinit, ρ ≥ ρinit and (ii) for time points ρ less than ρinit, ρ ≤ ρinit. Consequently, to simplify our notation, we only describe the process realized for ρ ≥ ρinit.

Figure 7
figure 7

sEMG signal with its shaping signal and its thresholds to detect limits t 1 and t 2 .

We consider ρ the processing time point, equal to ρinit at the first iteration, and V windows of index v = 1, 2, …, V with different sizes Θ(v). For ρ ≥ ρinit as depicted in Figure 8, all V windows are applied to the absolute value of sEMG sub-signal, s ( k , u , j , λ ( j ) ) [ n ] , and give V different sizes segments Z(v) of this sub-signal defined for all v as

Z v n = s k , u , j , λ j n for n = ρ , ρ + 1 , , ρ + Θ v 1 ,
(9)

where | · | is the absolute value function. The median value of all these segments is computed as M(v) = median(Z(v)[n]|n = ρ, ρ + 1, …, ρ + Θ(v) − 1) and v = 1, 2, …, V. The selected window referred as index ϖ corresponds to the lowest calculated median value, ϖ= argin f v Μ ( v ) | v = 1 , 2 , , V By choosing the window with the lowest median, we reduce the effect of local artefacts and ensure that only the major signal events are taken into account. It is a robust and simple method. The shaping signal of s ( k , u , j , λ ( j ) ) [ n ] , noted s ˜ ( k , u , j , λ ( j ) ) [ n ] , simply takes the value of the selected median M(ϖ) along the selected window interval Θ(ϖ):

s ˜ k , u , j , λ j n = Μ v = ϖ for n = ρ , ρ + 1 , , ρ + Θ v = ϖ 1.
(10)
Figure 8
figure 8

Process diagram of the iterative FRP shaping algorithm for ρρinit.

For the next iteration, the processing time point ρ becomes ρ = ρ + Θ(ϖ). And for ρ ≥ ρinit, the process continues and repeats itself until the end of signal, ρ + Θ(ϖ) − 1 = N(see Figure 8). As mentioned previously, the same iterative method is performed on the other side time of ρinit for ρ < ρinit. Once the shaping process performed, all shaping sub-signals s ˜ ( k , u , j , λ ( j ) ) [ n ] of all subject k and for their all U trails will be used to determine limits t1 and t2, as explained in the next section.

Limits detection of interval relaxation phase

The use of shaping signals simplifies the automated detection of time point t1 and t2. For each shaping signal s ˜ ( k , u , j , λ ( j ) ) [ n ] the time point detection process is similar for both side time of ρinit. For n < ρinit and n ≥ ρinit, the method is, respectively, looking for the time point t1 and for the time point t2. As illustrated in Figure 7, for both sides indexed 1 and 2 like for t1 and t2, distances between maximum and minimum amplitude values first are computed as follows

D 1 k , u , j , λ j = max s ˜ k , u , j , λ j n | n = ρ init 1 , ρ init 2 , , 1 min s ˜ k , u , j , λ j n | n = ρ init 1 , ρ init 2 , , 1
(11)
D 2 k , u , j , λ j = max s ˜ k , u , j , λ j n | n = ρ init , ρ init + 1 , , N min s ˜ k , u , j , λ j n | n = ρ init , ρ init + 1 , , N
(12)

Both side empirical threshold values, γ1 and γ2 (0 < γ ≤ 1) are imposed to maximum amplitude distances computed in (11) and (12). The time points t1 and t2 correspond to the position at which sub-signals amplitudes are equivalent to the percentage of distances γ 1 D 1 ( k , u , j , λ ( j ) ) and γ 2 D 2 ( k , u , j , λ ( j ) ) , as depicted in Figure 7. Therefore, time points t1 and t2 of the λ(j)th sEMG signal decomposition for the trial u of the subject k, noted t 1 ( k , u , j , λ ( j ) ) and t 2 ( k , u , j , λ ( j ) ) , are expressed as follow:

t 1 k , u , j , λ j = arg sup n s ˜ k , u , j , λ j n γ 1 D 1 k , u , j , λ j | n = ρ init 1 , ρ init 2 , , 1 Δ t ,
(13)
t 2 k , u , j , λ j = arg sup n s ˜ k , u , j , λ j n γ 2 D 2 k , u , j , λ j | n = ρ init , ρ init + 1 , , N Δ t .
(14)

In (13) and (14), Δt = 1/F S is employed in order to convert limits in time. At this point of the process, based on N S results of limits time point t1 and t2, we can make a decision to determinate time point interval of relaxation phase. We propose to apply two decision functions, one to compute the time point t ¯ 1 ( k , u ) by simply applying mean function, mean( · ):

t ¯ 1 k , u = mean t 1 k , u , j , λ j | j = 1 , 2 , , J and λ j = 1 , 2 , , 2 j ,
(15)

and a second to compute t ˜ 1 k , u by applying the median function median( · ):

t ˜ 1 k , u = median t 1 k , u , j , λ j | j = 1 , 2 , , J and λ j = 1 , 2 , , 2 j
(16)

Proposed evaluation method

As mentioned in Section “WTs application to sEMG signal”, we consider two sEMG sensors: one on the left side of L5 and another one on the right side of L5. Therefore, after the wavelet decomposition, the FRP shaping and the limits time point detection (see Figure 3), we determine limits time points t1(k,u)L(t1(k,u)R) and t2(k,u)L(t2(k,u)R) and for each trial u of each subject k for left (right) sEMG sensors.

In order to estimate limits time points t1 and t2, we report those estimations on kinematic signal, c(k,u)[n], which returns position angle in degree of the subject k during his u th trial. As shown in Figure 1, reporting limits time points t1 and t2 on c(k,u)[n], angles φ1 and φ2 are obtained, respectively. In the same manner, for each sEMG signal limits and based on limits time point t1(k,u)L(t1(k,u)R) and t2(k,u)L(t2(k,u)R) we determinate angles φ1(k,u)L(φ1(k,u)R) and φ2(k,u)L(φ2(k,u)R) for the left (right) sensor.

We are interested in the kinematics angles because the two physiological properties between muscles activities relaxation time point limits and kinematics angles permit us to determine two assessment criteria as follows:

  1. i.

    For each subject k, estimated kinematic angles must be similar between left and right sEMG sensors during the same trial u. The error between the angles of left and right sensors of the same trial u must tend to zero and permit to evaluate the accuracy criterion:

φ 1 k , u L φ 1 k , u R 0 , and φ 2 k , u L φ 2 k , u R 0
(17)

It is known that for a subject k during the flexion–relaxation task, the end of the flexion must occur at the same kinematic angle according to all U trials. The beginning of the extension must also appear at the same kinematic angle according to all trials. The standard deviation function, std( · ), of each angle φ1(k,u)L, φ1(k,u)R, φ2(k,u)L, and φ2(k,u)R in function of trials must tend to 0, with u = 1, 2, …, U, allows us to evaluate the repeatability criterion:

std φ 1 k , u L | u = 1 , 2 , , U 0 , std φ 1 k , u R | u = 1 , 2 , , U 0 , std φ 2 k , u L | u = 1 , 2 , , U 0 , and std φ 2 k , u R | u = 1 , 2 , , U 0.
(18)

Experimental results and analysis

Compared methods and parameters

Twenty-five participants, K = 25, with chronic LBP volunteered for the study: 14 men, 11 women with mean age of 28.3 years and standard deviation of ±8.9. All participants gave their written informed consent according to the research protocol approved by the local ethics committee. As previously depicted in Section “Experimental system and sEMG signal”, each subject completed five successive flexion and extension movements resulting (U = 5): five kinematics angles signals, five sEMG signals for the sensor localized to the left of L5 and five signals for the sEMG sensor to the right of L5.

In order to evaluate the interval limits of the relaxation phase, six different methods are compared based on the exact same sEMG experimental signals. Table 1 lists the main definitions of each studied method applying to left and right sEMG signals to find t1 and t2. Because it is the method currently used to analyze sEMG signal of FRP, the conventional method to determine the interval limit is based on the visual and manual techniques. This conventional method is compared to the proposed automated methods. Two experienced assessors performed the conventional method of FRP sEMG limits (t1 and t2). These two visual methods are referred as methods « visual1 » and « visual2 » in results. In order to observe the performances of different sizes of moving windows in the shaping method without WT, the method MonoWind is applied directly on the filtered sEMG signal based on the same process presented in Section “Proposed shaping algorithm” and Figure 8, but with only one moving window of size Θ(1) = 0.6 s with V = 1. The shaping method with V sizes of moving windows (Figure 8) is directly applied on the filtered sEMG signal (without WT) with Θ(1) = 0.6 s, Θ(2) = 0.8 s, Θ(3) = 1.6 s, and Θ(4) = 2 s (V = 4), and is referred as MultiWind method. Four different methods based on shaping with multi-moving windows and wavelet signal decomposition are also compared. The methods WPT/mean and WPT/med decompose the filtered sEMG signal with the wavelet per packet transform and, respectively, use a mean decision equation (15) and a median one Equation (16) to determine the two limits. Both based on DWT, DWT/mean and DWT/med, methods also estimate, respectively, limits with mean and median operations.

Table 1 Detection methods studied and their main properties

All parameters depicted in Table 2 have been chosen from parameters studies according to accuracy and repeatability criteria. Four levels of decomposition (J = 4, Table 2) represent the best compromise between obtained performances of criteria and processing time. Compared to other simulated wavelet family, the mother wavelet « Haar » shows good estimation of FRP limits. It results N S  = 4 decomposed signals in methods employing DWT and N S  = 30 sub-signals in WPT based methods. Following the same parameter selection, thresholds for limits detection (Section “Limits detection of interval relaxation phase”) equal to γ1 = γ2 = 0.2 have been chosen

Table 2 Parameters used for WT, shaping methods and limits detection

Statistical analyses

Accuracy and repeatability criteria for the limit t1 were analyzed using one-way repeated measures ANOVAs to determine the effect of time limit detection methods. Planned contrasts were used to test for a priori hypotheses. The MultiWind method was compared to both wavelet-type methods (based on WPT and DWT) in order to determine the additional contribution of wavelet-type methods to the detection. Planned contrasts were also used to compare the visual detection methods to the WPT wavelet-type methods. Finally, planned contrasts were used to assess the difference between the MonoWind and the MultiWind methods. Statistically significant results obtained for accuracy and repeatability criteria of t1 are specified in Figures 9a and 10a.

Figure 9
figure 9

Mean angles error between sensors right and left, as defined by the accuracy criterion (17), in function of limits detection methods for the limit t 1 (a) and limit t 2 (b) with corresponding standard error and statistically significant differences for t 1 .

Figure 10
figure 10

Mean for all subjects of standard deviation of angles through repeated trials, as defined by the repeatability criterion in (18), for sensors left and right in case of t 1 (a) and t 2 (b) in function of all methods with statistically significant differences for t 1 .

Accuracy and repeatability results

Figure 9 shows the ability of methods to find close boundaries between the left and right sensors for t1 and t2. This figure presents the means of accuracy criterion compared to all subjects and corresponding standard errors obtained for all methods. As explained in Section “Proposed evaluation method” with the criterion of accuracy (17), the mean angles error of right and left sensors for all trials and subjects must tend to 0.

First considering statistically significant differences for t1 in Figure 9a, methods based on WPT show an accuracy criterion significantly smaller than visual results (Visual1 and Visual2) and MonoWind method. It represents an interesting result for methods using WPT because no other method presents a significant difference with visual and MonoWind Methods. The use of WPT with the multi-windows shaping algorithm (methods WPT/med and WPT/mean) improves results compared to (i) using just the multi-windows shaping algorithm (method MultiWind) and (ii) using DWT with the multi-windows shaping algorithm (methods DWT/med and DWT/mean).

Based on average results, we also can make the following remarks:

  1. i.

    Both visual methods (Visual1 and Visual2) yielded equivalent mean error of approximately four degrees;

  2. ii.

    Unlike MultiWind method and all others, the MonoWind approach seemed less efficient for finding close limits between the right and left sensors. This technique shows the worse mean of the accuracy criterion and the highest variability. These bad results are mainly due to the presence of stronger local pulses and artifacts during the relaxation phase in case of LBP subjects. The shaping algorithm based on several moving windows is clearly more efficient for this kind of signals;

  3. iii.

    Even if there is no significant difference between those methods we observe that: (i) DWT-based methods offer better results than MultiWind and (ii) the mean error is the smallest for method based on WPT decomposition of sEMG signal;

  4. iv.

    In this case, no obvious difference appears between making a decision based on mean (15) or based on median (16).

Due to the nature of the flexion–extension task, looking for t2 (Figure 9b) is a smaller challenge than looking for t1. The extension phase uses more muscle activities and determining this limit from sEMG signal is easier. This fact explains the smaller mean angles errors obtained to find t2 compared to t1 for all methods. For this reason, no statistical analysis has been performed for the limit t2. Even if all methods return good right/left detection of t2, methods based on WPT presents the best average results (the smallest errors).

Figure 11 compares limits of all methods obtained with visual methods. If Figure 9 gives a good idea of the methods’ accuracy when they determinate limits from right and left sensors, this does not mean that limits are well localized. The goal of Figure 11 consists to inform us if estimated limits of methods are well localized. Our best limits localizations references are limits obtained by visual methods. Indeed, based on the confidence of experienced assessors, their results differ, but their limits localization definitively cannot be out of range. For each trial, subject, limit, and method, an error between visual estimates and all others methods has been computed and average results compared to subjects and trials are presented in Figure 11. Mean errors versus Visual2 and for limits t1 and t2 are exposed in this figure. Interesting remarks are given follow:

  1. i.

    First considering both visual methods (see paler gray color bar in each plot), the mean differences between the two visual methods are 4.7 degrees for t 1 Figure 11a and 5.4 degrees for t 2 Figure 11b. Based on the confidence of experienced assessors, we consider those results as errors references to analyze the others method. Smaller errors obtained with others methods will give a strong idea of their good limits localization;

  2. ii.

    Compared to Visual2 for t 1 and t 2 and except for MonoWind, automated methods present smaller differences with Visual2 than Visual2 versus Visual1.

Figure 11
figure 11

Mean angles error obtained for all methods compared to the Visual2 method for limit t 1 (a) and limit t 2 (b).

The angles differences of automated methods compared to visual methods are smaller than the differences between the two visual methods. This fact allows us to (i) be confident on automated methods limits localization and (ii) note the difference between visual estimations made by two different persons.

Based on the repeatability criterion defined in Section “Proposed evaluation method” by the expression (18), Figure 10 expresses the repeatability (in term of standard deviation) of limits estimation for t1 and t2 in function of all repetitive trials of a subject. Means of obtained standard deviations for all subjects is shown in Figure 10 for right and left sensors and for t1 and t2. As depicted by (18), results must tend to 0. According to significant differences indicated in Figure 10a for limits t1, we can make the following observations:

  1. i.

    Methods WPT/med and WPT/mean show a clear significant difference in term of repeatability compared to visual methods, MultiWind and MonoWind. The use of wavelet per packet transform clearly improves the repeatability criterion in the case of LBP subjects;

  2. ii.

    The multi-windows shaping algorithm employed in MultiWind provides a significant result of repeatability compared to the method MonoWind, which used a mono-window shaping technique;

  3. iii.

    MultiWind presents no significant difference with visual methods and with methods based on DWT considering the repeatability.

Based on average standard deviation, we can also note that

  1. i.

    MonoWind has the strongest difference results between trials limits estimation, especially in the case of right sensor. As explained in Figure 9, this shaping method is too sensitive to random artifacts present in sEMG signal;

  2. ii.

    Both visual methods (Visual1 and Visual2) have equivalent mean standard deviations in function of experimental trials repetitions of approximately 5 degree;

  3. iii.

    Once again, no meaningful difference is observed between taking a decision based on mean (15) or based on median (16).

In case of limits t2 (Figure 10b), all methods yield similar mean standard deviation in function of experimental trials repetitions. As already justified, it is easier to find this limit and all methods have a similar behavior.

Our analyses illustrate the difficulty to determine t1 compared to t2 in case of LBP subjects. In order to observe the distribution of obtained results in terms of accuracy and repeatability for t1, Figure 12 depicts box plot for each method with the smallest observation, the lower quartile, the median, upper observation, the largest observation and outliers. In the case of repeatability results (Figure 12b), we consider the mean between left and right sensors. Indeed, for all methods, a statistically study reveals no significant difference between left and right sensors. Significant differences presented in Figures 9a and 10a are based on data distributions exposed in Figure 12 for all methods.

Figure 12
figure 12

Box plot diagram of accuracy criterion (a) and repeatability criterion (b) in function of all methods showing the data distribution for limit t 1 .

The evaluation of accuracy and repeatability criteria demonstrates the interest of using wavelet decomposition to determine limits of the relaxation phase for LBP subjects. However, WPT methods which make their limits choices from 30 sub-signals show better performances than DWT based on 4 sub-signals. Methods with WPT yields to statically significant improvement compared to visual limits determination for both criteria. In addition, the proposed multi-windows shaping algorithm presents good abilities to smooth signals with strong local artifacts and impulses, as it is the case in sEMG response of LBP patients.

Computation times

It is also essential to mention the processing time to complete the detection of limits with all methods. In total, considering numbers of subjects, trials per subject, right and left sensors and 2 limits per signal to estimate, 500 limits were determined. It took in total 3 h for each experienced assessors and 8 min for all automated methods to estimate limits.

More precisely, we observe how many times each method takes to detect limits of all trials of one LBP subject. Methods MonoWind, MultiWind, respectively, realized the limits detection of one subject in 0.1 and 0.29 s. Wavelet-based methods need 20.8 s. Those methods are computed simultaneously with a Matlab wavelet function. It would be easy to reduce computing times of wavelet methods by optimizing it and compiling it in C language. Finally, the average processing time to visually determine limits of all trails of one subject is 7.2 min. We can clearly conclude that proposed methods offer a significant gain of time.

Conclusion

sEMG parameters based on timing characteristics within the sEMG signal and in relation to other biomechanical signals such as movement are often used in kinesiology and health-related research. From an experimental and technical standpoint, visual identification of discrete events embedded within sEMG signals is time consuming. Moreover, standardization between assessors and subjects has proven to be challenging and its evaluation is often neglected. On the other hand, several automated methods have been suggested to define the onset and offset of muscle activity. Among the proposed approaches, onset and offset detection based on calculation of the standard deviation range of the sEMG baseline before a certain activity is commonly used. However, the standard deviation-based approaches have several limitations. In addition, searching for limits of the relaxation phase of sEMG FRP response in case of LBP patients is a real challenge. For this kind of subjects, it is difficult to separate visually and automatically the three phases of FRP. WTs are employed here to improve the determination of FRP phases, by taking advantage of the time-domain sub-signals obtained with WT from the sEMG signal. The proposed method of limits detection based on wavelet per packet transform outperforms all others automated methods. Compared to visual identification, in addition to demonstrating an obvious saving of time, the use of WPT clearly improves the accuracy and repeatability in identifying onset and offset of the flexion relaxation response for LBP subjects. An optimization method is currently developed to determine automatically best parameters of the proposed method. Other future works should focus on the development and evaluation of FRP onset and offset detection in various sub-populations of LBP patients. Such tools should be integrated in kinesiology and health-related research and will eventually improve patients’ evaluation in various clinical spheres.

References

  1. Dagenais S, Caro J, Haldeman S: A systematic review of low back pain cost of illness studies in the United States and internationally. Spine J 2008, 8: 8-20. 10.1016/j.spinee.2007.10.005

    Article  Google Scholar 

  2. Dagenais S, Tricco AC, Haldeman S: Synthesis of recommendations for the assessment and management of low back pain from recent clinical practice guidelines. Spine J 2010, 10: 514-529. 10.1016/j.spinee.2010.03.032

    Article  Google Scholar 

  3. Hodges PW: Changes in motor planning of feedforward postural responses of the trunk muscles in low back pain. Exp Brain Res 2001, 141: 261-266. 10.1007/s002210100873

    Article  Google Scholar 

  4. Hodges PW, Richardson CA: Inefficient muscular stabilization of the lumbar spine associated with low back pain. A motor control evaluation of transversus abdominis. Spine 1996, 21: 2640-2650. 10.1097/00007632-199611150-00014

    Article  Google Scholar 

  5. Lariviere C, Arsenault AB, Gravel D, Gagnon D, Loisel P: Evaluation of measurement strategies to increase the reliability of EMG indices to assess back muscle fatigue and recovery. J Electromyogr Kinesiol 2002, 12: 91-102. 10.1016/S1050-6411(02)00011-1

    Article  Google Scholar 

  6. Esola MA, McClure PW, Fitzgerald GK, Siegler S: Analysis of lumbar spine and hip motion during forward bending in subjects with and without a history of low back pain. Spine (Phila Pa 1976) 1996, 21: 71-78. 10.1097/00007632-199601010-00017

    Article  Google Scholar 

  7. Colloca CJ, Hinrichs RN: The biomechanical and clinical significance of the lumbar erector spinae flexion-relaxation phenomenon: a review of literature. J Manipulative Physiol Ther 2005, 28: 623-631. 10.1016/j.jmpt.2005.08.005

    Article  Google Scholar 

  8. Holm S, Indahl A, Solomonow M: Sensorimotor control of the spine. J Electromyogr Kinesiol 2002, 12: 219-234. 10.1016/S1050-6411(02)00028-7

    Article  Google Scholar 

  9. Gupta A: Analyses of myo-electrical silence of erectors spinae. J Biomech 2001, 34: 491-496. 10.1016/S0021-9290(00)00213-X

    Article  Google Scholar 

  10. Callaghan JP, Dunk NM: Examination of the flexion relaxation phenomenon in erector spinae muscles during short duration slumped sitting. Clin Biomech (Bristol, Avon) 2002, 17: 353-360. 10.1016/S0268-0033(02)00023-2

    Article  Google Scholar 

  11. Sarti MA, Lison JF, Monfort M, Fuster MA: Response of the flexion-relaxation phenomenon relative to the lumbar motion to load and speed. Spine 2001, 26: 421-426. 10.1097/00007632-200109150-00019

    Article  Google Scholar 

  12. Olson MW, Li L, Solomonow M: Flexion-relaxation response to cyclic lumbar flexion. Clin Biomech (Bristol, Avon) 2004, 19: 769-776. 10.1016/j.clinbiomech.2004.05.007

    Article  Google Scholar 

  13. Descarreaux M, Lafond D, Jeffrey-Gauthier R, Centomo H, Cantin V: Changes in the flexion relaxation response induced by lumbar muscle fatigue. BMC Musculoskelet Disord 2008, 9: 1-9. 10.1186/1471-2474-9-1

    Article  Google Scholar 

  14. Descarreaux M, Lafond D, Cantin V: Changes in the flexion-relaxation response induced by hip extensor and erector spinae muscle fatigue. BMC Musculoskelet Disord 2010, 11: 1-7. 10.1186/1471-2474-11-1

    Article  Google Scholar 

  15. Lund JP, Donga R, Widmer CG, Stohler CS: The pain-adaptation model: a discussion of the relationship between chronic musculoskeletal pain and motor activity. Can J Physiol Pharmacol 1991, 69: 683-694. 10.1139/y91-102

    Article  Google Scholar 

  16. Mallat S: A Wavelet Tour of Signal Processing. Academic Press, San Diego, CA; 1998.

    Google Scholar 

  17. Tingting Y, Quan L, Qingsong A: Study on best wavelet packet based on independent threshold de-noising for MUAP. 2010 Second World Congress on Software Engineering (WCSE),Wuhan China 2010, 269-272.

    Chapter  Google Scholar 

  18. Kumar S, Prasad N: Torso muscle EMG profile differences between patients of back pain and control. J. Clin. Biomech. 2010, 25(2):103-109. 10.1016/j.clinbiomech.2009.10.013

    Article  MathSciNet  Google Scholar 

  19. Pope MH, Aleksiev A, Panagiotacopulos ND, Lee JS, Wilder DG, Friesen K, Stielau W, Goel VK: Evaluation of low back muscle surface EMG signals using wavelets. J. Clin. Biomech. 2000, 15(8):567-573. 10.1016/S0268-0033(00)00024-3

    Article  Google Scholar 

  20. Ibrahimy MI, Khalifa OO: Hand motion detection from EMG signals by using ANN based classifier for Human Computer Interaction. 2011 4th International Conference on Modeling, Simulation and Applied Optimization (ICMSAO), Lumpur, Malaysia 2011, 1-6.

    Chapter  Google Scholar 

  21. Murugappan M: Electromyogram signal based human emotion classification using KNN and LDA. 2011 IEEE International Conference on System Engineering and Technology (ICSET), Aura, Malaysia 2011, 106-110.

    Chapter  Google Scholar 

  22. Landry SC, Nigg BM, Tecante KE: Standing in unstable shoe increase postural sway and muscle activity of selected smaller extrinsic foot muscles. Gait Posture 2010, 32(2):215-219. 10.1016/j.gaitpost.2010.04.018

    Article  Google Scholar 

  23. Reaz MBI, Hussain MS, Mohd-Yasin F: EMG analysis using wavelet functions to determine muscle contraction. 8th International Conference on e-Health Networking, Applications and Services HEALTHCOM 2006 2006, 132-134.

    Google Scholar 

  24. Roy HS, Casavant D, Emley M, Gilmore LD, De Luca CJ: EMG analysis using wavelet functions to determine muscle contraction. Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society. vol 4, 1988, 1732-1733.

    Chapter  Google Scholar 

  25. Coorevits P, Danneels L, Cambier D, Ramon H, Druyts H, Karlsson JS, De Moor G, Vanderstraeten G: Correlations between short-time Fourier- and continuous wavelet transforms in the analysis of localized back and hip muscle fatigue during isometric contractions. J Electromyogr Kinesiol 2008, 18(4):637-644. 10.1016/j.jelekin.2007.01.006

    Article  Google Scholar 

  26. Legendre S, Massicotte D, Goyette J, Bose T: Wavelet-transform-based method of analysis for lamb-wave ultrasonic NDE signals. IEEE Trans. on Instrumentation and Measurement. 49(3) 2000, 524-530.

    Article  Google Scholar 

  27. Naït-Ali A: Advanced Biosignal Processing. Springer, Berlin; 2009.

    Book  Google Scholar 

  28. Chiodi R, Massicotte D: Voice activity detection based on wavelet packet transform in nonlinear space communication channel. Int. Conf. on Advances in Satellite and Space Communications (SPACOMM), Colmar, France 2009, 54-57.

    Google Scholar 

Download references

Acknowledgment

The authors thank Marie-Pierre Harvey (physical activity master degree student and DC), Emmanuelle Dion (chiropractic student) and Jean Daniel Dubois (psychology PhD student) for their contributions in data acquisition and statistical analyses. The authors wish to thank the Natural Sciences and Engineering Research Council of Canada and ReSMiQ for their financial support.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Daniel Massicotte.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Nougarou, F., Massicotte, D. & Descarreaux, M. Detection method of flexion relaxation phenomenon based on wavelets for patients with low back pain. EURASIP J. Adv. Signal Process. 2012, 151 (2012). https://doi.org/10.1186/1687-6180-2012-151

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1687-6180-2012-151

Keywords