Denoising of Magnetocardiography Based on Improved Variational Mode Decomposition and Interval Thresholding Method

Abstract

Recently, magnetocardiography (MCG) has attracted increasing attention as a non-invasive and non-contact technique for detecting electrocardioelectric functions. However, the severe background noise makes it difficult to extract information. Variational Mode Decomposition (VMD), which is an entirely non-recursive model, is used to decompose the non-stationary signal into the intrinsic mode functions (IMFs). Traditional VMD algorithms cannot control the bandwidth of each IMF, whose quadratic penalty lacks adaptivity. As a result, baseline drift noise is still present or medical information is lost. In this paper, to overcome the unadaptable quadratic penalty problem, an improved VMD model via correlation coefficient and new update formulas are proposed to decompose MCG signals. To improve the denoising precision, this algorithm is combined with the interval threshold algorithm. First, the correlation coefficient is calculated, to determine quadratic penalty, in order to extract the first IMF made up of baseline drift. Then, the new update formulas derived from the variance that describes the noise level are used, to perform decomposition on the rest signal. Finally, the Interval thresholding algorithm is performed on each IMF. Theoretical analysis and experimental results show that this algorithm can effectively improve the output signal-to-noise ratio and has superior performance.

1. Introduction

In recent years, the research on signal processing, modeling, imaging theory, and methods related to bio-electromagnetism has become a hot topic. With the efforts of many experts and scholars, this field already has high-level research results. The magnetocardiography [1] signal plays an increasingly important role in heart disease diagnosis, which is detected with Superconducting Quantum Interference Devices (SQUID) [2] and has considerable advantages over electrocardiography (ECG) [3]. As the detecting instrument of the magnetocardiography signal, the SQUID operates from low to high temperature, and changes the number of channels from the original single channel into multiple channels. Magnetocardiography signals transmitted to the human chest surface are incredibly helpful toward both cardiac model reconstruction and clinical application [4,5]. The relationship between heart function and heart disease is studied by researching the characteristics of magnetic field strength changes at different locations; this type of study can be called interdisciplinary basic research. Generally, such measurements are conducted in order to detect small magnetic field signals in the presence of large background noise [6,7]. Removing background noise and recovering useful signals are chief objectives. For periodic signals, adequate suppression of uncorrelated noise may often be achieved by the signal averaging method [8]. If the signal and noise have separate bandwidths, one can use conventional frequency domain filtering [9] techniques. Adaptive filtering techniques [10] measure the noise level of the measurement signal channels by measuring the noise channels. However, these simple preprocessing methods have limited effectiveness. The wavelet transform method [11] for signal denoising is based on the use of a set of predefined basis functions, in order to decompose the measured signals and remove components corresponding to noise. The main disadvantage of this method is that the selection of wavelet basis seriously affects the denoising results. Empirical Mode Decomposition (EMD) [12,13] is one of the decomposition methods of signal denoising, and is widely used to decompose a signal into different modes recursively. This method is, however, prone to mode mixing, and limited by sensitivity to noise and sampling [14]. The mode mixing is significantly reduced by a modified noise-assisted data analysis method known as the Ensemble Empirical Mode Decomposition (EEMD) method [14,15]. The denoising principle of magnetocardiography (MCG) signals by EEMD based methods was reported in [16,17]. However, the decomposition results were unsatisfactory because of the low signal-to-noise ratio. In addition, the decomposition results of EMD and EEMD heavily depend on the extremum seeking algorithm and the ending criterion. A lack of mathematical approach and predefined filter boundaries reduce the accuracies of such detections [18]. Lately, based on the definition of intrinsic mode functions (IMF), a new adaptive decomposition method called Variational Mode Decomposition (VMD) [19] has been proposed. Supporting documents [20] proposed using the VMD method to denoise ECG signals. However, the research results showed that the decomposition results lack adaptability. The results of studies [21,22] showed that baseline drift noise was not filtered out by the VMD method. In practice, it is not always possible to have the first IMF to be a noise-only IMF.

In order to overcome the problems above, we propose an improved VMD method that determines the bandwidth of modes adaptively via the optimized quadratic penalty. The proposed correlation coefficient, between the IMF obtained by VMD [19,23] and the baseline drift model, is calculated repeatedly until the criterion is satisfied and the baseline drift noise is extracted. The new IMFs are then obtained by using new proposed update formulas that can be deduced by the relationship between the penalty factor and noise. The interval threshold method is used for the subsequent processing of each component, which removes noise components.

The rest of this paper is organized as follows: Section 2 introduces the data model required for the VMD algorithm. In Section 3, a new VMD scheme based on the correlation coefficient and new updated formulas is proposed. The application for denoising methods of MCG is shown in Section 4. Conclusions are given in Section 5.

2. Data Model

In the expression of the traditional EMD and EEMD methods, the IMF is defined as a function where the difference between the number of zeros and poles does not exceed one [24]. In recent studies, the definition of the modality is changed to amplitude-modulated-frequency-modulated (AM-FM) signal, defined as follows:

$$u_k\left(t\right)=A_k\left(t\right)\cos \left({\varphi }_k\left(t\right)\right)$$

(1)

In the above equation, the phase ${\varphi }_k\left(t\right)$ is a non-decreasing function, whose first derivative is $\varphi {’}_k\left(t\right)>0$, where the envelope $A_k\left(t\right)$ is non-negative; both the envelope $A_k\left(t\right)$ and the instantaneous frequency ${\omega }_k\left(t\right)$ vary much slower than the phase ${\varphi }_k\left(t\right)$ [25,26].

The Hilbert transform [27] is the convolution of a real function and the corresponding impulse response of $h\left(t\right)=1/\pi t$ in time domain. It is an all-pass filter, characterized by the transfer function $H\left(\omega \right)=-j\mathrm{sgn}\left(\omega \right)=-j\omega /\left|\omega \right|$ in frequency domain. The Hilbert transform of a purely real IMF $u_k\left(t\right)$ can be expressed as ${\tilde{u}}_k\left(t\right)$, and the complex-valued analytic signal is now defined as:

$$\begin{array}{cc}{u}_{k,A}\left(t\right)=u_k\left(t\right)+j{\tilde{u}}_k\left(t\right)& =A_k\left(t\right)\left[\cos \left(\varphi \left(t\right)\right)-j\sin \left(\varphi \left(t\right)\right)\right]\hfill \\ & =A_k\left(t\right)e^{-j\varphi \left(t\right)}\hfill \end{array}$$

(2)

where $\varphi \left(t\right)$ is the phase, while the amplitude is governed by the real envelope. The expression $A_k\left(t\right)$. $\omega \left(t\right)=d\varphi \left(t\right)/dt$ is the instantaneous frequency. The amplitude $A_k$ for kth IMF signal changes slowly enough.

Research has shown, on a sufficiently long interval, that the mode can be considered to be a purely harmonic signal. In other words, the newer definition of signal components is slightly more restrictive than the original one, and the VMD mode is the particular case of the EMD mode.

3. Proposed New VMD Scheme

VMD as a new decomposition method, is a process to solve variational problems based on classic Wiener filtering and Hilbert transformation. We can use the VMD method to decompose a multi-component signal into several band-limited modes non-recursively, which are redefined as IMFs. However, the VMD algorithm cannot extract baseline drift noise when decomposing MCG signals. As such, a new VMD method framework is proposed in this paper.

3.1. Eliminate Baseline Drift Noise Using Proposed Formulas

To overcome the unadaptable quadratic penalty problem, we propose an improved VMD method with correlation coefficient and new update formulas. First, we need to extract the expected baseline drift noise that will be included in the first mode. The steps are given as follows:

• Compute the associated analytic signal of each mode $u_k$ by means of the Hilbert transform, that is:

  $$u_{k,A}\left(t\right)=\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)$$

  (3)
• Mix each mode with an exponential adjustment to the respective estimated center frequency in order to shift the mode spectrum to “baseband”.

  $${\widehat{u}}_{k,A}\left(t\right)=\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}$$

  (4)

  where ${\omega }_k$ is the center frequency of the kth IMF $u_k\left(t\right)$.
• Estimate the bandwidth through the squared $L^2$-norm of the gradient. The expression of the constrained variational problem is as follows:

  $$\{\begin{array}{c}\underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{\min }\left\{{\displaystyle {\displaystyle \sum }_k}\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}{\Vert }_2^2\right\}\\ s.t. {\displaystyle {\displaystyle \sum }_{k=1}^K}{u}_k\left(t\right)=f\left(t\right)\end{array}$$

  (5)

  where $\left\{u_k\right\}=\left\{u_1,\dots ,u_K\right\}$ and $\left\{{\omega }_k\right\}=\left\{{\omega }_1,\dots ,{\omega }_K\right\}$ are shorthand notations for the set of all modes and their center frequencies. In order to render the problem unconstrained, a quadratic penalty term $\mathsf{\alpha }$ and Lagrangian multiplier $\mathsf{\lambda }$ are brought in. The quadratic penalty can encourage reconstruction fidelity, typically in the presence of additive Gaussian noise. The Lagrange equation can enforce constraints strictly. Therefore, we introduce the augmented Lagrange equation $\mathcal{L}$ as follow [28]:

  $$\mathcal{L}\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \left(t\right)\right)=\alpha {\displaystyle \sum }_{k=1}^K\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}{\Vert }_2^2+\Vert f\left(t\right)-{\displaystyle \sum }_{k=1}^K{u}_k\left(t\right){\Vert }_2^2+\lambda \left(t\right),f\left(t\right)-{\displaystyle \sum }_{k=1}^K{u}_k\left(t\right)$$

  (6)

Alternate direction method of multipliers (ADMM) is brought to solve the original minimization problem [29,30,31]. To update the mode $u_k$, we can get the equivalent minimization problem as the following:

$$u_k^{n+1}=\underset{u_k\in X}{\mathrm{argmin}}\left\{\alpha {\displaystyle \sum }_{k=1}^K\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}{\Vert }_2^2+\Vert f\left(t\right)-{\displaystyle \sum }_{i=1}^K{u}_i\left(t\right)-\frac{\lambda \left(t\right)}{2}{\Vert }_2^2\right\}$$

(7)

where $n$ is the number of iterations. Now, making use of the Parseval/Plancherel Fourier isometry under the $L^2$ norm, this problem can be solved in spectral domain. Then, performing a change of variables $\mathsf{\omega }\leftarrow \mathsf{\omega }-{\omega }_k$ in the first term, we can get the following expression:

$${\widehat{u}}_k^{n+1}=\underset{\widehat{u_k,}{u}_k\in X}{\mathrm{argmin}}\left\{\alpha \Vert j\left(\omega -{\omega }_k\right)\left[\left(1+\mathrm{sgn}\left(\omega \right)\right){\widehat{u}}_k\left(\omega \right)\right]{\Vert }_2^2+\Vert \widehat{f}\left(\omega \right)-{\displaystyle \sum }_i{\widehat{u}}_i\left(\omega \right)+\frac{\widehat{\lambda }\left(\omega \right)}{2}{\Vert }_2^2\right\}$$

(8)

Exploiting the Hermitian symmetry of the real signals, we can write both terms as half-space integrals, then making the negative frequencies of the first variation disappeared as follows:

$${\widehat{u}}_k^{n+1}\left(\omega \right)=\frac{\widehat{f}\left(\omega \right)-{{\displaystyle \sum }}_{i<k}{\widehat{u}}_i^{n+1}\left(\omega \right)-{{\displaystyle \sum }}_{i>k}{\widehat{u}}_i^n\left(\omega \right)+\frac{{\widehat{\lambda }}^n\left(\omega \right)}{2}}{1+2\alpha {\left(\omega -{\omega }_k^n\right)}^2}$$

(9)

This is clearly identified as Wiener filtering of the current residual, with signal prior $1/{\left(\omega -{\omega }_k^n\right)}^2$; the time domain mode is obtained as the real part of the inverse Fourier transform of this filtered analytic signal. In order to obtain each component, the center frequency, corresponding to each component, needs to be solved. The center frequency appears in the bandwidth prior, but not in the reconstruction fidelity term. The relevant problem thus reads:

$${\omega }_k^{n+1}=\underset{{\omega }_k}{\mathrm{argmin}}\{\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}{\Vert }_2^2\}$$

(10)

As before, the optimization can be taken place in Fourier domain, and we end up optimizing:

$${\omega }_k^{n+1}=\frac{{{\displaystyle \int }}_0^{\infty }\omega {\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega }{{{\displaystyle \int }}_0^{\infty }{\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega }$$

(11)

In general, the baseline drift frequency is lower than the low frequency component of the MCG signal. The center frequency of the first IMF is approximately zero, and we need to reduce the bandwidth of the first IMF until the signal and baseline drift are separated. We need to know that the larger the penalty factor, the narrower the mode bandwidth. After completing the above iterative process, we can get the final ${\widehat{u}}_k\left(\omega \right)$ and ${\omega }_k$. In order to extract low-frequency baseline drift noise, we need to follow the above process to obtain the first IMF:

$$u_1\left(t\right)=e^{j{\mathsf{\omega }}_{1}t}\left(\frac{1}{2\pi }{{\displaystyle \int }}_{-\infty }^{\infty }{\widehat{u}}_1\left(\omega \right)e^{j\mathsf{\omega }t}d\omega \right)$$

(12)

To understand the relationship between the first IMF and baseline drift noise, we propose the correlation coefficient to estimate the relationship. Assuming that the baseline drift noise model is $u_1{}^{’}\left(t\right)$ and the correlation coefficient ${\rho }^{\prime }$ are calculated by Equation (13):

$${\rho }^{\prime }=\frac{{{\displaystyle \sum }}_i(u_{1i}{}^{’}-\overline{u_1{}^{’}})\left(u_{1i}-\overline{u_1}\right)}{\sqrt{{{\displaystyle \sum }}_i{(u_{1i}{}^{’}-\overline{u_1{}^{’}})}^2{{\displaystyle \sum }}_i{\left(u_{1i}-\overline{u_1}\right)}^2}}$$

(13)

In the above formula, the variable is simplified for convenience. Where $i$ represents the length of the data, $\overline{u_1{}^{’}}$ and $\overline{u_1}$ represent the mean of the baseline drift noise model and the first IMF, respectively. We believe that signal and baseline drift noise can be separated when the correlation coefficient reaches a certain threshold $\rho$. We increase $\alpha$ according to the proposed formulas to repeat the above process until satisfying ${\rho }^{\prime }>\rho$. This satisfies the following formulas:

$$\{\begin{array}{c} im{f}_1=u_1\left(t\right) if{\rho }^{\prime }>\rho \hfill \\ {\alpha }^{m+1}={\alpha }^m+c_1,{\widehat{u}}_1^{n+1}\left(\omega \right)=\frac{\widehat{f}\left(\omega \right)-{{\displaystyle \sum }}_{i<1}{\widehat{u}}_i^{n+1}\left(\omega \right)-{{\displaystyle \sum }}_{i>1}{\widehat{u}}_i^n\left(\omega \right)+\frac{{\widehat{\lambda }}^n\left(\omega \right)}{2}}{1+2{\alpha }^{m+1}{\left(\omega -{\omega }_1^n\right)}^2} if {\rho }^{\prime }<\rho \end{array}$$

(14)

where $c_1$ is a constant and $m$ is the number of loop decomposition. After each updating of modes and center frequencies, the Lagrange multiplier $\widehat{\lambda }$ is also updated by Equation (15):

$${\widehat{\lambda }}^{n+1}\leftarrow {\widehat{\lambda }}^n+\tau \left(\widehat{f}-{\displaystyle \sum }_k{\widehat{u}}_k^{n+1}\right)$$

(15)

The updating stops until following equation is set up,

$${\displaystyle \sum }_k\Vert {\widehat{u}}_k^{n+1}-{\widehat{u}}_k^n{\Vert }_2^2/\Vert {\widehat{u}}_k^n{\Vert }_2^2<\epsilon$$

(16)

From the above, the first IMF $im{f}_1$ is the baseline drift noise.

3.2. Proposed Adaptive Decomposition

In order to get more reasonable decomposition results, we define $f_{new}=f-im{f}_1$. Document 19 proposed that the penalty factor introduced by the traditional VMD method is inversely proportional to the noise level in the signal. In a limited high frequency noise environment, we can assume that the weight of penalty is directly proportional to the power (which may be obtained by calculating variance) of each IMF. In order to further improve the adaptability of the penalty factor, we propose ${\alpha }_k=c_2\cdot D\left[u_k\left(t\right)\right]$. From the foregoing description, we can see that the low-frequency signal component has a large penalty factor, which can achieve low-frequency refinement and degrade the noise component in each signal mode. The original augmented Lagrange equation $\mathcal{L}$ becomes:

$$\begin{gathered} \mathcal{L}\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \left(t\right)\right)={\alpha }_k{\displaystyle \sum }_{k=1}^K\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]e^{-j{\omega }_{k}t}{\Vert }_2^2+ \\ \Vert f_{new}\left(t\right)-{\displaystyle \sum }_{k=1}^K{u}_k\left(t\right){\Vert }_2^2+\langle \lambda \left(t\right),f_{new}\left(t\right)-{\displaystyle \sum }_{k=1}^K{u}_k\left(t\right)\rangle \end{gathered}$$

(17)

Since the iterative solution process for each IMF component is performed in the frequency domain, the derivation of the penalty factor requires the use of a time-domain representation of each component. In order to obtain accurate decomposition results, we need to solve the penalty factors for each iteration, and then substitute the new penalty factors into the next iteration, which leads to the loop, which is very time-consuming. As such, we need to unify the derivation process into the frequency domain. Expanding the formula to solve the variance creates the following:

$$\begin{gathered} {\alpha }_k=\frac{c_2}{N}[{\left(u_{k1}-\overline{u_k}\right)}^2+{\left(u_{k2}-\overline{u_k}\right)}^2+\cdots +{\left(u_{kN}-\overline{u_k}\right)}^2] \\ =\frac{c_2}{N}\left(u_{k1}^2+u_{k2}^2+\cdots +u_{kN}^2\right)-\frac{1}{N^2}{\left(u_{k1}+u_{k2}+\cdots +u_{kN}\right)}^2 \end{gathered}$$

(18)

According to Parseval’s Theorem, we know that $\int {\left|u_k\left(t\right)\right|}^{2}dt=\int {\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega$ and the first item in the above formula can thus be converted into a frequency representation. By applying the Fourier transform of the second term to the above formula, we can get the frequency–domain description of penalty factor:

$${\alpha }_k=c_2\frac{N{{\displaystyle \int }}_0^{\infty }{\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega -2\pi {\widehat{u}}_k{}^2\left(0\right)}{2\pi N^2}$$

(19)

where $c_2$ is a constant and N is the length of the data. Combine the penalty factor into the previous update formulas, we can get the new update formulas for ${\widehat{u}}_k\left(\omega \right)$, ${\omega }_k$ and ${\alpha }_k$:

$$\{\begin{array}{c}{\widehat{u}}_k^{n+1}\left(\omega \right)=\frac{{\widehat{f}}_{new}\left(\omega \right)-{{\displaystyle \sum }}_{i<k}{\widehat{u}}_i^{n+1}\left(\omega \right)-{{\displaystyle \sum }}_{i>k}{\widehat{u}}_i^n\left(\omega \right)+\frac{{\widehat{\lambda }}^n\left(\omega \right)}{2}}{1+2{\alpha }_k^n{\left(\omega -{\omega }_k^n\right)}^2}\\ {\omega }_k^{n+1}=\frac{{{\displaystyle \int }}_0^{\infty }\omega {\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega }{{{\displaystyle \int }}_0^{\infty }{\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega }\\ {\alpha }_k^{n+1}=c_2\frac{N{{\displaystyle \int }}_0^{\infty }{\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega -2\pi {\widehat{u}}_k{}^2\left(0\right)}{2\pi N^2}\end{array}$$

(20)

The Lagrange multiplier update formula and iteration stop criterion remain unchanged.

3.3. Iterative Thresholding and Improved VMD Method

In this section, we obtain a new denoising process by combining the interval threshold and the improved variational mode decomposition. In this paper, the hard and soft thresholding methods, with multiples of Donoho-Johnstone threshold known as universal threshold parameters, are proposed, in order to cut off each IMF after performing the improved VMD. The threshold parameters are defined as:

$$T_k=C\frac{median\left(\left|u_k\left(t\right)-median\left(u_k\left(t\right)\right)\right|\right)}{0.6745}\sqrt{2\ln N}$$

(21)

where C is a constant and N is the length of the data. It is necessary to adopt a scale dependent threshold for each IMF instead of an identical universal threshold for all IMFs. The values of the scaling factor C are the range [0.6, 1.2]. We use the interval thresholding (IT) method to alleviate the catastrophic consequences caused by the direct application of thresholding. In this method, the first step is to find the zero points of each IMF. The second step is to compare the threshold and the extremum between two zero crossing intervals; if the extremum exceeds the threshold, it will allow all of the samples within the interval to be retained. The interval thresholding can be represented as:

$${\widehat{c}}_k\left(Z_i^k\right)=\{\begin{array}{c}{c}_k\left(Z_i^k\right), \left|c_k\left(r_i^k\right)\right|>T_k\\ 0 \left|c_k\left(r_i^k\right)\right|\le T_k\end{array}$$

(22)

where $i$ varies from 1 to $N_k$, $N_k$ indicates the number of zero crossings of the kth IMF, $c_k\left(r_i^k\right)$ indicates the sample at the time instance $r_i^k$ between the two successive zero points at $Z_i^k$ and $Z_{i+1}^k$, and $c_k\left(Z_i^k\right)$ refers to the samples from the instant $Z_i^k$ to $Z_{i+1}^k$ of the kth IMF.

Summarizing this denoising process as follows: First, the MCG signal with noise is decomposed into corresponding IMFs by improved VMD. The first order IMF, including the baseline drift noise, is eliminated. Second, the interval thresholding is performed on the rest IMFs. Clearly, the higher the order, the greater the frequency of the IMF. The detailed procedures are as follows:

• Acquire the MCG signal with noise and initialize the number of modes k, the default of the penalty factor $\alpha$ is 2000, the default of the bandwidth $\tau$ is 0.
• Determine the value of penalty factor by performing the improved VMD method based on the correlation coefficient and obtain the first IMF.
• Eliminate the first IMF that contains baseline drift noise. Decompose the rest signal with the final update formulas and obtain several IMFs.
• Perform the interval thresholding operation on the IMFs obtained from step 3.
• Add all the processed IMFs together and refactor the MCG signal.

In Figure 1 we show the flow chart of the improved VMD algorithm:

4. Results and Discussion

Three different types of noise have been added to the MCG signal, in order to investigate the effectiveness of denoising by an improved VMD and the interval thresholding method. The types of noise include a low frequency (0.3 Hz) sinusoidal signal for simulating the baseline drift, 50 Hz sinusoidal signal for simulating the interference at power line frequency, and high frequency random noise. In Figure 2, we compare the waveforms of an original simulation signal and a mixed signal. It is clearly seen that background noise affects signal analysis. It is very significant in the process of measuring MCG signals to detect heart disease.

The signal is denoised by different algorithms. In Figure 3a,b, we show the signal decomposition results obtained by the EEMD method. In Figure 3, according to the length of the data, the MCG signal is decomposed into ten IMFs, and the frequency of IMFs reduces as the order increases. In general, the low frequency smooth variation of the baseline is expected to be contained in the residue of the higher order IMFs. However, we cannot determine the low-frequency component that contains baseline drift noise. It is evident that we can extract signal characteristics from IMF2 to IMF6. Unfortunately, other IMFs may also contain useful signal components that are invisible to the human eye. Given this, we apply thresholds only to a few low order IMFs that contain contributions from the high frequency noise components, then exclude a few high order IMFs that contain the low frequency contents with a view, to ensure that the low frequency contents of the MCG signal are not affected or distorted by thresholding.

In Figure 4, we compare the original signal with the reconstructed signals obtained from EEMD based denoising methods using soft and hard thresholding. The results show that the hard threshold processing can reconstruct the QRS peak waves, but there are obvious errors in the reconstruction of other signal parts. The signal obtained by soft thresholding processing has errors compared with the original signal, especially the QRS peaks. It is seen from Figure 4 that EEMD algorithm is difficult to distinguish noise components from signal components. And the number of low frequency IMFs is too many to result in waveform distortion after performing thresholding operation. It may also be noted that, hard-IT (hard thresholding subsequently interval thresholding) is adopted for EEMD based denoising of the experimental data. At the same time, we can see that the signal waveform is not smooth and slightly distorted.

Based on the bandwidth of the measurement signal and multiple tests, the number of modes decomposed by VMD is assigned to 6. The original algorithm proposer studied some of the convergence characteristics of the VMD algorithm and its sensitivity to the initial conditions, then got relatively suitable initialization parameters. The initial value of quadratic penalty $\alpha$ is assigned to 2000, and the default of the bandwidth $\tau$ is 0. With this method, the MCG signal is divided into several frequency bands centering on respective center frequency, which benefits subsequent operation. As shown in Figure 5a, the MCG signal is divided into six IMFs. The baseline drift noise is found in IMF1. In Figure 5b, we compare the original signal with the reconstructed signals obtained from VMD methods with soft and hard thresholding. The results show that the baseline of the reconstructed signal is not uniform with the original signal. To make matters worse, there is serious distortion in the reconstructed signal from the soft thresholding processing. It is obvious that the algorithm cannot effectively remove baseline drift noise. The decomposition process of the algorithm lacks of adaptability. Improper decomposition results can easily lead to waveform distortion after threshold processing.

To get a better result, we use the proposed method to decompose the signal. According to many experiments, we choose the parameter value with better effect as the next simulation initial values. The threshold value $\rho$ of correlation coefficient is assigned as 0.95. The value of the scaling factor $c_1$ is the range 1500 and 2000. The value of $c_2$ can be set according to the detailed data features and is assigned as 3.5. As Figure 6 shown, the baseline noise is extracted. After removing the noise, the signal is divided into six modes adaptively. According to iterative formulas for multiple iterations, the penalty factors for six components are [1790, 987, 956, 218, 223, 208]. The decomposition results are shown in Figure 7a and the denoising results are shown in Figure 7b.

The range of the value of the penalty factors can support the adaptive decomposition result. From this result, it can be seen that the low-frequency IMF’s penalty factors are large, and the penalty factor increases approximately with the increasing of the center frequency, achieving the purpose of the meticulous decomposition of low-frequency signal component. In addition, we can obtain a set of different penalty factors by adjusting the size of $c_2$ to adjust the decomposition result. By comparing the reconstructed signal with the original signal, we can find that the fitting degree between the reconstructed signal and the original signal is good, and all three kinds of noise in the signal are effectively removed. It may also be noted that, signal to noise ratio improvement is much better of the hard interval thresholding (IT) method compared to the soft interval thresholding.

In order to better compare the performance of the algorithms, we use the root-mean-square error (the square root of the mean of the sum of squared residuals, RMSE) to characterize the fitting degree of the reconstructed signal and the original signal. In Figure 8, we compare the root-mean-square deviation (RMSE) of the three methods with the input signal-to-noise ratio (SNR). In the case of low input SNR, the RMSE of the improved VMD method is significantly less than the other two methods, and the method has better denoising performance even with the low input SNR.

In Figure 9, we compare the performance of the original VMD method, the improved VMD method, and the EEMD methods, using soft and hard thresholding in each case. For computing SNR, the logarithmic ratio of variance of a signal (from the beginning of P-wave to the end of T-wave for one cardiac cycle) to the variance of noise (from the end of T-wave to the beginning of P-wave, i.e., in the TP interval) has been taken. The reduction in the SNR for the soft thresholding method is due to the reduction of signal components of the lower frequency IMFs. Hereafter, hard-IT (hard thresholding is subsequently applied to the interval thresholds) is appropriate for both EEMD and VMD based denoising of the experimental data. It is seen from Figure 9 that the improved VMD method is capable of achieving better SNR when compared with EEMD and VMD methods.

From the above simulation results, we can see that the proposed algorithms have better denoising performance compared with EEMD and VMD methods. It should be noted, however, that the waveform is still distorted, even by the proposed method; this is a crucial detail in electrocardiographic-like signals. One encouraging factor, is that results of medical research have shown that the QRS spikes and S-T waves contain information on the main electrical function parameters of the heart. The MCG signal is filtered by the algorithm proposed in this paper, in order to obtain QRS spikes and S-T waves that approximate the original signal. Although there are slight disturbances in the QRS spikes and S-T waves obtained by filtering, it does not affect the calculation of heart related parameters! Measurements of magnetic field energy and current density remain accurate). In the case of the low input signal-to-noise ratio used in this paper, the SNR improvement of the proposed algorithm can be up to 20 dB. The algorithm filtering results can support feature extraction of MCG and detection of heart disease.

Despite these positive results, the MCG signals measured in a clinical non-magnetic shielding environment would still contain a large amount of background noise, which would still cause the waveform, denoised by the improved VMD method proposed in this paper, to produce severe distortion. In order to solve this problem, we need to study the spatial filtering technology in order to further denoise the MCG signal, in efforts to achieve the purpose of joint filtering in the time domain, frequency domain, and space domain.

5. Conclusions

The proposed method in this paper overcomes the unadaptable quadratic penalty problem of VMD, which improves the availability and precision of denoising of the MCG signal. This method adaptively adjusts the bandwidth of modes by repeatedly executing the VMD method with different quadratic penalties. The low-frequency noise is eliminated, according to the correlation coefficient for baseline drift noise, and the first mode. The new update formulas are used to decompose residual signals adaptively. Then, threshold processing is performed on each IMF to eliminate other noise. The simulation experiments show the superiority of the improved VMD in denoising performance of the MCG signal. The acceleration of the proposed method, and the suitable signal preprocessing methods, should be considered in future research and applications.