# A Novel Noise Suppression and Artifact Removal Method of Mechanomyography Based on RLS, IGWO-VMD, and CEEMDAN

Authors: Zebin Li; Lifu Gao; Wei Lu; Daqing Wang; Huibin Cao; Gang Zhang
Journal: Journal of Sensors
 (2022)

---

## Abstract

Mechanomyography (MMG) signals have extensive applications in muscle function assessment and human intention recognition. However, during signal acquisition, MMG signals are easily contaminated by noise and artifacts, which seriously affects the recognition of their characteristics. To address these issues, a novel noise suppression and artifact removal method based on recursive least square (RLS), improved Gray Wolf Optimizer-optimized variable mode decomposition (IGWO-VMD), and complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) is proposed. In this paper, the RLS algorithm is first applied to adaptively filter out the power line interference (PLI). Then, IGWO is designed to select the appropriate VMD parameters and use the VMD to decompose the noisy signal into band-limited intrinsic mode functions (BLIMFs). In addition, the BLIMFs are classified into the low-frequency part and high-frequency part according to the given correlation coefficient (CC) threshold value. The effective components of the low-frequency part are identified by the center frequency. Meanwhile, the high-frequency part is decomposed by CEEMDAN, and its effective components are obtained according to the proposed sample entropy threshold range. Finally, the effective components of the low and high-frequency parts are reconstructed to obtain the denoised signal to realize the extraction of useful signals. Simulation experiment results demonstrate that the proposed method outperforms the classical methods and the designed IGWO-VMD method in terms of denoising performance. The effectiveness of the proposed method is verified through the measured MMG signal experiments. The proposed method not only effectively suppresses noise and artifacts but also overcomes the limitations of VMD and CCEMDAN.


---

## Body

## 1. Introduction

In biomedical engineering, mechanomyography (MMG) is an important surface biosignal that is commonly used to record mechanical vibrations characteristics generated from gross lateral movement of the muscle fiber at the initiation of a contraction [1]. Since MMG signals contain rich information on muscle activity, it has received a lot of attention in recent years and has been successfully applied in many fields, such as human-machine interface, muscle fatigue assessment, rehabilitation exercise, and muscle disease diagnosis [2–5]. However, MMG signals are often disturbed by various noises during acquisition, and their useful components can be masked by strong background noise, such as power line interference (PLI) and noise generated by electronic devices. In addition, MMG signals acquired using acceleration sensors are disturbed by artifacts such as gravitational acceleration information, motion caused by limb movement, and tremors generated by muscle contractions. These seriously affect the quality of the MMG signal and further increase the difficulty of signal detection, feature extraction, and identification. Therefore, it is of great significance to effectively remove the noise and correct artifacts from the measured MMG signals for further MMG signal identification and application.Many methods are used to denoise and correct artifacts in various signals. The main common methods include Fourier transform-based methods, wavelet transform (WT), empirical mode decomposition (EMD), and variable mode decomposition (VMD). MMG signals are not only nonlinear, nonstationary, and non-Gaussian but also chaotic and fractal [5]. Krueger et al. [3] collated 43 MMG studies published from 1987 to 2013, showing that it appears to be a consensus on the use of Butterworth filters to obtain useful MMG signals. However, the use of Fourier transform-based band-pass filters can filter out noise as well as important high-frequency information in the signal and cannot analyze local frequencies; in addition, useful low-frequency information in MMG signal may also be lost, when removing low-frequency artifacts. The energy of MMG signals is mainly distributed in the low-frequency band, which can overlap with artifacts at low frequencies, making it difficult for the band-pass filter to distinguish them in the frequency domain. As a result, many Fourier transform-based denoising methods do not always perform well and have little effect in correcting artifacts when applied to nonlinear and nonstationary signals similar to MMG signals [6].In the 1980s, a new time-frequency interaction analysis method, namely, wavelet analysis which is widely used in time-series signals and image noise reduction [7–9], emerged, which decomposes the signal into different frequencies by multiscale analysis and is suitable for nonlinear signals with good time-frequency localization properties. However, wavelet analysis that is developed from Fourier transform is limited in removing noise because of problems in selecting the appropriate wavelet type, the number of wavelet decomposition layers, and the threshold selection. In addition, a certain degree of distortion may occur when reconstructing the signal [9]. To overcome the drawbacks of the conventional wavelet transform, a second-generation wavelet transform (SGWT), which does not depend on the Fourier transform, is proposed. SGWT consists of three phases: splitting, predicting, and updating, in which the wavelet coefficients generated in the decomposition process need to set a threshold for processing, and the processed wavelet coefficients are inverted to obtain the denoised signal. Usually, the threshold processing methods include the soft threshold method and the hard threshold method. Since the hard thresholding method would produce oscillations at discontinuity points, the signal reconstructed using it does not have the smoothness of the original signal. The signal processed using the soft thresholding method has better continuity, which is also a commonly used method, but it produces bias and thus reduces the approximation to the original signal. Therefore, in practice, it is necessary to improve the soft threshold function, such as introducing Savitzky-Golay smoothing algorithm [10] and semisoft thresholding [8], to reduce the deviation between the estimated wavelet coefficients and the original wavelet coefficients, making the reconstructed signal approximate the real signal. However, setting the appropriate thresholding function is still a nontrivial task.With the development of novel techniques, empirical mode decomposition (EMD) [11] and VMD [12] have been proposed for the analysis of nonlinear and nonstationary signals. Unlike wavelet analysis, EMD and VMD have significant adaptability by avoiding the dependence on basis functions. EMD overcomes the shortcomings of signal processing that require prior knowledge. However, due to its theoretical defects, the decomposition process is prone to problems such as mode aliasing and endpoint effects [13], and the denoising effect is not satisfactory. Subsequently, ensemble empirical mode decomposition (EEMD) [14] was introduced, which is essentially a multiple empirical mode decomposition with superimposed Gaussian white noise. The overall average of the corresponding intrinsic mode function (IMF) obtained by multiple EMDs is used to eliminate the added white noise and suppress mode aliasing. However, the reconstructed component still contains residual noise of a certain amplitude. Although the reconstruction error can be reduced by increasing the number of integrations, it increases the computational scale. To overcome this problem, Torres et al. [15] proposed complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) based on EEMD. It adds adaptive white noise at each stage of decomposition and obtains each mode component by calculating a unique residual. Compared with EEMD, the reconstruction error is almost zero regardless of the number of integration, and the decomposition process is complete, overcoming the problem of low efficiency of EEMD decomposition. Due to these important properties, CEEMDAN has been widely used in biomedical engineering [16], seismology [17], etc. However, effective IMFs require appropriate methods for identification and selection.In 2014, Dragomiretskiy and Zosso [12] proposed a nonrecursive mode decomposition method, namely, variational mode decomposition (VMD), according to the constrained variational problem, which has a solid theoretical foundation. VMD can decompose any signal into an ensemble of band-limited intrinsic mode functions (BLIMFs). It has significant advantages in processing nonrecursive signals and not only can overcome the mode aliasing problem in EMD but also can use its own Wiener filtering characteristics to obtain a better filtering effect. VMD has received a lot of attention from researchers and has been successfully applied to many fields such as mechanical diagnosis, biomedical science, and other signal processing [18–20]. Although VMD obtains excellent results in signal denoising, these effects are determined by two critical parameters, namely, the number of modes K and the value of the penalty factor α, which are usually selected within a certain range. Selecting these parameters by trial-and-error would require a large number of operations and would waste a lot of time. Therefore, these critical parameters are usually determined based on empirical methods, which greatly limits the performance of VMD and may lead to inaccurate decomposition results. Hence, appropriate methods are needed to obtain the optimal values of these parameters. A popular approach is to use an intelligent optimization algorithm to adaptively determine the combination of parameters. There are many intelligent optimization algorithms such as genetic algorithm [21], particle swarm algorithm [22], artificial fish swarm algorithm [23], and gray wolf algorithm [24]. For example, Zhang et al. [25] proposed a parameter-adaptive VMD method based on the grasshopper optimization algorithm (GOA) to analyze vibration signals from rotating machinery. Similarly, to address the limitations of traditional VMD methods, Ni et al. [26] proposed a fault information-guided variational mode decomposition (FIVMD) method to improve the VMD method for extracting weak bearing repetitive transients and obtained significant fault-related frequencies.Grey Wolf Optimizer (GWO) is a new swarm intelligence optimization algorithm proposed by Mirjalili et al. [27] in 2014. The algorithm simulates the predation behavior of gray wolves and seeks the global optimal solution by tracking, encircling, chasing, and attacking prey. The GWO algorithm has the advantages of fewer input parameters, fast solution speed, and high accuracy. Among these optimization algorithms, the GWO algorithm has stronger competitiveness [28]. Combined with the previous work [29], the improved GWO (IGWO) has a better global optimization-seeking performance. Therefore, in this paper, the improved GWO is used for VMD parameter optimization.VMD is one of the most popular techniques in the field of biomedical signal processing. However, Lahmiri and Boukadoum [30] found that it is not safe to use VMD alone for denoising. To solve the problems caused by using CEEMDAN and VMD alone, some researchers have used EMD-related methods in combination with VMD methods in recent years and achieved better results [31–34]. From these literatures, it can be seen that the combination of the two methods can significantly improve the decomposition process and has a better noise suppression effect than a single one.In the specific application of denoising methods, Maji and Pal [35] found some differences between EMD and VMD in the range of signal decomposition levels and in the ability to extract low and high frequencies from the signal; i.e., the VMD method can better process high-frequency signals, while EMD has better effects on low-frequency signals. Therefore, this paper proposes a novel denoising method based on the above-mentioned research work. Firstly, VMD is used to decompose the measured MMG signal to obtain the low and high-frequency parts of the signal, then, CEEMDAN is used to extract the effective components of this high-frequency part, and finally, the useful information in the measured MMG signal is obtained by reconstructing the effective components in the low and high-frequency parts. The main contributions of the proposed method are presented as follows:
(1)
To suppress the 50 Hz PIL, the RLS algorithm is first applied to the measured MMG signal(2)
To obtain the optimal combination of the VMD parametersK,a, an IGWO algorithm is proposed in which the energy entropy is selected as the fitness function(3)
To improve the VMD denoising performance, the CEEMDAN algorithm is used to decompose the high-frequency BLIMFs. In the CEEMDAN algorithm, a sample entropy threshold range is proposed to identify the effective components of the high-frequency part(4)
To evaluate the effectiveness of the proposed method on the denoised signals reconstructed from the low and high-frequency effective components, the power spectral density (PSD) analysis method is used in addition to the frequency spectrum analysis methodThe rest of this paper is arranged as follows. The related theories, including RLS, IGWO-VMD, and CEEMDAN, are introduced in Section2. The proposed method is presented in Section 3. In Section 4, the proposed method is compared with other methods through simulation experiments, and the superiority of the proposed method is proved. In Section 5, the proposed method is applied to the actual measured MMG signals. Finally, the conclusions are drawn in Section 6.


## 2. The Basic Methods

### 2.1. RLS Algorithm

PLI is an electromagnetic noise from electronic equipment and transmission lines during biosignal measurements, which reduces the quality of the signal [36]. Therefore, removing PLI from the original signal is a problem that cannot be ignored. PLI is usually described as additive noise, having a sinusoidal signal at a fixed frequency (50 Hz or 60 Hz) with an unknown phase and amplitude. Many conventional filtering techniques are used to solve this problem [37]. However, these filters will fail in the case of PLI with frequency drift. Ahmed et al. [38] compared LMS, NLMS, and RLS algorithms and found that RLS has the best performance in removing PLI. The recursive least square (RLS) method has been widely used because of its advantages of the easy numerical solution and fast parameter convergence [39]. The RLS algorithm minimizes the weighted linear least square cost function by calculating the filter coefficients recursively. Therefore, in this paper, RLS is used to remove the PLI from the measured MMG signals. The specific description of the RLS algorithm is as follows [38, 40]:
(1)ek=dk−yk,(2)yk=ωkxk,(3)Kk=Pk−1∙xkλ+xHk∙Pk−1∙xk,(4)Pk=λ−1Pk−1−λ−1KkxHkPk−1,(5)ωk+1=ωk+Kk∙ek,where ek is the error signal, dk is desired signal, yk is the output signal from the adaptive filter, xk is the filter input vector, ωk is the filter coefficients vector, Kk denotes the gain vector, and Pk is a correlation matrix of the input signal. In this paper, P0=I/c, where I is the unit matrix and c is set to 0.01; λ is the forgetting factor which is set to 0.99. In addition, filter order M is set to 2.For PLI in the measured MMG signals, it is difficult for a hardware to filter it out completely, but the adaptive filtering algorithm RLS cannot only suppress it well, but PLI-free MMG signals can pass through the filter unchanged.


### 2.2. IGWO-VMD Algorithm

VMD is a new time-frequency analysis method, which iteratively searches for the optimal solution of the variational modes, continuously updates each mode function and the center frequency, and obtains an ensemble of BLIMFs, avoiding the endpoint effect and spurious component problems during the iterative process. The constrained variational optimization problem using the alternate direction method of multipliers can be described as follows:(6)minφk,ωk∑i=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kut=st,where K represents the number of modes, k=1,2,⋯,K; uk and ωk are the modes and their center frequencies, respectively, uk= u1,⋯,uK, ωk=ω1,⋯,ωK; st is the VMD processed signal; δt is the Dirac distribution function; ∗ and ∂t denotes the convolution and partial differential operators, respectively; j is the imaginary unit, and ∙2 represents L2 norm.To solve the above constrained variational model, a quadratic penalty factorα and Lagrangian multiplier λt are introduced to transform the constrained variational problem into an unconstrained variational problem as follows:
(7)uk,ωk,λ=α∑k∂tδt+jπte−jωkt22+st−∑kKukt22+λt,st−∑kKukt.The specific implementation process of VMD is referred to [20–22, 41, 42]. By applying the alternating direction method of multipliers to iteratively update uk, ωk, and λt, the narrowband components and their associated center frequency can be calculated.Although VMD has a strong decomposition ability in processing the noisy signal, it cannot effectively suppress the pattern aliasing phenomenon if the parameters of VMD, such as the number of modesK and the value of the penalty factor α, are not properly set. [43]. Thus, parameters K,α determine the quality of the decomposed modes, i.e., whether the same frequency band signal can be decomposed onto a single mode; and how much signal loss and noise is retained in the bandwidth of each mode after signal decomposition.Compared with traditional intelligent optimization algorithms, the GWO algorithm provides unparalleled advantages [44, 45]. Nonetheless, the GWO algorithm also has some disadvantages, such as poor stability and easy falling into local optimization [46]. Therefore, an improved Grey Wolf Optimization (IGWO) algorithm is applied to hunt the optimal VMD parameter K,α. Based on the previous improved GWO work [29], tent chaotic mapping is further introduced to make the initial population with uniformly distributed diversity and enhance the global convergence speed of the IGWO algorithm. The optimization process of IGWO is led by alpha, beta, and delta wolves. When the maximum iteration is reached, the optimal VMD parameter K,α would be found by IGWO. In order to ensure that the reconstructed MMG signal after IGWO-VMD decomposition retains as many muscle activity feature as possible and removes as much noise as possible, the energy entropy is used as the fitness function of IGWO in this paper. The energy entropy reflects the uncertainty and complexity of the signal and is used to represent the distribution of the signal energy. A small energy entropy indicates the increased significance of the corresponding signal in the total energy [24, 43]. It can be described as follows:
(8)Huk=−PklogPk,where Huk is the energy entropy of the component uk, Pk=Ek/E is the percentage of the energy of the kth mode in the total signal energy, Ek=∑i=1Nuk2t and E=∑i=1Kui2t.
(9)Fitness=minKHuk.When the optimal fitness is obtained, the corresponding parametersK,α are optimal, and further, the signal can be properly decomposed into a series of modes by VMD.


### 2.3. CEEMDAN Algorithm

CEEMDAN, as an improved algorithm of EMD and EEMD, can adaptively decompose a complex signal into a series of IMFs. Therefore, CEEMDAN is suitable for analyzing nonlinear, nonstationary, and non-Gaussian signals [47]. Since CEEMDAN overcomes the shortcomings of EMD and EEMD, its decomposition process can effectively overcome the mode aliasing problem, the reconstruction error is almost zero, and the computational cost is greatly reduced.Ifxt is the target signal, the procedure of CEEMDAN is summarized as follows:DefineEj· as the operator that produces the jth IMF obtained by EMD decomposition, and let the jth IMF obtained by CEEMDAN decomposition be IMFjt¯.Step 1.
A signalxt+ε0wit is obtained by adding positive and negative paired Gaussian white noise to xt, where wit is the noise sequence added in the ith experiment and ε0 is the noise amplitude. The signal is decomposed by EMD I times to obtain the first mode IMF1t¯.
(10)IMF1t¯=1I∑i=1IIMF1it.
Then, the first residue can be obtained in(11)r1t=xt−IMF1t¯.Step 2.
A new signalr1t+ε1E1wit is obtained by adding the above white noise and continuing to implement decomposition to obtain the second IMF2t¯ component, as in
(12)IMF2t¯=1I∑i=1nE1r1t+ε1E1wit.
Then, the second residue can be obtained in(13)r2t=r1t−IMF2t¯.Step 3.
Repeat the above steps until the extremum points of the margin do not exceed two; in the case that this condition is met, the decomposition is terminated. Assuming thatkIMF¯s are obtained, the final residue is rkt. The original signal xt is decomposed as
(14)xt=∑k=1KIMFkt¯+rkt.
After the above steps, the original signal is finally decomposed into a series ofIMF¯ components and a residual component. A good decomposition effect can be obtained by adjusting various parameters, such as the noise standard deviation (Nstd), number of realizations (NR), and maximum number of iterations (MaxIter). In this paper, the parameters of CEEMDAN are set to Nstd=0.2,NR=100,andMaxIter=1000.


## 2.1. RLS Algorithm

PLI is an electromagnetic noise from electronic equipment and transmission lines during biosignal measurements, which reduces the quality of the signal [36]. Therefore, removing PLI from the original signal is a problem that cannot be ignored. PLI is usually described as additive noise, having a sinusoidal signal at a fixed frequency (50 Hz or 60 Hz) with an unknown phase and amplitude. Many conventional filtering techniques are used to solve this problem [37]. However, these filters will fail in the case of PLI with frequency drift. Ahmed et al. [38] compared LMS, NLMS, and RLS algorithms and found that RLS has the best performance in removing PLI. The recursive least square (RLS) method has been widely used because of its advantages of the easy numerical solution and fast parameter convergence [39]. The RLS algorithm minimizes the weighted linear least square cost function by calculating the filter coefficients recursively. Therefore, in this paper, RLS is used to remove the PLI from the measured MMG signals. The specific description of the RLS algorithm is as follows [38, 40]:
(1)ek=dk−yk,(2)yk=ωkxk,(3)Kk=Pk−1∙xkλ+xHk∙Pk−1∙xk,(4)Pk=λ−1Pk−1−λ−1KkxHkPk−1,(5)ωk+1=ωk+Kk∙ek,where ek is the error signal, dk is desired signal, yk is the output signal from the adaptive filter, xk is the filter input vector, ωk is the filter coefficients vector, Kk denotes the gain vector, and Pk is a correlation matrix of the input signal. In this paper, P0=I/c, where I is the unit matrix and c is set to 0.01; λ is the forgetting factor which is set to 0.99. In addition, filter order M is set to 2.For PLI in the measured MMG signals, it is difficult for a hardware to filter it out completely, but the adaptive filtering algorithm RLS cannot only suppress it well, but PLI-free MMG signals can pass through the filter unchanged.


## 2.2. IGWO-VMD Algorithm

VMD is a new time-frequency analysis method, which iteratively searches for the optimal solution of the variational modes, continuously updates each mode function and the center frequency, and obtains an ensemble of BLIMFs, avoiding the endpoint effect and spurious component problems during the iterative process. The constrained variational optimization problem using the alternate direction method of multipliers can be described as follows:(6)minφk,ωk∑i=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kut=st,where K represents the number of modes, k=1,2,⋯,K; uk and ωk are the modes and their center frequencies, respectively, uk= u1,⋯,uK, ωk=ω1,⋯,ωK; st is the VMD processed signal; δt is the Dirac distribution function; ∗ and ∂t denotes the convolution and partial differential operators, respectively; j is the imaginary unit, and ∙2 represents L2 norm.To solve the above constrained variational model, a quadratic penalty factorα and Lagrangian multiplier λt are introduced to transform the constrained variational problem into an unconstrained variational problem as follows:
(7)uk,ωk,λ=α∑k∂tδt+jπte−jωkt22+st−∑kKukt22+λt,st−∑kKukt.The specific implementation process of VMD is referred to [20–22, 41, 42]. By applying the alternating direction method of multipliers to iteratively update uk, ωk, and λt, the narrowband components and their associated center frequency can be calculated.Although VMD has a strong decomposition ability in processing the noisy signal, it cannot effectively suppress the pattern aliasing phenomenon if the parameters of VMD, such as the number of modesK and the value of the penalty factor α, are not properly set. [43]. Thus, parameters K,α determine the quality of the decomposed modes, i.e., whether the same frequency band signal can be decomposed onto a single mode; and how much signal loss and noise is retained in the bandwidth of each mode after signal decomposition.Compared with traditional intelligent optimization algorithms, the GWO algorithm provides unparalleled advantages [44, 45]. Nonetheless, the GWO algorithm also has some disadvantages, such as poor stability and easy falling into local optimization [46]. Therefore, an improved Grey Wolf Optimization (IGWO) algorithm is applied to hunt the optimal VMD parameter K,α. Based on the previous improved GWO work [29], tent chaotic mapping is further introduced to make the initial population with uniformly distributed diversity and enhance the global convergence speed of the IGWO algorithm. The optimization process of IGWO is led by alpha, beta, and delta wolves. When the maximum iteration is reached, the optimal VMD parameter K,α would be found by IGWO. In order to ensure that the reconstructed MMG signal after IGWO-VMD decomposition retains as many muscle activity feature as possible and removes as much noise as possible, the energy entropy is used as the fitness function of IGWO in this paper. The energy entropy reflects the uncertainty and complexity of the signal and is used to represent the distribution of the signal energy. A small energy entropy indicates the increased significance of the corresponding signal in the total energy [24, 43]. It can be described as follows:
(8)Huk=−PklogPk,where Huk is the energy entropy of the component uk, Pk=Ek/E is the percentage of the energy of the kth mode in the total signal energy, Ek=∑i=1Nuk2t and E=∑i=1Kui2t.
(9)Fitness=minKHuk.When the optimal fitness is obtained, the corresponding parametersK,α are optimal, and further, the signal can be properly decomposed into a series of modes by VMD.


## 2.3. CEEMDAN Algorithm

CEEMDAN, as an improved algorithm of EMD and EEMD, can adaptively decompose a complex signal into a series of IMFs. Therefore, CEEMDAN is suitable for analyzing nonlinear, nonstationary, and non-Gaussian signals [47]. Since CEEMDAN overcomes the shortcomings of EMD and EEMD, its decomposition process can effectively overcome the mode aliasing problem, the reconstruction error is almost zero, and the computational cost is greatly reduced.Ifxt is the target signal, the procedure of CEEMDAN is summarized as follows:DefineEj· as the operator that produces the jth IMF obtained by EMD decomposition, and let the jth IMF obtained by CEEMDAN decomposition be IMFjt¯.Step 1.
A signalxt+ε0wit is obtained by adding positive and negative paired Gaussian white noise to xt, where wit is the noise sequence added in the ith experiment and ε0 is the noise amplitude. The signal is decomposed by EMD I times to obtain the first mode IMF1t¯.
(10)IMF1t¯=1I∑i=1IIMF1it.
Then, the first residue can be obtained in(11)r1t=xt−IMF1t¯.Step 2.
A new signalr1t+ε1E1wit is obtained by adding the above white noise and continuing to implement decomposition to obtain the second IMF2t¯ component, as in
(12)IMF2t¯=1I∑i=1nE1r1t+ε1E1wit.
Then, the second residue can be obtained in(13)r2t=r1t−IMF2t¯.Step 3.
Repeat the above steps until the extremum points of the margin do not exceed two; in the case that this condition is met, the decomposition is terminated. Assuming thatkIMF¯s are obtained, the final residue is rkt. The original signal xt is decomposed as
(14)xt=∑k=1KIMFkt¯+rkt.
After the above steps, the original signal is finally decomposed into a series ofIMF¯ components and a residual component. A good decomposition effect can be obtained by adjusting various parameters, such as the noise standard deviation (Nstd), number of realizations (NR), and maximum number of iterations (MaxIter). In this paper, the parameters of CEEMDAN are set to Nstd=0.2,NR=100,andMaxIter=1000.


## 3. The Proposed Method Based on RLS, IGWO-VMD, and CEEMDAN

Combining the above analysis and theoretical basis, a novel noise suppression and artifact removal method based on RLS, IGWO-VMD, and CEEMDAN is proposed. The specific procedures are summarized as follows:Step 1.
Collect MMG signals during isometric muscle contractions in the laboratory.Step 2.
Remove the 50 Hz PLI from the MMG signal using RLS.Step 3.
Using the IGWO-VMD method to decompose the PLI-free MMG signal, a series of BLIMF components are obtained; in this paper, the parameters of IGWO-VMD are configured as shown in Table1.Step 4.
To extract the useful and effective BLIMF components, the correlation coefficient (CC) between each BLIMF and the PLI-free signals is calculated. By setting the correlation coefficient threshold (Ct) [48], the corresponding BLIMFs are selected as the low-frequency part and the high-frequency part, respectively. If CC is greater than Ct, the BLIMF will be selected as the low-frequency part. Otherwise, the BLIMF will be selected as the high-frequency part.
(15)Ct=maxCC10maxCC−3,where Ct represents the threshold value and maxCC denotes the maximum value of the correlation coefficient.Step 5.
Identify the effective components in the low-frequency part. A center frequency threshold of 5 Hz is set to identify the low-frequency part. Specifically, the low-frequency component less than 5 Hz is identified as an artifact component and will be abandoned, while the low-frequency component greater than 5 Hz is identified as the effective component and will be retained.Step 6.
Decompose the high-frequency part of the signal using CEEMDAN, and extract the effective components of the high-frequency part. SE [49] is a new approach for the measure of time-series complexity, similar to approximate entropy, but with better accuracy, less dependence on data length, and better computational consistency. Generally, the larger the SE value, the higher the complexity and randomness of the signal, and correspondingly, the more complex the dynamic system and signal generation mechanism. Therefore, in this step, the sample entropy (SE) of each IMF obtained by CEEMDAN is calculated. The SE threshold (St) range is proposed to select the effective IMF components.
(16)St=∑i=1LSEi−SE¯2L,where L represents the number of IMFs; SEi denotes sample entropy value of the ith IMF component, and SE¯ denotes the average of SE.
If SE is within the range1/5St,St, the relevant IMF components are maintained as the effective IMF components. Otherwise, the relevant IMF components are removed as useless components.Step 7.
Reconstruct the effective BLIMF components of the low-frequency part and the effective IMF components of the high-frequency part to realize useful signal extraction.Table 1
Parameter configuration of IGWO-VMD.
ParametersParameter values or rangesPopulation number10Maximum iterations number10Range ofK(2,15)Range ofα(500,10000)The flow chart of the proposed method in this paper is shown in Figure1. The pseudocode of IGWO-VMD algorithm is shown in Algorithm 1. The experiments in this study are carried out on a workstation with a 2.30 GHz Intel(R) i7 Core (TM), 8 GB RAM.Figure 1
Flowchart of the proposed method for MMG signal denoising.Algorithm 1: The pseudocode of IGWO-VMD algorithm.
Input: MMG signal segmentOutput: The global optimum (VMD optimal parameters [K, α])StartInitializing VMD parameters (tau=0; DC=0; init=1; tol=1e-7)Initializing IGWO parameters (dim=2; Maxiter=10; N=10, N wolves was initialized using Tent chaotic mapping)Initializing 𝑎, 𝐴 and 𝐶Xδ = the best search agentXβ = the second best search agentXδ = the third best search agentWhile (t <maximum iteration number)for each search agentVMD decomposition, and using Eq. (9) to calculate the fitness of each search agentUpdate the position of the current search agentend forUpdate a, A, and C using Eq.(10),(11),(16) in literature 29UpdateXδ, Xβ and XδUpdated optimal solutionX using Eq.(17)–(20) in literature 29t = t + 1end whilereturn the global optimum XδK,αEnd


## 4. Simulation Analysis Experiments with Synthetic Signals

To verify the effectiveness and advantages of the proposed method in noise suppression and artifact removal of the measured MMG signal, this section will simulate the proposed method through two simulated synthetic signals and compare with classical methods, such as EMD, VMD, wavelet, CEEMDAN, and IGWO-VMD. Both signals have some characteristics of MMG signal, such as simulated synthesized signal 1 with MMG-like nonlinear and nonstationary characteristics and simulated synthesized signal 2 with MMG-like frequency bandwidth. In the simulation experiments, the simulated synthesized signals are firstly suppressed with RLS for PLI and then processed with IGWO-VMD and CEEMDAN.The mode selection in EMD and CEEMDAN is also done by the SE range (1/5St,St). The parameter K in VMD is selected according to the number of IMFs in EMD, the parameter α is usually selected as 2000, and the BLIMFs are selected according to the central frequency and SE less than St. In addition, wavelet basis function db6, 6-layer decomposition, and soft threshold functions are used in wavelet processing, where the threshold values are as follows:
(17)λj=σ2lnNlnj+1,where σ=medianWj,k/0.6745, j represents the decomposition scale, N is the signal length, σ is a rough estimate of noise level, and Wj,k is the kth detail coefficient at scale j.To evaluate the performance of these denoising methods, three evaluation indicators, namely, signal-to-noise ratio (SNR), mean square error (MSE), and correlation coefficient (CC), are adopted. The calculation formulas of the three evaluation indicators are as follows:(18)SNR=10log∑i=1Nx2i∑i=1nxi−x∧i2,(19)MSE=1N∑i=1Nxi−x∧i2,(20)CC=∑i=1Nxi−x¯x^i−x^¯∑i=1Nxi−x¯2∙∑i=1Nx∧i−x∧¯i2,where xi is the original target signal, x^i is the denoised signal, and N is the signal length. x¯ represents the average value of the original target signal; x^¯ represents the average value of the denoised signal.If the SNR is larger, the denoised signal is closer to the original target signal, and the noise suppression effect is better. If the MSE is smaller, the deviation between the denoised signal and the original target signal is smaller, and the smoothness of the signal is better. If the value of CC is closer to 1, the similarity between the denoised signal and the original target signal is higher.


### 4.1. Simulation Experiment 1

#### 4.1.1. Construction of Simulation Signal

Since MMG signals have nonlinear and nonstationary characteristics, a similar synthetic signal is designed to simulate the noisy signal measured by the acceleration sensors. The synthetic signalss1t is made up of four components: s1t, sjt, spt, and gst, as shown in equation (21). s1t is the target signal which is a nonlinear, nonstationary waveform based on the coupled Van-der Pol oscillators [50]. sjt is baseline drift which is modulated by a DC component and a low-frequency sine signal with a frequency of 1 Hz. In this paper, baseline drift is used to simulate artifacts. spt is PLI which is a sine signal with a frequency of 50 Hz. gst is Gaussian white noise. The sampling frequency is 360 Hz, and signal length N is 2000.
(21)sjt=sin2πt+5,spt=0.95sin100πt+π4,ss1t=s1t+sjt+spt+gst.Figures2(a)–2(b) shows the time domain and frequency spectrum of the target signal s1t and the synthetic signal ss1t with a SNRin of 5 dB.Figure 2
The signal and its corresponding frequency spectrum: (a) the target signal; (b) the synthetic signal ss1; (c) the PLI-free signal.
(a)(b)(c)


#### 4.1.2. Denoising of Simulation Signal

From Figure2(b), it can be seen that the target signal s1t is contaminated by noise, PLI, and baseline drift. To better suppress noise and remove baseline drift, the synthetic signal ss1t is first processed using RLS for PLI, the removal PLI result is shown in Figure 2(c). It can be seen that PLI is completely invisible in the frequency spectrum, and the target signal is not changed after filtering, which better preserves the original target signal.IGWO is used to search for the optimal value of VMD parametersK,α by minimizing the fitness function. After iterations, the optimization results are obtained as K=13andα=9800. Then, VMD is carried out for the PLI-free signal with these parameters, and the decomposition of VMD is shown in Figure 3. The correlation coefficient (CC) and the center frequency (Fc) of the above 13 BLIMFs (abbreviated as Bs) are calculated, and the results are shown in Table 2. It can be seen that the 13 Bs have different center frequencies, and the mode aliasing problem is not obvious. According to formula (15), the calculated Ct value is 0.171, and the PLI-free signal is divided into the high-frequency part (B4-B13) and the low-frequency part (B1-B3). Then, in the low-frequency part, according to the Fc, less than 5 Hz is judged as the baseline drift, and higher than 5 Hz is judged as the low-frequency effective components. Then, the B2 and B3 are reconstructed to get the denoised signal by IGWO-VMD.Figure 3
VMD decomposition results of the PLI-free signal.Table 2
Correlation coefficients and the center frequency of the Bs.
ModeB1B2B3B4B5B6B7Fc0.01034.55019.683136.242754.983169.183386.9504CC0.41180.72240.55970.15270.14410.14470.1393ModeB8B9B10B11B12B13Fc102.5289112.7258122.9406137.5556148.7139163.3797CC0.14210.12640.13010.12610.13310.1400Further, CEEMDAN is conducted for the high-frequency part (the noise-dominant Bs), and the decomposition results are shown in Figure4. According to the formula (16), the SE of each IMF component is calculated as shown in Table 3. The SE threshold range is obtained as 0.0761–0.3803. Thus, IMF5, IMF6, IMF7, and IMF8 are selected as the effective IMF components, and the rest of the IMFs are considered as high-frequency noise components and useless components and then discarded.Figure 4
CEEMDAN decomposition results of the high-frequency part.Table 3
The SE of CEEMDAN decomposition results.
ModeIMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF8IMF9IMF10IMF11SE1.26180.61950.86620.46490.26680.24560.19810.11810.05580.01680.0015Finally, the effective components of the low-frequency part and the high-frequency part are reconstructed to obtain the signal processed by the proposed method.


#### 4.1.3. Comparison with Other Methods

To verify the effectiveness of the proposed method, the proposed method, EMD, VMD, wavelet, CEEMDAN, and IGWO-VMD are employed to denoise the same noisy signalss1t with a SNRin of 5 dB for comparison, and the results are recorded in Table 4.Table 4
The denoising performance indicators of various methods with a SNRin of 5 dB.
IndicatorEMDVMDWaveletCEEMDANIGWO-VMDThe proposed methodMSE0.31371.00241.00240.41230.22240.1991SNR10.83965.79515.79499.653012.334412.8140CC0.95840.86390.89670.94650.97250.9736As can be seen from Table4, these methods have some denoising effect. Nevertheless, the IGWO-VMD method is far superior to VMD in terms of performance indicators, which further indicates that using IGWO to search the VMD parameters can effectively improve the decomposition efficiency of VMD. The denoising results of both IGWO-VMD and the proposed method are significantly better than the other methods, but the proposed method in this paper has lower MSE and higher SNR and CC.Further, the denoising results of the above-mentioned methods for the synthetic signalss1t with SNRin values varying from 1 dB to 11 dB are presented in Figure 5. It can be seen that the proposed method obtained the best results, as expected, with different SNRin values. In addition, the designed IGWO-VMD method obtains suboptimal results. Specifically, the comparison of performance indicators at different SNRin values indicates that the proposed method obtains the optimal results with the smallest MSE, the largest SNR, and the largest CC in reconstructing signal s1t. In particular, when the SNRin value is 9 dB, the MSE is reduced to 0.1498, the SNRout is improved to 13.6085 dB, and the CC is improved to 0.9806 for the proposed method. Through these experimental comparisons, it is shown that the performance of the proposed algorithm in this paper outperforms other algorithms under different decibel noises. The above analysis shows that the proposed method is most suitable for denoising and baseline drift removal of the nonlinear and nonstationary signal.Figure 5
Denoising performance with different SNRin on the synthetic signal ss1t.


### 4.2. Simulation Experiment 2

#### 4.2.1. Construction of Simulation Signal

When human muscles are active, the internal muscle vibrations would produce MMG signals in a certain frequency range. The main components of MMG signals are distributed between 10 and 40 Hz and are contaminated by noise, PLI, and artifacts. A synthetic signalss2t are used to simulate this situation. The synthetic signal ss2t is also made up of four components: s2t, sjt, spt, and gst. s2t is the target signal which is a coupling frequency component that is used to simulate MMG signal with the frequency of 5-100 Hz. sjt, spt, and gst refer to Simulation Experiment 1. The sampling frequency is 1000 Hz, and signal length N is 2000.
(22)s2t=21+0.6sin20πtcos60πt+1.5sin8πt,ss2t=st+sjt+spt+gst.The synthetic signalss2t is used as an example to further illustrate the effectiveness of the proposed method. Figures 6(a)–6(b) show the time domain and frequency spectrum of the target signal s2t and the synthetic signal ss2t with a SNRin of 5 dB.Figure 6
The signal and its corresponding frequency spectrum: (a) the target signal; (b) the synthetic signal ss2; (c) the PLI-free signal.
(a)(b)(c)


#### 4.2.2. Denoising of Simulation Signal

The process of denoising and baseline drift removal is similar to Simulation Experiment 1. The PLI removal result by adaptive RLS filter is shown in Figure6(c). The decomposition results of IGWO-VMD and CEEMDAN are shown in Figures 7 and 8, respectively. The CC and Fc of each BLIMF are shown in Table 5. The SE of each IMF is shown in Table 6. In Table 5, the calculated CC threshold value is 0.1606, and the low-frequency effective components (B2 and B3) are obtained. Correspondingly, the B2 and B3 are reconstructed to get the denoised signal by IGWO-VMD. In Table 6, the threshold range of the SE is calculated as 0.0819–0.4096, and the high-frequency effective components (IMF5-IMF8) are selected. Finally, the effective components of the low-frequency part and the high-frequency part are reconstructed to obtain the signal processed by the proposed method.Figure 7
VMD decomposition results of the PLI-free signal with parameters obtained by IGWO,K is 9, and α is 9665.Figure 8
CEEMDAN decomposition results of the high-frequency part.Table 5
The CC and Fc of the IGWO-VMD decomposition results.
ModeB1B2B3B4B5B6B7B8B9Fc0.012330.5882104.1649177.8528230.7533275.0756319.2273380.4707446.7987CC0.42570.79540.17350.15070.14800.15100.14320.14580.1370Table 6
The SE of the CEEMDAN decomposition results.
ModeIMF1IMF2IMF3IMF4IMF5IMF6SE1.32180.59891.03870.46500.29510.2415ModeIMF7IMF8IMF9IMF10IMF11IMF12SE0.21410.12630.04520.02360.00780.0003


#### 4.2.3. Comparison with Other Methods

To verify the effectiveness of the proposed method, the same comparison method as Simulation Experiment 1 is adopted, and the performance indicators of different methods are recorded in Table7. It can be seen that for the processing of the synthetic signal ss2t, the proposed method significantly outperforms the other methods in terms of MSE, SNR, and CC.Table 7
The denoising performance indicators of various methods with a SNRin of 5 dB.
IndicatorEMDVMDWaveletCEEMDANIGWO-VMDThe proposed methodMSE0.61740.58610.87860.35240.22940.1234SNR5.82336.04944.29118.259010.123812.4949CC0.89050.88520.83360.92370.95560.9736Similar to Experiment 1, the synthetic signalss2t with different decibel noises is denoised by the above methods, and the results are shown in Figure 9. As expected, the proposed method obtains optimal results. Moreover, it can be observed from Figure 9 that all performance indicators of the proposed method show an almost linear change with increasing SNRin values; i.e., MES gradually decreases with increasing SNRin, SNRout gradually increases, and CC also gradually improves. In particular, when the SNRin value is 11 dB, the MSE is reduced to 0.0491, the SNRout is improved to 16.6394 dB, and the CC is improved to 0.9896 using the proposed method. It also further shows that the proposed method is suitable for processing strong noisy background signals with a certain bandwidth. The superiority and reliability of the proposed method are further proved by these experiments.Figure 9
Denoising performance with different SNRin on the synthetic signal ss2t.


### 4.3. Experimental Results

From the denoising results of the above two synthesized signals, it can be seen that EMD, VMD, wavelet, and CEEMDAN are not ideal for suppressing 50 Hz PLI and noise, as well as for removing artifacts. With the appropriate VMD parameters selected by the IGWO designed in this paper, the signal can be effectively decomposed and the denoising effect can be improved. However, using IGWO-VMD to reconstruct the signal by directly removing the high-frequency part causes some information to be lost in the denoised signal. Therefore, in this paper, the high-frequency part is further decomposed by CEEMDAN to obtain the effective components, and then the effective components are reconstructed with the effective components of the low-frequency part obtained by IGWO-VMD to obtain a higher quality processed signal, which is almost consistent with the target signal. All of these further demonstrate that the proposed method has a better denoising effect on noisy signals with nonlinearity, nonstationarity, and certain bandwidth, which proves the effectiveness of the proposed method in this paper.


## 4.1. Simulation Experiment 1

### 4.1.1. Construction of Simulation Signal

Since MMG signals have nonlinear and nonstationary characteristics, a similar synthetic signal is designed to simulate the noisy signal measured by the acceleration sensors. The synthetic signalss1t is made up of four components: s1t, sjt, spt, and gst, as shown in equation (21). s1t is the target signal which is a nonlinear, nonstationary waveform based on the coupled Van-der Pol oscillators [50]. sjt is baseline drift which is modulated by a DC component and a low-frequency sine signal with a frequency of 1 Hz. In this paper, baseline drift is used to simulate artifacts. spt is PLI which is a sine signal with a frequency of 50 Hz. gst is Gaussian white noise. The sampling frequency is 360 Hz, and signal length N is 2000.
(21)sjt=sin2πt+5,spt=0.95sin100πt+π4,ss1t=s1t+sjt+spt+gst.Figures2(a)–2(b) shows the time domain and frequency spectrum of the target signal s1t and the synthetic signal ss1t with a SNRin of 5 dB.Figure 2
The signal and its corresponding frequency spectrum: (a) the target signal; (b) the synthetic signal ss1; (c) the PLI-free signal.
(a)(b)(c)


### 4.1.2. Denoising of Simulation Signal

From Figure2(b), it can be seen that the target signal s1t is contaminated by noise, PLI, and baseline drift. To better suppress noise and remove baseline drift, the synthetic signal ss1t is first processed using RLS for PLI, the removal PLI result is shown in Figure 2(c). It can be seen that PLI is completely invisible in the frequency spectrum, and the target signal is not changed after filtering, which better preserves the original target signal.IGWO is used to search for the optimal value of VMD parametersK,α by minimizing the fitness function. After iterations, the optimization results are obtained as K=13andα=9800. Then, VMD is carried out for the PLI-free signal with these parameters, and the decomposition of VMD is shown in Figure 3. The correlation coefficient (CC) and the center frequency (Fc) of the above 13 BLIMFs (abbreviated as Bs) are calculated, and the results are shown in Table 2. It can be seen that the 13 Bs have different center frequencies, and the mode aliasing problem is not obvious. According to formula (15), the calculated Ct value is 0.171, and the PLI-free signal is divided into the high-frequency part (B4-B13) and the low-frequency part (B1-B3). Then, in the low-frequency part, according to the Fc, less than 5 Hz is judged as the baseline drift, and higher than 5 Hz is judged as the low-frequency effective components. Then, the B2 and B3 are reconstructed to get the denoised signal by IGWO-VMD.Figure 3
VMD decomposition results of the PLI-free signal.Table 2
Correlation coefficients and the center frequency of the Bs.
ModeB1B2B3B4B5B6B7Fc0.01034.55019.683136.242754.983169.183386.9504CC0.41180.72240.55970.15270.14410.14470.1393ModeB8B9B10B11B12B13Fc102.5289112.7258122.9406137.5556148.7139163.3797CC0.14210.12640.13010.12610.13310.1400Further, CEEMDAN is conducted for the high-frequency part (the noise-dominant Bs), and the decomposition results are shown in Figure4. According to the formula (16), the SE of each IMF component is calculated as shown in Table 3. The SE threshold range is obtained as 0.0761–0.3803. Thus, IMF5, IMF6, IMF7, and IMF8 are selected as the effective IMF components, and the rest of the IMFs are considered as high-frequency noise components and useless components and then discarded.Figure 4
CEEMDAN decomposition results of the high-frequency part.Table 3
The SE of CEEMDAN decomposition results.
ModeIMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF8IMF9IMF10IMF11SE1.26180.61950.86620.46490.26680.24560.19810.11810.05580.01680.0015Finally, the effective components of the low-frequency part and the high-frequency part are reconstructed to obtain the signal processed by the proposed method.


### 4.1.3. Comparison with Other Methods

To verify the effectiveness of the proposed method, the proposed method, EMD, VMD, wavelet, CEEMDAN, and IGWO-VMD are employed to denoise the same noisy signalss1t with a SNRin of 5 dB for comparison, and the results are recorded in Table 4.Table 4
The denoising performance indicators of various methods with a SNRin of 5 dB.
IndicatorEMDVMDWaveletCEEMDANIGWO-VMDThe proposed methodMSE0.31371.00241.00240.41230.22240.1991SNR10.83965.79515.79499.653012.334412.8140CC0.95840.86390.89670.94650.97250.9736As can be seen from Table4, these methods have some denoising effect. Nevertheless, the IGWO-VMD method is far superior to VMD in terms of performance indicators, which further indicates that using IGWO to search the VMD parameters can effectively improve the decomposition efficiency of VMD. The denoising results of both IGWO-VMD and the proposed method are significantly better than the other methods, but the proposed method in this paper has lower MSE and higher SNR and CC.Further, the denoising results of the above-mentioned methods for the synthetic signalss1t with SNRin values varying from 1 dB to 11 dB are presented in Figure 5. It can be seen that the proposed method obtained the best results, as expected, with different SNRin values. In addition, the designed IGWO-VMD method obtains suboptimal results. Specifically, the comparison of performance indicators at different SNRin values indicates that the proposed method obtains the optimal results with the smallest MSE, the largest SNR, and the largest CC in reconstructing signal s1t. In particular, when the SNRin value is 9 dB, the MSE is reduced to 0.1498, the SNRout is improved to 13.6085 dB, and the CC is improved to 0.9806 for the proposed method. Through these experimental comparisons, it is shown that the performance of the proposed algorithm in this paper outperforms other algorithms under different decibel noises. The above analysis shows that the proposed method is most suitable for denoising and baseline drift removal of the nonlinear and nonstationary signal.Figure 5
Denoising performance with different SNRin on the synthetic signal ss1t.


## 4.1.1. Construction of Simulation Signal

Since MMG signals have nonlinear and nonstationary characteristics, a similar synthetic signal is designed to simulate the noisy signal measured by the acceleration sensors. The synthetic signalss1t is made up of four components: s1t, sjt, spt, and gst, as shown in equation (21). s1t is the target signal which is a nonlinear, nonstationary waveform based on the coupled Van-der Pol oscillators [50]. sjt is baseline drift which is modulated by a DC component and a low-frequency sine signal with a frequency of 1 Hz. In this paper, baseline drift is used to simulate artifacts. spt is PLI which is a sine signal with a frequency of 50 Hz. gst is Gaussian white noise. The sampling frequency is 360 Hz, and signal length N is 2000.
(21)sjt=sin2πt+5,spt=0.95sin100πt+π4,ss1t=s1t+sjt+spt+gst.Figures2(a)–2(b) shows the time domain and frequency spectrum of the target signal s1t and the synthetic signal ss1t with a SNRin of 5 dB.Figure 2
The signal and its corresponding frequency spectrum: (a) the target signal; (b) the synthetic signal ss1; (c) the PLI-free signal.
(a)(b)(c)


## 4.1.2. Denoising of Simulation Signal

From Figure2(b), it can be seen that the target signal s1t is contaminated by noise, PLI, and baseline drift. To better suppress noise and remove baseline drift, the synthetic signal ss1t is first processed using RLS for PLI, the removal PLI result is shown in Figure 2(c). It can be seen that PLI is completely invisible in the frequency spectrum, and the target signal is not changed after filtering, which better preserves the original target signal.IGWO is used to search for the optimal value of VMD parametersK,α by minimizing the fitness function. After iterations, the optimization results are obtained as K=13andα=9800. Then, VMD is carried out for the PLI-free signal with these parameters, and the decomposition of VMD is shown in Figure 3. The correlation coefficient (CC) and the center frequency (Fc) of the above 13 BLIMFs (abbreviated as Bs) are calculated, and the results are shown in Table 2. It can be seen that the 13 Bs have different center frequencies, and the mode aliasing problem is not obvious. According to formula (15), the calculated Ct value is 0.171, and the PLI-free signal is divided into the high-frequency part (B4-B13) and the low-frequency part (B1-B3). Then, in the low-frequency part, according to the Fc, less than 5 Hz is judged as the baseline drift, and higher than 5 Hz is judged as the low-frequency effective components. Then, the B2 and B3 are reconstructed to get the denoised signal by IGWO-VMD.Figure 3
VMD decomposition results of the PLI-free signal.Table 2
Correlation coefficients and the center frequency of the Bs.
ModeB1B2B3B4B5B6B7Fc0.01034.55019.683136.242754.983169.183386.9504CC0.41180.72240.55970.15270.14410.14470.1393ModeB8B9B10B11B12B13Fc102.5289112.7258122.9406137.5556148.7139163.3797CC0.14210.12640.13010.12610.13310.1400Further, CEEMDAN is conducted for the high-frequency part (the noise-dominant Bs), and the decomposition results are shown in Figure4. According to the formula (16), the SE of each IMF component is calculated as shown in Table 3. The SE threshold range is obtained as 0.0761–0.3803. Thus, IMF5, IMF6, IMF7, and IMF8 are selected as the effective IMF components, and the rest of the IMFs are considered as high-frequency noise components and useless components and then discarded.Figure 4
CEEMDAN decomposition results of the high-frequency part.Table 3
The SE of CEEMDAN decomposition results.
ModeIMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF8IMF9IMF10IMF11SE1.26180.61950.86620.46490.26680.24560.19810.11810.05580.01680.0015Finally, the effective components of the low-frequency part and the high-frequency part are reconstructed to obtain the signal processed by the proposed method.


## 4.1.3. Comparison with Other Methods

To verify the effectiveness of the proposed method, the proposed method, EMD, VMD, wavelet, CEEMDAN, and IGWO-VMD are employed to denoise the same noisy signalss1t with a SNRin of 5 dB for comparison, and the results are recorded in Table 4.Table 4
The denoising performance indicators of various methods with a SNRin of 5 dB.
IndicatorEMDVMDWaveletCEEMDANIGWO-VMDThe proposed methodMSE0.31371.00241.00240.41230.22240.1991SNR10.83965.79515.79499.653012.334412.8140CC0.95840.86390.89670.94650.97250.9736As can be seen from Table4, these methods have some denoising effect. Nevertheless, the IGWO-VMD method is far superior to VMD in terms of performance indicators, which further indicates that using IGWO to search the VMD parameters can effectively improve the decomposition efficiency of VMD. The denoising results of both IGWO-VMD and the proposed method are significantly better than the other methods, but the proposed method in this paper has lower MSE and higher SNR and CC.Further, the denoising results of the above-mentioned methods for the synthetic signalss1t with SNRin values varying from 1 dB to 11 dB are presented in Figure 5. It can be seen that the proposed method obtained the best results, as expected, with different SNRin values. In addition, the designed IGWO-VMD method obtains suboptimal results. Specifically, the comparison of performance indicators at different SNRin values indicates that the proposed method obtains the optimal results with the smallest MSE, the largest SNR, and the largest CC in reconstructing signal s1t. In particular, when the SNRin value is 9 dB, the MSE is reduced to 0.1498, the SNRout is improved to 13.6085 dB, and the CC is improved to 0.9806 for the proposed method. Through these experimental comparisons, it is shown that the performance of the proposed algorithm in this paper outperforms other algorithms under different decibel noises. The above analysis shows that the proposed method is most suitable for denoising and baseline drift removal of the nonlinear and nonstationary signal.Figure 5
Denoising performance with different SNRin on the synthetic signal ss1t.


## 4.2. Simulation Experiment 2

### 4.2.1. Construction of Simulation Signal

When human muscles are active, the internal muscle vibrations would produce MMG signals in a certain frequency range. The main components of MMG signals are distributed between 10 and 40 Hz and are contaminated by noise, PLI, and artifacts. A synthetic signalss2t are used to simulate this situation. The synthetic signal ss2t is also made up of four components: s2t, sjt, spt, and gst. s2t is the target signal which is a coupling frequency component that is used to simulate MMG signal with the frequency of 5-100 Hz. sjt, spt, and gst refer to Simulation Experiment 1. The sampling frequency is 1000 Hz, and signal length N is 2000.
(22)s2t=21+0.6sin20πtcos60πt+1.5sin8πt,ss2t=st+sjt+spt+gst.The synthetic signalss2t is used as an example to further illustrate the effectiveness of the proposed method. Figures 6(a)–6(b) show the time domain and frequency spectrum of the target signal s2t and the synthetic signal ss2t with a SNRin of 5 dB.Figure 6
The signal and its corresponding frequency spectrum: (a) the target signal; (b) the synthetic signal ss2; (c) the PLI-free signal.
(a)(b)(c)


### 4.2.2. Denoising of Simulation Signal

The process of denoising and baseline drift removal is similar to Simulation Experiment 1. The PLI removal result by adaptive RLS filter is shown in Figure6(c). The decomposition results of IGWO-VMD and CEEMDAN are shown in Figures 7 and 8, respectively. The CC and Fc of each BLIMF are shown in Table 5. The SE of each IMF is shown in Table 6. In Table 5, the calculated CC threshold value is 0.1606, and the low-frequency effective components (B2 and B3) are obtained. Correspondingly, the B2 and B3 are reconstructed to get the denoised signal by IGWO-VMD. In Table 6, the threshold range of the SE is calculated as 0.0819–0.4096, and the high-frequency effective components (IMF5-IMF8) are selected. Finally, the effective components of the low-frequency part and the high-frequency part are reconstructed to obtain the signal processed by the proposed method.Figure 7
VMD decomposition results of the PLI-free signal with parameters obtained by IGWO,K is 9, and α is 9665.Figure 8
CEEMDAN decomposition results of the high-frequency part.Table 5
The CC and Fc of the IGWO-VMD decomposition results.
ModeB1B2B3B4B5B6B7B8B9Fc0.012330.5882104.1649177.8528230.7533275.0756319.2273380.4707446.7987CC0.42570.79540.17350.15070.14800.15100.14320.14580.1370Table 6
The SE of the CEEMDAN decomposition results.
ModeIMF1IMF2IMF3IMF4IMF5IMF6SE1.32180.59891.03870.46500.29510.2415ModeIMF7IMF8IMF9IMF10IMF11IMF12SE0.21410.12630.04520.02360.00780.0003


### 4.2.3. Comparison with Other Methods

To verify the effectiveness of the proposed method, the same comparison method as Simulation Experiment 1 is adopted, and the performance indicators of different methods are recorded in Table7. It can be seen that for the processing of the synthetic signal ss2t, the proposed method significantly outperforms the other methods in terms of MSE, SNR, and CC.Table 7
The denoising performance indicators of various methods with a SNRin of 5 dB.
IndicatorEMDVMDWaveletCEEMDANIGWO-VMDThe proposed methodMSE0.61740.58610.87860.35240.22940.1234SNR5.82336.04944.29118.259010.123812.4949CC0.89050.88520.83360.92370.95560.9736Similar to Experiment 1, the synthetic signalss2t with different decibel noises is denoised by the above methods, and the results are shown in Figure 9. As expected, the proposed method obtains optimal results. Moreover, it can be observed from Figure 9 that all performance indicators of the proposed method show an almost linear change with increasing SNRin values; i.e., MES gradually decreases with increasing SNRin, SNRout gradually increases, and CC also gradually improves. In particular, when the SNRin value is 11 dB, the MSE is reduced to 0.0491, the SNRout is improved to 16.6394 dB, and the CC is improved to 0.9896 using the proposed method. It also further shows that the proposed method is suitable for processing strong noisy background signals with a certain bandwidth. The superiority and reliability of the proposed method are further proved by these experiments.Figure 9
Denoising performance with different SNRin on the synthetic signal ss2t.


## 4.2.1. Construction of Simulation Signal

When human muscles are active, the internal muscle vibrations would produce MMG signals in a certain frequency range. The main components of MMG signals are distributed between 10 and 40 Hz and are contaminated by noise, PLI, and artifacts. A synthetic signalss2t are used to simulate this situation. The synthetic signal ss2t is also made up of four components: s2t, sjt, spt, and gst. s2t is the target signal which is a coupling frequency component that is used to simulate MMG signal with the frequency of 5-100 Hz. sjt, spt, and gst refer to Simulation Experiment 1. The sampling frequency is 1000 Hz, and signal length N is 2000.
(22)s2t=21+0.6sin20πtcos60πt+1.5sin8πt,ss2t=st+sjt+spt+gst.The synthetic signalss2t is used as an example to further illustrate the effectiveness of the proposed method. Figures 6(a)–6(b) show the time domain and frequency spectrum of the target signal s2t and the synthetic signal ss2t with a SNRin of 5 dB.Figure 6
The signal and its corresponding frequency spectrum: (a) the target signal; (b) the synthetic signal ss2; (c) the PLI-free signal.
(a)(b)(c)


## 4.2.2. Denoising of Simulation Signal

The process of denoising and baseline drift removal is similar to Simulation Experiment 1. The PLI removal result by adaptive RLS filter is shown in Figure6(c). The decomposition results of IGWO-VMD and CEEMDAN are shown in Figures 7 and 8, respectively. The CC and Fc of each BLIMF are shown in Table 5. The SE of each IMF is shown in Table 6. In Table 5, the calculated CC threshold value is 0.1606, and the low-frequency effective components (B2 and B3) are obtained. Correspondingly, the B2 and B3 are reconstructed to get the denoised signal by IGWO-VMD. In Table 6, the threshold range of the SE is calculated as 0.0819–0.4096, and the high-frequency effective components (IMF5-IMF8) are selected. Finally, the effective components of the low-frequency part and the high-frequency part are reconstructed to obtain the signal processed by the proposed method.Figure 7
VMD decomposition results of the PLI-free signal with parameters obtained by IGWO,K is 9, and α is 9665.Figure 8
CEEMDAN decomposition results of the high-frequency part.Table 5
The CC and Fc of the IGWO-VMD decomposition results.
ModeB1B2B3B4B5B6B7B8B9Fc0.012330.5882104.1649177.8528230.7533275.0756319.2273380.4707446.7987CC0.42570.79540.17350.15070.14800.15100.14320.14580.1370Table 6
The SE of the CEEMDAN decomposition results.
ModeIMF1IMF2IMF3IMF4IMF5IMF6SE1.32180.59891.03870.46500.29510.2415ModeIMF7IMF8IMF9IMF10IMF11IMF12SE0.21410.12630.04520.02360.00780.0003


## 4.2.3. Comparison with Other Methods

To verify the effectiveness of the proposed method, the same comparison method as Simulation Experiment 1 is adopted, and the performance indicators of different methods are recorded in Table7. It can be seen that for the processing of the synthetic signal ss2t, the proposed method significantly outperforms the other methods in terms of MSE, SNR, and CC.Table 7
The denoising performance indicators of various methods with a SNRin of 5 dB.
IndicatorEMDVMDWaveletCEEMDANIGWO-VMDThe proposed methodMSE0.61740.58610.87860.35240.22940.1234SNR5.82336.04944.29118.259010.123812.4949CC0.89050.88520.83360.92370.95560.9736Similar to Experiment 1, the synthetic signalss2t with different decibel noises is denoised by the above methods, and the results are shown in Figure 9. As expected, the proposed method obtains optimal results. Moreover, it can be observed from Figure 9 that all performance indicators of the proposed method show an almost linear change with increasing SNRin values; i.e., MES gradually decreases with increasing SNRin, SNRout gradually increases, and CC also gradually improves. In particular, when the SNRin value is 11 dB, the MSE is reduced to 0.0491, the SNRout is improved to 16.6394 dB, and the CC is improved to 0.9896 using the proposed method. It also further shows that the proposed method is suitable for processing strong noisy background signals with a certain bandwidth. The superiority and reliability of the proposed method are further proved by these experiments.Figure 9
Denoising performance with different SNRin on the synthetic signal ss2t.


## 4.3. Experimental Results

From the denoising results of the above two synthesized signals, it can be seen that EMD, VMD, wavelet, and CEEMDAN are not ideal for suppressing 50 Hz PLI and noise, as well as for removing artifacts. With the appropriate VMD parameters selected by the IGWO designed in this paper, the signal can be effectively decomposed and the denoising effect can be improved. However, using IGWO-VMD to reconstruct the signal by directly removing the high-frequency part causes some information to be lost in the denoised signal. Therefore, in this paper, the high-frequency part is further decomposed by CEEMDAN to obtain the effective components, and then the effective components are reconstructed with the effective components of the low-frequency part obtained by IGWO-VMD to obtain a higher quality processed signal, which is almost consistent with the target signal. All of these further demonstrate that the proposed method has a better denoising effect on noisy signals with nonlinearity, nonstationarity, and certain bandwidth, which proves the effectiveness of the proposed method in this paper.


## 5. Application to MMG Signals

To verify the application effect of the proposed method in practical work, MMG signals measured by the acceleration sensor ADXL335 are analyzed in this section. The MMG signal acquisition method refers to the previous work [29]. MMG signal segments are randomly selected from two subjects with different force situations. During isometric muscle contraction, two types of MMG signals can be collected: one is static muscle force MMG signal (SMMG) which is measured at 60% of the maximal voluntary contraction (MVC); the other is dynamic muscle force MMG signal (DMMG) which is measured at 10-60% MVC. Two healthy male subjects aged 23 and 43 years are free of neuromuscular and musculoskeletal diseases.


### 5.1. Experiment 1: The Measured MMG Signal Segments Are Selected from Subject A

The measured MMG signal segments of SMMG and DMMG from subject A are shown in Figure10. It can be seen that both signals have PLI, noise, and artifacts, which completely obscure the effective MMG information and seriously affect the interpretation and application of the MMG signals.Figure 10
The measured MMG signals during isometric muscle contraction: (a) SMMG; (b) DMMG.
(a)(b)Figure11 shows the denoised signals and their frequency spectrum processed by the proposed method, IGWO-VMD, and the classical methods (EMD, VMD, wavelet, and CEEMDAN). The artifacts are all well corrected to zero level. However, the classical methods are not ideal for removing noise, and these methods do not completely remove the 50 Hz PLI. Compared with the classical methods, the proposed method and IGWO-VMD show outstanding results, i.e., 50 Hz PLI and artifacts are significantly removed, noise is obviously suppressed, and the effective components of the MMG signals are almost unchanged in the frequency spectrum.Figure 11
The processed MMG signals with different methods: (a) SMMG; (b) DMMG.
(a)(b)Since the real MMG is not directly available, the power spectral density (PSD) is used to further evaluate the denoising effects of the proposed method and IGWO-VMD. PSD is an excellent performance measurement tool that can be used to understand the performance of filtering techniques for noise reduction. In contrast to the frequency spectrum, PSD emphasizes the analysis of the average energy in the frequency range, which is essentially a representation of the energy distribution in the frequency domain [51]. In this paper, the Welch method is used to calculate the PSD of MMG signals, so as to obtain the energy distribution at different frequencies.In Figure12, the blue line is the PSD curve of the original MMG signal; the red line is the PSD curve of the MMG signal obtained by the proposed method; the green line is the PSD curve of the MMG signal obtained by IGWO-VMD. As can be seen from Figure 12, the original MMG signal has a large energy distribution in the low-frequency band, which is mainly caused by DC and artifacts. For MMG signals between 5 and 100 Hz, MMG energy is mainly distributed in the range of 10-40 Hz [52]. It can be seen that using the proposed method and IGWO-VMD, the PSD of MMG signals decreases significantly below 5 Hz and above 40 Hz, artifacts and noise are well suppressed, and the PSD values are close to those of the original MMG in the frequency band of 10-40 Hz. The analysis shows that the proposed method and IGWO-VMD can retain useful information in the original signal. Nevertheless, comparing the PSD curves of the proposed method and IGWO-VMD, it can be found that the proposed method is optimal.Figure 12
PSD of the original MMG signals and processed MMG signals from subject A: (a) SMMG; (b) DMMG.
(a)(b)


### 5.2. Experiment 2: The Measured MMG Signal Segments Are Selected from Subject B

To further prove the effectiveness of the proposed method in the measured signals, the measured MMG signals of subject B are processed, and the denoising results are shown in Figures13 and 14. There are some differences between Figures 12 and 14, mainly because different subjects, with different physical qualities, induce different recruitment and firing rates of MUs during muscle contraction. Figure 13 presents the denoising effect of three SMMG signal segments and three DMMG signal processing, respectively. It can be seen from Figure 13 that the proposed method is very successful in noise suppression and artifact removal. In addition, the proposed method has better denoising performance than IGWO-VMD. Specifically, the proposed method has a slightly larger signal amplitude than IGWO-VMD, mainly because the IGWO-VMD reconstructed signal ignores some useful signals in the high-frequency part. Furthermore, it can also be seen from Figure 14 that the proposed method preserves the effective components of the signal better than IGWO-VMD, i.e., the PSD values of the proposed method are closer to those of the original MMG in the frequency band of 10-40 Hz. This experiment confirms again that the proposed method is effective in MMG signal denoising, not only extracting the effective components of the actual measured MMG signals but also maintaining more signal details.Figure 13
The processed MMG signal segments of subject B: (a) SMMG; (b) DMMG.
(a)(b)Figure 14
PSD of the original MMG signals and processed MMG signals from subject B: (a) SMMG; (b) DMMG.
(a)(b)


### 5.3. Experimental Results

The above two experimental results show that the proposed method in this paper outperforms the classical methods and IGWO-VMD in terms of noise suppression and artifact removal. In addition, the proposed method better maintains the main energy components of the original signal in the range of 10-40 Hz, which further indicates that the proposed method is most appropriate and effective in extracting the effective components of the measured MMG signals. Therefore, the proposed method could lay a good foundation for further MMG signal identification and application.


## 5.1. Experiment 1: The Measured MMG Signal Segments Are Selected from Subject A

The measured MMG signal segments of SMMG and DMMG from subject A are shown in Figure10. It can be seen that both signals have PLI, noise, and artifacts, which completely obscure the effective MMG information and seriously affect the interpretation and application of the MMG signals.Figure 10
The measured MMG signals during isometric muscle contraction: (a) SMMG; (b) DMMG.
(a)(b)Figure11 shows the denoised signals and their frequency spectrum processed by the proposed method, IGWO-VMD, and the classical methods (EMD, VMD, wavelet, and CEEMDAN). The artifacts are all well corrected to zero level. However, the classical methods are not ideal for removing noise, and these methods do not completely remove the 50 Hz PLI. Compared with the classical methods, the proposed method and IGWO-VMD show outstanding results, i.e., 50 Hz PLI and artifacts are significantly removed, noise is obviously suppressed, and the effective components of the MMG signals are almost unchanged in the frequency spectrum.Figure 11
The processed MMG signals with different methods: (a) SMMG; (b) DMMG.
(a)(b)Since the real MMG is not directly available, the power spectral density (PSD) is used to further evaluate the denoising effects of the proposed method and IGWO-VMD. PSD is an excellent performance measurement tool that can be used to understand the performance of filtering techniques for noise reduction. In contrast to the frequency spectrum, PSD emphasizes the analysis of the average energy in the frequency range, which is essentially a representation of the energy distribution in the frequency domain [51]. In this paper, the Welch method is used to calculate the PSD of MMG signals, so as to obtain the energy distribution at different frequencies.In Figure12, the blue line is the PSD curve of the original MMG signal; the red line is the PSD curve of the MMG signal obtained by the proposed method; the green line is the PSD curve of the MMG signal obtained by IGWO-VMD. As can be seen from Figure 12, the original MMG signal has a large energy distribution in the low-frequency band, which is mainly caused by DC and artifacts. For MMG signals between 5 and 100 Hz, MMG energy is mainly distributed in the range of 10-40 Hz [52]. It can be seen that using the proposed method and IGWO-VMD, the PSD of MMG signals decreases significantly below 5 Hz and above 40 Hz, artifacts and noise are well suppressed, and the PSD values are close to those of the original MMG in the frequency band of 10-40 Hz. The analysis shows that the proposed method and IGWO-VMD can retain useful information in the original signal. Nevertheless, comparing the PSD curves of the proposed method and IGWO-VMD, it can be found that the proposed method is optimal.Figure 12
PSD of the original MMG signals and processed MMG signals from subject A: (a) SMMG; (b) DMMG.
(a)(b)


## 5.2. Experiment 2: The Measured MMG Signal Segments Are Selected from Subject B

To further prove the effectiveness of the proposed method in the measured signals, the measured MMG signals of subject B are processed, and the denoising results are shown in Figures13 and 14. There are some differences between Figures 12 and 14, mainly because different subjects, with different physical qualities, induce different recruitment and firing rates of MUs during muscle contraction. Figure 13 presents the denoising effect of three SMMG signal segments and three DMMG signal processing, respectively. It can be seen from Figure 13 that the proposed method is very successful in noise suppression and artifact removal. In addition, the proposed method has better denoising performance than IGWO-VMD. Specifically, the proposed method has a slightly larger signal amplitude than IGWO-VMD, mainly because the IGWO-VMD reconstructed signal ignores some useful signals in the high-frequency part. Furthermore, it can also be seen from Figure 14 that the proposed method preserves the effective components of the signal better than IGWO-VMD, i.e., the PSD values of the proposed method are closer to those of the original MMG in the frequency band of 10-40 Hz. This experiment confirms again that the proposed method is effective in MMG signal denoising, not only extracting the effective components of the actual measured MMG signals but also maintaining more signal details.Figure 13
The processed MMG signal segments of subject B: (a) SMMG; (b) DMMG.
(a)(b)Figure 14
PSD of the original MMG signals and processed MMG signals from subject B: (a) SMMG; (b) DMMG.
(a)(b)


## 5.3. Experimental Results

The above two experimental results show that the proposed method in this paper outperforms the classical methods and IGWO-VMD in terms of noise suppression and artifact removal. In addition, the proposed method better maintains the main energy components of the original signal in the range of 10-40 Hz, which further indicates that the proposed method is most appropriate and effective in extracting the effective components of the measured MMG signals. Therefore, the proposed method could lay a good foundation for further MMG signal identification and application.


## 6. Conclusions

To improve the denoising performance of the measured MMG signals, a novel noise suppression and artifact removal method based on RLS, IGWO-VMD, and CEEMDAN is proposed in this paper. The proposed method is easy to use and can effectively remove noise from the signal and correct artifacts. The proposed method is compared with the classical methods and the IGWO-VMD method by a large number of repeating simulation experiments. The results show that the proposed method is superior to the classical methods and the IGWO-VMD method in terms of quantitative denoising performance indexes. In the actual MMG signal processing experiments, The proposed method not only effectively eliminated the noise and PLI of the measured signal but also well corrected the artifacts to zero level. In addition, compared with other methods, the internal mechanical vibration components of muscles are effectively extracted using the proposed method, maintaining the main energy components of the original MMG signals with almost no energy loss. Therefore, it is concluded that the proposed method in this paper is effective and feasible.Although the effectiveness of the proposed method has been verified by the denoising results of the synthetic and measured signals, the results still have some limitations; for example, the synthetic signal types are not comprehensive enough, and the measured signals are limited to healthy subjects. Therefore, the universality of the proposed method needs to be further investigated. In addition, the feasibility of the proposed method in practical engineering applications such as muscle function assessment and human intention recognition needs to be further tested. A further work will be carried out to investigate the universality of the proposed method and apply the denoised MMG signals to practical engineering applications.



---

*Source: 4239211-2022-09-06.xml*