The proposed algorithm is implemented under the MATLAB 2018a environment and the computer numerical simulations are conducted using the Core i7-6700 3.41 GHz CPU with 8 GB RAM. In this paper, the electroencephalograms are acquired by a MUSE2 headband. The MUSE2 headband acquires the electroencephalograms at four different locations of the head. Two are at the forehead and two are behind the ear. In the following computer numerical simulations, the electroencephalograms are taken from the left forehead. The sampling rate of the electroencephalograms is 100 Hz and the acquired electroencephalograms are downloaded via the mobile application unit built in the MUSE2 headband. It is worth noting that the acquired electroencephalograms have significant changes in the amplitudes in the short periods of times at the locations where the spikes occur, and the locations of the spikes are random.
4.1. Computer Numerical Simulation Results
Figure 2a,b show a realization of a section of an acquired electroencephalogram in both the time domain and the frequency domain, respectively. That is,
Figure 2a shows a section of the entire electroencephalogram in the time domain.
Figure 2b shows the same section of the entire electroencephalogram in the frequency domain. It can be seen from
Figure 2a that there are some spikes corrupted in the acquired electroencephalogram. Also, it can be seen from
Figure 2b that the acquired electroencephalogram is a wide spectrum signal.
Figure 3 shows the locations of the spikes that need to be suppressed in the time domain.
Figure 3 also shows that there are 11 spikes in the acquired electroencephalogram and the acquired electroencephalogram has significant changes in the amplitudes in the short periods of times at the locations where the spikes occur.
Figure 4a shows the same section of the entire electroencephalogram after applying the ideal lowpass filtering in the time domain.
Figure 4b shows the same section of the entire electroencephalogram after applying the ideal lowpass filtering in the frequency domain. Here, since most of the information of the electroencephalogram is found in the 0–50 Hz frequency band,
is chosen. It can be seen from
Figure 4b that the noise outside the signal bands is removed after performing the ideal lowpass filtering. Although it can be seen from
Figure 4a that these 11 spikes are still there, the background noise is significantly suppressed.
Figure 5 shows the singular spectrum analysis components obtained by performing the singular spectrum analysis on the filtered electroencephalogram based on the discussion presented in
Section 2.1. It is worth noting that a too large value of
requires a heavy computational power while a too small value of
does not contain enough of the singular spectrum analysis components to be selected for performing the further processing. Therefore,
is chosen in the following computer numerical simulation results. This value is a good tradeoff between the above two factors. Also, as it is 0.5% of the length of the section of an acquired electroencephalogram which is less than half of the length of the section of an acquired electroencephalogram, it satisfies the property of the singular spectrum analysis.
Figure 6 shows the sums of the singular spectrum analysis components starting from the last singular spectrum analysis component. That is, the singular spectrum analysis component shown in
Figure 6a is the same as that shown in
Figure 5j. The singular spectrum analysis component shown in
Figure 6b is the sum of the singular spectrum analysis components shown in
Figure 5i,j. The singular spectrum analysis component shown in
Figure 6c is the sum of the singular spectrum analysis components shown in
Figure 5h–j. The rest of the singular spectrum analysis components shown in
Figure 6 are the sums of the singular spectrum analysis components shown in corresponding subfigures in
Figure 5. Finally, the singular spectrum analysis component shown in
Figure 6j is the sum of the singular spectrum analysis components in all the subfigures in
Figure 5. In fact,
Figure 6j is the filtered electroencephalogram.
Table 1 lists the values of the variances of the sums of the singular spectrum analysis components under the unit energy normalization. When comparing
Figure 6j with other figures in
Figure 6, it can be seen that there are spikes in the sums of the singular spectrum analysis components if the summations added from the last spectrum analysis component to the singular spectrum analysis components corresponding to the indices larger than or equal to four. Therefore, the threshold value defined on the variances of the sums of the singular spectrum analysis components under the unit energy normalization is set at
such that only the sum of the fifth singular spectrum analysis component to the last singular spectrum analysis component is used to generate the first scale of the electroencephalogram.
Figure 7a,b show the first scale of the electroencephalogram and the residue of the electroencephalogram, respectively. It can be seen from
Figure 7a that there is no spike in the first scale of the electroencephalogram, while all the spikes are found in the residue of the electroencephalogram.
Figure 8a,b show the low-rank component and the sparse component after performing the low-rank decomposition on the residue of the electroencephalogram. Here, since both the sparse component and the low-rank component have the same effects on the results, let
. Similarly, it can be seen from
Figure 8a that there is no spike in the low-rank component, while all the spikes are found in the sparse component. By adding the low-rank component to the first scale of the electroencephalogram, we obtain the second scale of the electroencephalogram which is shown in
Figure 9.
To further dig out the useful information of the sparse component, the singular spectrum analysis is applied to the sparse component to obtain a new set of the singular spectrum analysis components. The new set of the singular spectrum analysis components corresponding to the small eigenvalues are added to the current scale of the electroencephalogram to obtain the next scale of the electroencephalogram. Also, the low-rank decomposition is applied to the new residue of the electroencephalogram. The new low-rank components are also added to the current scale of the electroencephalogram to obtain the next scale of the electroencephalogram. Here, these scales of the electroencephalogram do not contain the spikes, but small amounts of useful information are lost. In order to recover the discarded useful information, the above procedures are repeated on the sparse component. From here, it is worth noting that the spikes are not gradually suppressed in each iteration. On the other hand, the useful information is added to the scales of the electroencephalogram.
Figure 10 and
Figure 11 show the 2
nth scales of the electroencephalogram in the time domain and the frequency domain, respectively. It can be seen from
Figure 10c that there is no spike after performing three iterations. On the other hand, it can be seen from
Figure 10d that there are spikes after performing four iterations. Compared with
Figure 11b,c, it can be seen that the electroencephalogram after performing two iterations lost some low-frequency information. Compared with
Figure 11c,d, the low-frequency information of the electroencephalogram after performing four iterations does not increase significantly.
Table 2 shows the variances of the 2
nth scales of the electroencephalogram under the unit energy normalization. Here, the threshold value defined on the variances of the 2
nth scales of the electroencephalogram under the unit energy normalization is set at 10 such that only three iterations are performed.
Figure 12 shows the despiked electroencephalograms based on the same set of the threshold values for the new realizations of an electroencephalogram. It can be seen from
Figure 12 that our proposed method still achieves the good despiked performances. This is because the same type of electroencephalograms is acquired by the same device. Therefore, the acquired electroencephalograms have consistent characteristics. After the threshold values are predefined, it is not required to change the predefined threshold values for the new realization of the electroencephalogram.
4.2. Comparisons to the Existing Methods
In order to demonstrate the effectiveness of our proposed method, our proposed method is compared to the following three states of the art methods. They are the empirical mode decomposition based denoising method [
13,
14], the singular spectrum analysis based denoising method [
15,
16,
17] and the low-rank decomposition based denoising method [
18,
19,
20].
Since the spikes contain the high-frequency contents of the electroencephalogram, the empirical mode decomposition based denoising method reconstructs the despiked electroencephalogram using the intrinsic mode functions corresponding to the low-frequency components [
13,
14]. Here, the intrinsic mode functions are obtained by applying the empirical mode decomposition on the original electroencephalogram. It is worth noting that if too many intrinsic mode functions are used for reconstructing the despiked electroencephalogram, then the despiked electroencephalogram will still contain the spikes. On the other hand, if very few intrinsic mode functions are used for reconstructing the despiked electroencephalogram, then the despiked electroencephalogram will be very smooth and some useful information will be lost. To determine which intrinsic mode functions are used for reconstructing the despiked electroencephalogram,
Figure 13 and
Figure 14 show the despiked electroencephalogram reconstructed using the last five intrinsic mode functions, the last seven intrinsic mode functions and the last nine intrinsic mode functions shown in the time domain and in the frequency domain, respectively. On the other hand,
Figure 15 shows the despiked electroencephalogram in the time domain obtained using our proposed approach. Our proposed approach can significantly suppress the spikes while most of the information is retained in the despiked electroencephalogram.
Figure 16a shows the despiked electroencephalogram in the frequency domain obtained using our proposed approach. The main lobe of the despiked electroencephalogram obtained using our proposed approach is mainly localized within 0–10 Hz while there is a large second lobe component localized between 10 Hz and 50 Hz. From here, we can reconstruct the despiked electroencephalogram by summing up the intrinsic mode functions from the one with the lowest frequency component to the one such that the cutoff frequency of the summed intrinsic mode functions is located at 10 Hz. From
Figure 14, we can see that the cutoff frequency of the despiked electroencephalogram reconstructed using the last five intrinsic mode functions is around 5 Hz, which is too low, resulting in a large amount of information being lost. On the other hand, the despiked electroencephalogram reconstructed using the last nine intrinsic mode functions contain too many high-frequency noises. Whereas, the despiked electroencephalogram reconstructed using the last seven intrinsic mode functions achieves the best tradeoff results. Therefore, the despiked electroencephalograms are reconstructed using the last seven intrinsic mode functions and it is shown in
Figure 14. However, it can be seen from
Figure 15 that the spikes are still found in the despiked electroencephalogram obtained based on the conventional empirical mode decomposition-based method. Compared to our obtained despiked electroencephalogram, it can be concluded that our proposed method outperforms the empirical mode decomposition-based method in terms of suppressing the spikes.
The reconstructed electroencephalograms obtained based on our proposed method and the conventional singular spectrum analysis-based method are shown in
Figure 17. Here, the conventional singular spectrum analysis-based method reconstructs the electroencephalogram by summing up the fifth singular spectrum analysis component to the last singular spectrum analysis component of the filtered electroencephalogram. On the other hand, as our proposed method also sums up the fifth singular spectrum analysis component to the last singular spectrum analysis component of the filtered electroencephalogram to generate the first scale of the electroencephalogram, the reconstructed electroencephalogram obtained based on our proposed method is similar to that based on the conventional singular spectrum analysis based method [
15,
16,
17]. However, it can be seen from
Figure 17 that the reconstructed electroencephalogram obtained based on the conventional singular spectrum analysis-based method is over-smoothen. This is because the information in the first four singular spectrum analysis components is lost. Compared with our reconstructed despiked electroencephalogram, it can be concluded that our proposed method outperforms the conventional singular spectrum analysis-based method in terms of retaining the characteristics of the electroencephalogram.
Finally, the low-rank decomposition-based method [
18,
19,
20,
21,
22] is also compared. Here, the low-rank component obtained by applying the low-rank decomposition to the filtered electroencephalogram is taken as the despiked electroencephalogram. This is because the spikes are found in the sparse component. The electroencephalograms obtained based on our proposed method and the conventional low-rank decomposition-based method are shown in
Figure 16. Similar to the conventional singular spectrum analysis approach [
15,
16,
17], it can be understood from
Figure 18 that the reconstructed electroencephalogram based on the low-rank decomposition-based method is over-smoothen. This is because the information in the sparse component is lost. Compared with our reconstructed despiked electroencephalogram, it can be concluded that our proposed method outperforms the low-rank decomposition-based method in terms of retaining the characteristics of the electroencephalogram.
It can be seen in
Figure 3 that there are 11 spikes in the electroencephalogram. Here, the objective is to suppress the magnitudes of these 11 spikes. However, it is worth noting that the waveforms used in the computer numerical simulations are the realizations of a practical electroencephalogram. There is no external noise added to the electroencephalogram. Since there is no clean electroencephalogram for performing the comparison, neither the noise level nor the signal to noise ratio cannot be measured. Therefore, it is difficult to have a comparison of quantitative performance. To address this issue,
Table 3 shows the differences in the magnitudes of the spikes at these 11 points before and after applying the above four despiked methods. It can be understood from
Table 3 that the empirical mode decomposition-based method does not effectively suppress the magnitudes of the spikes.
Figure 16 shows the reconstructed electroencephalograms in the frequency domain obtained based on the above four methods. It can be seen in
Figure 16 that our proposed method does not lose too much information in the low-frequency content of the electroencephalogram. On the other hand, the singular spectrum analysis-based method loses the information in the low-frequency content of the electroencephalogram, while the low-rank decomposition-based method loses the information in the high-frequency content of the electroencephalogram. Hence, it can be concluded that our proposed method outperforms the other three methods in terms of suppressing the spikes as well as retaining the information of the electroencephalogram.
To further demonstrate the outperformance of our proposed method, a synthetic signal is illustrated. The clean signal
is given by
Here,
,
,
and
. Also, the spike signal
composed of the convolution of two square waves with the amplitude and the position of the spike in the spike signal being random is corrupted to
. That is, the synthetic signal is
and it is shown in
Figure 19c. To perform the quantitative evaluation, the signal to noise ratio criterion defined as
is evaluated. By changing the range of the amplitude and the distribution of the position of the spike, the synthesized signals with different signal to noise ratios are obtained.
Table 4 shows the signal to noise ratios of the despiked signals obtained by the above four despiked methods. For our proposed method,
and
are chosen. Moreover, the threshold value for selecting the singular spectrum analysis components is set in such a way that the first four singular spectrum analysis components are chosen. Furthermore, the threshold value for terminating the algorithm is set in such a way that the eighth scale of the synthetic signal is the final despiked signal. For the singular spectrum analysis based denoising method, the sum of the first fourth singular spectrum analysis components is chosen as the despiked signal. For the empirical mode decomposition based denoising method, the sum of the last seven intrinsic mode functions is chosen as the despiked signal. For the low-rank decomposition based denoising method,
is chosen. All of these parameters are chosen as the same as that before. As shown in
Table 4, our proposed method outperforms the other three denoising methods in terms of suppressing the spikes. This is because the signal to noise ratio achieved by our proposed method is 20 dB larger than the other three denoising methods.