Variational Mode Decomposition for Raman Spectral Denoising

1. Introduction

Raman spectroscopy is a vibrational spectrum that excites molecular motion by light and used to chemical analysis of the sample by interprets the interaction, which is used for non-destructive analytical tool [1,2,3]. It is extensively used for qualitative and quantitative characterization in the field of food [4,5], diseases [6,7,8], biochemistry [9,10] and materials science [11,12,13]. Such as, Dai et al. [4] identified the age of pericarpium citri reticulatae (PCR) products by Raman spectroscopy. Zhang et al. [7] recorded the comprehensive molecular information about tumors in situ by Raman spectroscopy. However, due to the influence of thermal noise and dark current, noise inevitably exists in Raman spectrum. Moreover, Raman signal is inherently weak, making useful information susceptible to noise [14]. Therefore, it is important for Raman spectrum to remove noise.

Several methods have been used for spectral denoising, such as Savitzky-Golay (SG) smoothing [15,16], Whittaker smoother (WS) [17], deep learning (DL) [18,19], discrete wavelet transformation (DWT) [20,21], and empirical mode decomposition (EMD) [22,23]. SG smoothing, also known as convolution smoothing, is a signal smoothing method based on polynomial approximation of local least squares proposed by Savitzky and Golay [15]. It is a weighted average method, emphasizing the role of center point. The advantage of SG smoothing is that not has delay, and it can handle the data missed at a short period of time well. Luo et al. [3] applied SG smoothing to the Raman spectra denoising of cefuroxime axetil tablets and mixture of pyrene and cyclodextrin, as well as the Raman spectra denoising of cefuroxime axetil and cervical cancer cells. The intensities of the peaks in the Raman spectrum denoised by SG smoothing partially decreases with the retention noise. Meanwhile, due to the difference in the change rate of background and Raman spectrum, there is an irregular deviation between background and Raman spectrum. For SG smoothing, it is only performed once for denoise of signal. Moreover, it is important to choose the window length for the denoising performance. The different widths of the moving window could produce different denoising effects for SG smoothing. Large window width may result in the loss of valid information while small window width may remain noise in spectrum [24]. WS is a signal smoothing method based on least squares penalties proposed by Whittaker. The program of WS is very concise and it has the ability of adaptive to boundaries automatically. Compared with SG smoothing, WS has fast computation speed and easy to implement [25]. However, WS is not widely used because it requires the intervals between adjacent collection points of the data must be equal. Due to the highly flexible characteristic, DL has become a hot topic and also can be used for signal denoising [18,19,26]. Horgan et al. [18] proposed a comprehensive deep learning framework for high-throughput molecular imaging through Raman spectroscopy that supports deep learning. It can effectively reconstruct Raman characteristics from low signal-to-noise ratio (SNR) Raman spectra. In the denoising process of DL, a model is built between the noised signals and the pure signals of the same samples in training set. Then the new noised signal is taken as the input to the model and the denoised signal is obtained. Although excellent denoising ability was demonstrated, DL denoising requires a lot of training samples. Furthermore, sometimes it is difficult to obtain the pure signals.

DWT and EMD are both decomposition-based denoising methods, which decompose the original signal into a certain number of components with different frequencies. Compared with DL, the two methods do not need training samples and pure signals. For DWT denoising, the original signal was decomposed into a series of details and approximations. The high frequency details are deleted and the denoised signal is reconstructed by inverse wavelet transform [27]. Dou et al. [28] applied DWT to the research of intelligent transportation system data denoising and compression. However, it is difficult to select the optimal parameters for DWT because the large number of wavelet functions and decomposition scales. An adaptive method that does not require a defined basis as DWT, EMD was proposed by Huang et al. in 1998 [22,23]. EMD has been applied to Raman spectral denoising. Leon-Bejarano et al. [29] applied EMD as an adaptive parameterless signal processing method to Raman spectrum denoising of biological samples. However, EMD has mode mixing and endpoint effect in the decomposition process.

To overcome the drawbacks of EMD, variational mode decomposition (VMD) was proposed by Dragomiretskiy et al. in 2014 [30]. VMD can determine the relevant frequency band adaptively and estimate the correspond modes uk simultaneously in the decomposition process [31]. Compared with EMD, VMD has solid mathematical theory and can improve mode mixing and endpoint effect [32]. It has been extensively used for bearing fault diagnosis [33,34,35], pipeline leak detection [36,37], tidal analysis [38], and wind speed forecasting [39,40]. However, few research uses VMD for spectral denoising in analytical chemistry.

In this research, VMD is introduced to Raman spectral denoising. The original spectrum is decomposed into a series of uk by VMD. The noise modes are deleted and the remaining modes are reconstructed to obtain the denoised spectrum. The feasibility and effectiveness of the method are verified by two artificial noised signals and two actual Raman spectra. EMD, SG smoothing and DWT are used to compare with the proposed method, and the denoising effect is verified by visualization and SNR.

2. Methods

Denoising by variational mode decomposition (VMD)

VMD is a decomposition method for adaptive non-recursive signals and a process of solving variational problems. The advantage of this method is that it has solid mathematical theory foundation and strong noise suppression ability. The original signal is decomposed into a series of modes uk, each mode is tight around a central frequency and with limited bandwidth. The sum of all modes is equal to the original signal. The bandwidth problem of modes is transformed into a constrained optimization problem for solving each mode. The above problem is solved by introducing a quadratic penalty term $\alpha$ and Lagrangian multiplier $\lambda$. The saddle point of the augmented Lagrangian function is calculated by alternate direction method of multipliers (ADMM). Finally, the final mode $u_k$ is obtained by Fourier transform of the analysis signal. The detailed iteration process is shown in Figure 1. The detailed algorithm please refers to reference [30].

In this research, VMD is introduced for denoising of Raman spectra. The spectrum is decomposed into a series of uk by VMD, all the frequcencies of the modes are sorted from low to high. The high frequency modes represent noise and low frequency modes represent useful information. After the high frequency modes are removed, the low frequency modes are reconstructed to obtain the denoised spectrum.

2.2. Denoising by empirical mode decomposition (EMD)

As comparision, EMD is also used for the two artificial signals and two actual Raman spectra. Based on the inherent characteristics of the signal, the original signal is decomposed into a series of intrinsic mode functions (IMFs) and a residual (r) by EMD. The oscillation frequencies of the modes decrease as the order of IMFs increase. After the decomposition, the high frequency modes are removed as noise and the remaining modes are reconstructed to obtain the denoised signal.

2.3. Denoising by SG smoothing

SG smoothing is a widely used denoising method in spectra analysis. Thus, the method is also used to compare with VMD. SG smoothing is polynomials least square fitting of the data in the moving window, which is essentially a weighted average with more emphasis on the role of the central point. A symmetric window with length of i=2w+1 is used, where i is window size and w is the half window width. The output value of the filter is equal to the value of the polynomial at the center point. The next output value of the filter is calculated by shift the window sample and the step is repeated to obtain the denoised signal.

2.4. Denoising by DWT

Discrete wavelet transformation (DWT) implements the wavelet transformation by discrete translation factor and telescopic factor. The addition of these wavelets can be regarded as the decomposed signal. First, DWT is applied to the original signal, a series of approximations and details are obtained. Second, after determining the appropriate threshold, the noise in the signal is removed by processing the details with the hard or soft thresholding strategy. Finally, the processed wavelet coefficients are reconstructed to obtain the denoised signal.

2.5. Evaluation of denoising performance

In order to evaluate the denoising performance of the method, the figure of spectra before and after denoising and SNR are used. SNR is defined as

$$\mathrm{SNR}=10\lg \left\{{\displaystyle \sum }_s{f}^2\left(x\right)/{\displaystyle \sum }_s{\left[\widetilde{f(}x)-f\left(x\right)\right]}^2\right\}$$

(1)

where $f\left(x\right)$ is the original signal, $\widetilde{f(}x)$ is the denoised signal and s is the number of variable points [41]. A larger SNR means less noise mixed into the signal, and a smaller SNR means more noise mixed into the signal. Therefore, the SNR should be as high as possible.

3. Datasets

In order to verify the effectiveness of the proposed method, two artificial noised signals and two actual Raman spectra were employed.

3.1. The first artificial noised signal

The first artificial signal was constructed by three different frequency functions to verify the effect of VMD in decomposes the signal into different frequencies. It contains three cosine functions of different frequencies with a random noise, the three cosine functions are expressed by Equations (2)–(4). Figure 2a–e show the first artificial noised signal, functions of f1, f2, f3 and random noise, respectively.

$${\mathrm{f}}_1=\cos \left(2\mathsf{\pi }\times 2x-3\right)$$

(2)

$${\mathrm{f}}_2=\frac{1}{4}\cos \left(2\mathsf{\pi }\times 20x\right)$$

(3)

$${\mathrm{f}}_3=\frac{1}{6}\cos \left(2\mathsf{\pi }\times 50x\right)$$

(4)

3.2. The second artificial noised signal

Due to the characteristic peaks in Raman spectra are the Lorentz peaks, the second artificial signal is constructed by the Lorentz function. The pure signal of artificial signal is constructed by three Lorentz functions which have different peak widths and peak heights and a background function. The three Lorentz functions and the background function are expressed by Equations (5)–(8). The white Gaussian noise with 20dB is added to the pure signal to obtain the second artifical signal. Figure 3a,b show the curves of f4, f5, f6, f7 and the artificial noised signal, respectively. The artificial noised signal can be considered as artificial Raman noised spectrum, which is composed of 551 variables with a digitization interval of 0.1. As shown in Figure 3b, it is obvious that the artificial signal contains a lot of noise.

$${\mathrm{f}}_4=\frac{2}{\left({\left(5x-70\right)}^2+3\right)}$$

(5)

$${\mathrm{f}}_5=\frac{5}{\left({\left(x-3\right)}^2+3\right)}$$

(6)

$${\mathrm{f}}_6=\frac{4}{\left({\left(x+2\right)}^2+4\right)}$$

(7)

$${\mathrm{f}}_7=0.03\times x^{1/2}$$

(8)

3.3. The first actual Raman spectrum

The first actual Raman spectrum of MnCo single atom catalyst anchored on carbon and nitrogen materials (MnCo-ISAs/CN) was measured on a Raman spectrometer (XploRA PLUS, Japan). The spectrum was recorded with a digitization interval of 2.9 cm^-1 in the range of 750-2600 cm^-1, in total of 715 variables. Figure 4 shows the original Raman spectrum of MnCo-ISAs/CN, from which a large amount of noise can be seen. Hence, it is necessary to remove noise and retain valid information.

3.4. The second actual Raman spectrum

The second actual Raman spectrum of the iron-loaded carbon nanotubes (Fe-NCNT) synthesized by a convenient thermal treatment method was measured on a Raman spectrometer (XploRA PLUS, Japan). The spectrum was recorded with a digitization interval of 3.2 cm^-1 in the range of 100-3500 cm^-1, in total of 1351 variables. The Raman spectrum was shown in Figure 5.

4. Results and Discussion

4.1. Denoising of the first artificial signal

Figure 6a shows the decomposition result of VMD for the first artificial noised signal, which is decomposed into five modes u1-u5. It can be seen from Figure 6a that the frequency distribution of modes increases from low to high. The decomposition of VMD starts at the coarsest level. With the mode order increases, the oscillation frequencies of the modes also increase. The curve shapes of the three modes are almost as same as the three curves in Figure 2b, respectively. Obviously, the modes u1, u2 and u3 obtained by VMD correspond to the three functions f1, f2 and f3 that constitute the first artificial signal in Figure 2a. The modes u4 and u5 are random noise. Hence, the different frequency modes in the first artificial signal is separated completely by VMD.

As comparison, the first artificial noised signal is also denoised by EMD, SG smoothing and DWT. For the decomposition result of EMD, IMF1-IMF5 and a r are obtained, which are shown in Figure 6b. It is clear that the oscillation frequencies of the modes decrease as the order of IMFs increase. IMF2 and IMF3 are high frequency in the middle and low frequency at both ends, which are both high frequency mixing modes. The pure signal of first artificial signal is consisted by three different frequency curves. However, it is decomposed to four low frequency modes by EMD. Therefore, mode mixing and over-decomposition are produced for the first artificial signal by EMD.

Figure 7 shows the denoised signals by SG smoothing, DWT, EMD and VMD. For SG smoothing, with window size 16, most noise has been removed. However, some peaks of the SG smoothing signal are deviated or lower than those of the original signal. For the result of DWT, the signal still contains noise, which makes each peak not smooth enough. Visual inspection the decomposition result of EMD for the first artificial noised signal, IMF1 and IMF2 are the high frequency modes, which are removed as noise. IMF3-IMF5 and r are low frequency modes, which are reconstructed to obtain the denoised signal. Due to IMF2 is a mixing mode, the signal denoised by EMD has endpoint effect when IMF2 is removed, resulting in signal distortion at both ends. Furthermore, some information of the signal after denoising is lost due to mode mixing. For VMD, it is obvious that u1-u3 are low frequency modes with useful information, u4 and u5 are high frequency modes with little useful information. Hence, the artificial denoised signal is obtained by reconstruct u1-u3, while u4 and u5 are removed as noise. Obviously, the denoising result of VMD improves mode mixing and endpoint effects generated by EMD. The denoised signal by VMD almost overlaps with the original signal. Hence, the denoising effect of VMD is ideal. The SNR values of SG smoothing, DWT, EMD and VMD are 22.11, 23.11, 19.95 and 23.31, respectively. Therefore, for the first artificial signal, VMD has superiority in visual inspection and SNR compared with SG smoothing, DWT and EMD.

4.2. Denoising of the second artificial signal

Figure 8a shows the decomposition result of VMD for the second artificial noised signal, which is decomposed into five modes u1-u5. It can be seen from Figure 8a that u1 and u2 are low frequency modes, u3-u5 are the high frequency modes. The first artificial signal is effectively decomposed into different frequencies without mode mixing by VMD.

As comparison, the second artificial noised signal was also denoised by SG smoothing, DWT and EMD. Figure 8b shows the decomposition result of EMD, obtaining IMF1-IMF7 and a r. IMF1-IMF3 are the high frequency modes. IMF5-r are the low frequency modes. However, IMF4 is low frequency in the middle and high frequency at both ends. Thus, it indicates that IMF4 is a mixing mode. Hence, mode mixing is produced in the decomposition result of EMD.

Figure 9 shows the denoised signals by SG smoothing, DWT, EMD and VMD. For SG smoothing, the window size is selected as 33. Visual inspection the decomposition result of EMD for the second artificial noised signal, IMF3-r are used to reconstruct the denoised signal. For the result of VMD, u1 and u2 are reconstructed to obtain the artificial denoised signal and u3-u5 are removed as noise. It is obvious that most noise has been removed by SG smoothing, DWT, EMD and VMD. The signal denoised by VMD is smoother than that of SG smoothing and EMD. For EMD, some useful information of the signal after denoising is lost due to mode mixing, resulting in signal distortion in the range of 110th-190th points. As shown in Figure 9, the denoised signal has endpoint effect at left end. Furthermore, the peaks at 120th and 190th points of the signal denoised by SG smoothing, EMD and VMD are basically coincident. However, the signal after SG smoothing still contains noise at the two peaks mentioned above. For the second artificial signal denoised by DWT, the peaks located at 120th and 190th points are lower than those of the original signal. Moreover, the valley around 150th point after SG smoothing, DWT and EMD deviate significantly from the original signal position. For the small sharp peak in the range of 290-310 points, some information is lost after denoising by the four denoising methods. However, SG smoothing, DWT and EMD lost more information than that of VMD denoising. Thus, VMD shows better denoising performance than that of SG smoothing, DWT and EMD.

In order to evaluate the denoising performance of the four methods more intuitively, the SNR as a criterion is calculated. The SNR values of SG smoothing, DWT, EMD and VMD are 15.43, 15.36, 16.3 and 16.8, respectively. Thus, VMD has superiority in SNR compared with SG smoothing, DWT and EMD.

4.3. Denoising of the Raman spectrum of MnCo-ISAs/CN

To evaluate the denoising effect of the proposed method on the actual spectrum, the Raman spectrum of material MnCo-ISAs/CN was used for denoising. Figure 10a shows the decomposition result of VMD, which composes of six modes u1-u6. Obviously, u1 and u2 are the low frequency modes and u3-u6 are the high frequency modes by visual inspection. Furthermore, the decomposition result of Raman spectrum by VMD is easier to distinguish the boundary of noise modes and useful information modes.

For comparison, SG smoothing, DWT and EMD are also applied to the Raman spectrum denoising. The decomposition result of EMD is shown in Figure 10b, IMF1-IMF8 and a r are obtained. It is clear that IMF1-IMF3 are the high frequency modes and IMF5-r are the low frequency modes. However, IMF4 is low frequency in the range of 1400-1700 cm^-1 and high frequency in other ranges. Thus, IMF4 is a mixing mode. Hence, mode mixing is produced in the decomposition result of EMD.

Figure 11 shows the denoised signals by SG smoothing, DWT, EMD and VMD. For SG smoothing, the window size is selected as 33. The spectrum located at 1500cm^-1 still contains noise. For the denoising result of DWT, the denoised spectrum is smoother than that of SG smoothing. However, the peak located at 1350cm^-1 is not smooth enough. Moreover, the spectrum located at 1500cm^-1 still contains noise after DWT denoising. For the result of EMD, IMF5-IMF8 and r are reconstructed to obtain the Raman denoised spectrum and IMF1-IMF3 are removed as noise. However, IMF4 can be served as noise mode and also be served as useful information mode due to it is a mixing mode. If IMF4 is used to reconstruct the denoised spectrum as a useful information mode, the result is labeled as EMD1 and shown in Figure 11. The denoised spectrum still contains much noise and the endpoint effect has appeared at the right end. Moreover, for the result of EMD2, the peaks at 1350cm^-1 and 1600cm^-1 of the denoised spectrum are higher than that of other methods. If IMF4 is removed as noise mode, IMF5-r are used to reconstruct the denoised spectrum for EMD and labeled as EMD2 in Figure 11. It is obvious that the denoised spectrum is smoother than that of when IMF4 is reserved as useful information mode. However, the weak peak in the range of 1400-1700 cm^-1 is removed. As shown in Figure 11, the spectrum denoised by EMD has endpoint effect, resulting in spectrum distortion at right end. Furthermore, some useful information after denoising is lost due to mode mixing. For the result of VMD, u1 and u2 are reconstructed to obtain the Raman denoised spectrum, u3-u6 are removed as noise. VMD shows excellent denoising capability and almost all useful information is retained. Moreover, the spectrum denoised by VMD is smoother than that of EMD. Hence, it can be concluded that the denoising effect of VMD is better than that of SG smoothing, DWT, and EMD for the Raman spectrum.

4.4. Denoising of the Raman spectrum of Fe-NCNT

The Raman spectrum of Fe-NCNT was denoised by VMD. For the decomposition result of VMD, u1-u6 are obtained, which are shown in Figure 12a. It is clear that u1 and u2 are the low frequency modes and u3-u6 are the high frequency modes. The different frequency modes are obtained without mode mixing by VMD. Hence, it is easier to distinguish the boundary between noise modes and useful information modes.

In order to compare the denoising effect of the proposed method, SG smoothing, DWT and EMD are also used for the Raman spectrum denoising. For the decomposition result of EMD, IMF1-IMF9 and a r are obtained by EMD, which are shown in Figure 12b. IMF1-IMF3 are high frequency modes obviously and IMF5-r are low frequency modes. IMF4 is low frequency in the range of 1100-1700 cm^-1 and high frequency in other ranges. Thus IMF4 is a mixing mode. Therefor, the modes obtained by EMD have mode mixing obviously.

The denoised spectra by SG smoothing, DWT, EMD and VMD are shown in Figure 13. For SG smoothing, the window size is selected as 41. For EMD, since IMF4 is a mixing mode, it can be regarded as noise mode and also be regarded as useful information mode. If IMF4 is used to reconstruct the denoised spectrum, the result is labeled as EMD1 and shown in Figure 13. If IMF4 is removed as noise, the result is labeled as EMD2 in Figure 13. For the decomposition result of VMD, u1 and u2 are used to reconstruct the denoised spectrum and u3-u6 are removed as noise. It is obvious that most noise has been removed by SG smoothing, DWT, EMD1, EMD2 and VMD. The spectrum denoised by VMD is smoother than that of SG smoothing and EMD1. Morever, the denoising result of VMD improves mode mixing generated by EMD. For the result of EMD2, the denoised spectrum is distorted in the range of 1100-1750cm^-1 since the existence of mode mixing. For the denoised spectrum of DWT, the valley around 1500cm^-1 still contains noise. For the small peak located at 400 cm^-1, the spectrum retains the shape of the peak after VMD, EMD1 and EMD2. However, the small peak deviates from the original position after SG smoothing and DWT. Although the peak after VMD denoising is slightly wider than that of the original signal, the peak after DWT is wider than that of other denoising methods. The peaks at 1300cm^-1 and 1600cm^-1 of the spectrum by SG smoothing, DWT and VMD are basically coincident. However, for the denoising results of EMD1, the above two peaks are significantly higher than other methods. Therefore, the denoising effect of the proposed method is more ideal than that of SG smoothing, DWT and EMD.

5. Conclusions

VMD method was firstly introduced for denoising of Raman spectrum with noise. Two artificial noised signals and two actual Raman spectra were used to verify the denoising performance of the proposed method. The spectrum is decomposed by VMD into a series of modes, which are sorted from low frequency to high frequency. The high frequency modes are removed and the low frequency modes are reconstructed to obtained denoised spectrum. In order to verify the denoising effect of this method, SG smoothing, DWT and EMD are used as comparisons. Results show that VMD has highest SNR among the four denoising methods for the artificial denoised signal. Compared with the denoised result of EMD, VMD can decompose the signal into different modes with different frequencies, which improves the mode mixing and endpoint effects. Moreover, it is easier to distinguish the noise modes and useful information modes by VMD than that of EMD. The peaks of the spectrum denoised by VMD are smoother than that of the denoising spectra by SG smoothing and DWT. For the small sharp peaks, some information is lost after VMD denoising. However, SG smoothing, DWT and EMD lost more information than that of VMD. Hence, the proposed methods shows effectiveness both in SNR and denoising performance compared with SG smoothing, DWT and EMD. Therefore, VMD shows a high potential for signal denoising for Raman spectral analysis.