Nonstationary signal extraction based on BatOMP sparse decomposition technique

Sparse decomposition technique is a new method for nonstationary signal extraction in a noise background. To solve the problem of accuracy and efficiency exclusive in sparse decomposition, the bat algorithm combined with Orthogonal Matching Pursuits (BatOMP) was proposed to improve sparse decomposition, which can realize adaptive recognition and extraction of nonstationary signal containing random noise. Two general atoms were designed for typical signals, and dictionary training method based on correlation detection and Hilbert transform was developed. The sparse decomposition was turned into an optimizing problem by introducing bat algorithm with optimized fitness function. By contrast with several relevant methods, it was indicated that BatOMP can improve convergence speed and extraction accuracy efficiently as well as decrease the hardware requirement, which is cost effective and helps broadening the applications.


Background information.
Signal extraction and signal-noise separation are always two of the research focuses in modern signal processing 1 , which are commonly used in biomedical signal features extraction, vibration signal analysis, seismic signal detection 2 , sound signals recognition 3 . In practical applications, such signals are often submerged in a variety of vibration or electromagnetic noise, and the occurrence times of the target signals are random, which are typical nonstationary signals. Fourier transform is one of the most classic signal analysis and extraction method, but it cannot accurately describe nonstationary signals 4 . In recent years, new theories and technologies continue to appear in signal extraction. For example, wavelet decomposition (WD) 5,6 , empirical mode decomposition 7 , Hilbert-Huang Transform (HHT), variational mode decomposition (VMD) algorithm 8 . These methods need to meet certain conditions to work, for example the decomposition levels, modal number, and termination thresholds.
To achieve a more flexible, concise and adaptive signal decomposition, researchers proposed sparse decomposition. This method represents the signal with as few atoms as possible in a given redundant dictionary by matching pursuit (MP) algorithm 9 , which is a greedy algorithm for sparse decomposition. Various new evaluation criteria and basis pursuit, orthogonal matching pursuit algorithm (OMP) 10 , and time-frequency spectrum segmentation methods 11 were generated to select a set of optimal atoms from the constructed over-complete dictionary. In principle, if the dictionary redundancy is high enough and the iterations is large enough, the target signal can be perfectly extracted by OMP. On this basis, some general improved algorithms were proposed for example Regularized Orthogonal MP (ROMP) 12 and Compressive Sampling MP (CoSaMP) 13 . These methods require the signal Sparsity K for efficient execution, but K is generally unknown in practice. Sparsity Adaptive MP (SAMP) was proposed for signal reconstruction without prior information of the sparsity, but it is more complex than other greedy algorithms under large sparsity level 14 . And improper initial step size will lead to excessive decomposition for SAMP. The accuracy of signal sparse decomposition mainly depends on the redundancy and refinement accuracy of the redundant dictionary. Over or under estimation as well as long-time running will appear in these algorithms under the condition of large sparsity. Generally, the greater the redundancy and refinement, the greater the probability of accurate signal decomposition. However, for the greedy algorithm mentioned above, these are at the cost of algorithm efficiency. The accuracy and efficiency are exclusive.
Aiming at two main research hotspots including sparse decomposition algorithm and over-complete atom dictionary of signal sparse decomposition 15 , we designed two typical universal atoms, and proposed an adaptive feature-based atom construction method for the extraction of non-stationary signals with unknown sparsity. Redundant dictionary is obtained by extending the feature-based atoms, which can balance the completeness and redundancy. A signal matching tracking extraction algorithm was developed based on the bat algorithm and OMP, which successfully combined the accuracy and efficiency and could effectively realize nonstationary time domain signal extraction.
Classical signal sparse decomposition algorithms. Signal sparse decomposition represents a signal by specific combinations of some atoms in a dictionary. For a given dictionary, the optimal combination can be accurately determined when all possible combinations were calculated. However, exhausting all combinations in a dictionary is a non-deterministic polynomial problem that is almost impossible to achieve for large dictionary bases. So, the requirement was changed to finding a suboptimal combination from the dictionary with the lowest possible number of atoms and the smallest possible extraction error. This will reduce the computational complexity significantly, and the MP algorithm is one of the algorithms that can achieve this requirement.
Assume that the represented signal is x with length of N. Let ℝ denote the Hilbert space in which a dictionary matrix D composed of a set of vectors {g 1 , g 2 ,…,g n }. Each vector is an atom with the same length N and these vectors have been treated as normalized as g 2 = 1.
With ξ 1 = x, the MP algorithm selects one atom at a time from the dictionary matrix D that best matches x, satisfying (1), where i best is the index of the best matching atom in D. �·� is the inner product function.
The signal x is then decomposed into two parts, a sparse approximation ⌢ x and an approximation residual ξ 2 : Continues to select the atoms that best matches ξ 2 , iterating repeatedly and eventually the signal x can be approximated as a linear sum of these atoms: For MP algorithm, the non-orthogonality between the vertical projection of the signal (or residuals) on the selected atoms and the residuals will lead to suboptimal iterative results instead of the best optimal, and convergence requires many iterations. The OMP algorithm is the orthogonalization of all selected atoms at each step of the decomposition, which makes the convergence faster with the same accuracy requirement. The convergence process of MP and OMP are described by a dictionary D with length of three, as shown in Fig. 1. However, although the OMP algorithm reduced iterations to some extent, it had to calculate the current residual and the inner product of all atoms within the current dictionary during each iteration, resulting in unsatisfied effectiveness. Therefore, this paper introduced the bat algorithm (BA) to optimize the matching tracking algorithm.
Bat algorithm presentation. The basic flow of bat algorithm is as follows: (1) c i = x · g i best = max i∈(1,...n) x · g i , www.nature.com/scientificreports/ (2) The best position P n best is determined by the fitness function Fit n .
(3) Update the velocity and position of the individual bat: where r 1 was a random number, satisfying r 1 ∈[0,1]; f i was the search pulse frequency of the i-th bat; v i n denoted the velocity of the i-th bat in the n-th igeneration, P i n denoted the position of the i-th bat in the n-th igeneration; and P n best is the current global optimal solution. (4) Generate a random number r 2i ∈ [0,1] for each bat and update bat position according to (7).
where: η was a random number, satisfying η ∈ [− 1; 1] and Ā n was the mean fitness of the bat population. where, λ ∈ (0,1), γ > 0, when n → ∞ , A n i → 0 , r n i → r 0 . (8) Find the current matching atom based on the optimal solution. (9) The random perturbation of the current optimal solution in step 4 can effectively avoid the iterative result from falling into a local optimal solution, which helped to find the global optimal solution fast and accurate. (10) The Ackley function iss used to test the BA. The expression of the Ackley function is as follows: In this study, n = 2, c 1 = 20, e = 2.71289. the Ackley function was taken as the fitness function and the global minimum of this function was searched by the above methods. The particle swarm optimization (PSO) 16 , artificial fish school algorithm (AFSA) 17 and Cuckoo Search (CS) 18 are used for comparison. The population size and iteration numbers of these intelligent algorithms are the same to ensure rigorous comparison. The search paths and results are shown in Fig. 2.
The detailed values are shown in Table 1. The comparison of the tracking trajectory and the optimization results show that BA has advantage of high convergence speed and computational accuracy because the gradient of the optimization deviation is the largest and the optimization results are closest to the true value.

Methods
BatOMP sparse decomposition. General atomics designed for typical signals. For sinusoidal-like and one-sided decaying oscillatory signals, g-atoms are constructed: where c is the normalization factor to ensure that the original signal has the same energy as its sparse decomposition results; d: the attenuation factor; t: the sampling time; t 1 : start point of atomic appearance; t 2 : the ending point; f: signal frequency and ϕ: phase. The time domain waveforms of g-atoms with different parameters are shown in Fig. 3.  www.nature.com/scientificreports/ When the attenuation factor d = 0, the g-atom degenerates to standard sine wave; when d increases, the g-atom performs sinusoidal damped oscillation. Therefore, this atom has a strong match with sinusoidal signals, and single-sided oscillatory decay signals.
For kind of triangle waves, charge-discharge waves, and bilateral decay oscillation signals, t r atoms are constructed: where d 1 and d 2 are the bilateral damping factors; t 0 is the bilateral boundary of the atomic; [t 1 , t 2 ] is the atomic time range; and η is the bilateral scaling factor.
The time domain waveforms of t r -atoms with different parameters are shown in Fig. 4. When the bilateral scaling factor η = 0, the t r -atom degenerates to single-sided oscillating atom (reverse-order g-atom); when 0 < η < 1 and the atomic frequency is low enough, the t r -atom behaves as a charge-discharge triangle wave; when η = 1, the t r -atom with low-frequency behaves as a triangle-like wave, and behaves as bilateral oscillating decay signal with high-frequency.
The above analysis shows that the constructed g-atoms and t r -atoms are very flexible and could match almost typical testing signals by parameter adjustment.
, Figure 2. The optimal trajectory of different methods. (a,b) Show the 3D view and contour attempt of the optimal trajectory, respectively. Colors: black: the merit-seeking trajectories of PSO, blue: the merit-seeking trajectories of AFSA, green: the merit-seeking trajectories of CS, red: the merit-seeking trajectories of BA.  where: h(t) is the Hitch transform factor. Complex analytic signal as follows is constructed: where, A(t) is the amplitude function: and, ϕ(t) is the phase function: The instantaneous frequency of E i is given by (18): The base-atom is obtained with the time information gained by the correlation detection and localization algorithm and the time-frequency parameter information obtained by Hilbert transform as the reference. And the redundant dictionary of this feature atom is constructed by performing equal-step discrete expansion of the time-frequency parameters on both sides of the reference values. BatOMP improved sparse decomposition algorithm. The optimization-seeking process can be viewed as a global optimization problem. In order to solve the problems of large computation and low efficiency of existing matching tracking algorithms, the adaptive matching tracking algorithm, BatOMP, with fast convergence and accurate approximation is studied by combining BA into the OMP algorithm. For BatOMP, the bat individual positions P i represent the atoms column index in the redundant dictionary D, thus: g i = D(:,P i ). And for noise-containing signal extraction, the fitness function of the traditional sparse decomposition is improved to take the ratio of the ℓ-2 norm of the residual and the inner product as the fitness function. The target signal tends to be regular signals and most random noises obeys Gaussian distributions with zero mean error. So, the ℓ-2 norm of the former is greater than the latter. In consequence, the smaller the fitness indicates that the residual sequence www.nature.com/scientificreports/ contains smaller effective signal components and higher signal-to-noise separation. In addition, the larger the inner product, the better matchs between the atom and the residual. So, the optimal individual bat position P b is determined and saved according to (20). and where A is the matched dictionary, composed by the selected best matching atoms. The flow chart of BatOMP is as follows:

Experiments.
We constructed a nonstationary signal x to test the methods described above: where s is nonstationary target signal including pulse signal s 1 and partial discharge signal s 2 distributed in different regions, and ns is background noise subjecting to Gauss distribution. The sampling rate f s = 1500 Hz, sampling time T = 1 s, and SNR is 7.402 dB. Thus, the sequence length N is 1500.
The PC for the testing features a i7-8550U CPU Core(TM) @ 1.80 GHz with 16.0 GB RAM, 4 cores and 8 Logic processors, running the 64 bit operating system.

Time-frequency parameter calculation. First, time and frequency analysis was performed, and results
were shown in Fig. 6. From Fig. 6, after 0.4 s, s 1 decreases to zeros with the action of attenuation term.
Secondly, time information of the target signal was calculated by the correlation detection and localization algorithm described above, shown in Fig. 7. L 1 and L 2 are the calculated start and end indexes of the target segments, and L is the length of the segmentations. Accordingly, t 1 = L 1 /fs ≈ 0 s, t 2 = L 2 /fs ≈ 0.415 s for ŝ 1 , and t 1 = L 1 / fs ≈ 0.501 s, t 2 = L 2 /fs ≈ 0.671 s for ŝ 2 .
Different dictionaries construction and testing. After determining the key parameters of s, the redundant dictionary GT consisting of the two new atoms was created by dictionary training algorithm. And two (25)  Table 2.
Since the OMP method will iterate over the whole dictionary repeatedly, the extraction results are relatively accurate. In view of this, we take the results of OMP from statistical analyses to illustrate the performance of the different dictionaries, as shown in Fig. 8. The vertical axis represents the deviation between the extracted signal ŝ and actual signal s: We quantified the errors by RMS, the root mean square value of Amp err . The time-frequency parameters were not considered when building the DCT dictionary by MATLAB, so there was the maximum deviation in the correlative results. The atom expression of Gabor dictionary is as follows: We can see the lack of unilateral oscillation atoms by contrast with g-atom and t r -atom. Relatively, GT dictionary is completer and more accurate, so OMP based on GT dictionary gave minimum errors and the shortest optimizing time.
Algorithm performance testing. The algorithms involved in the article including MP, OMP, SAMP, CoSaMP, BatMP and the proposed BatOMP were carried out for performance comparison.
The extraction results and corresponding errors obtained by different methods were shown in Fig. 9. Figure 9a,c,e display the extracted signals ŝ by different methods based on DCT, Gabor and GT dictionary respectively. And Fig. 9b,d,f presente corresponding errors obtained by Eq. (26).
The efficiency analysis of different algorithms and dictionaries are shown in Fig. 10 and Table 3.
Because that the MP and OMP method traversed through the whole dictionary, the accuracies are relatively high and the latter is superior to the former.
The step length l of SAMP and sparsity q of CoSaMP were determined by expert experience. Because a certain amount atoms have to choose every time, there are clearly overextraction for these two methods. It's important to note that better parameters may be obtained by trial and error, but it is not suitable for real-time data processing.
For BatOMP, the best match atoms are determined by bat colony optimization. Every time before the searching, the bat individuals randomly scattered over the whole dictionary, and then gradually gather to the optimum solution through local optimization and global optimization. The optimal trajectories of ten bat individuals were randomly selected to show the convergence process, as shown in Fig. 11. The optimum solution is the index of the optimum matching atom in the redundant dictionary.
The difference between BatMP and BatOMP is similar to MP and OMP. BatOMP based on GT dictionary occupied the highest precision, probably because that MP and OMP took the inner product as the fitness function which leading to suboptimum for signal extraction. So, the results reflect the availability of the new fitness in some extent.
Moreover, the first four methods executed vast and complex matrix computations many times during the optimizing period, so they are time-consuming and require very high CPU occupancy rate compared with BatOMP.
(26) Amp err = s − s.  In other words, BatOMP can be widely used even on low lever machines. This is important for the occasions without algorithmic workstation and high-performance computer, i.e. field data processing or low cost testing.

Conclusion
For nonstationary signal extraction, the dictionary training algorithm based on feature parameters is firstly used to determine the key parameter range of feature atoms, which can effectively reduce the redundancy while ensuring the completeness of the redundant dictionary; the bat algorithm combined with OMP is proposed to transform the signal sparse decomposition problem into an optimization problem with ratio of the ℓ-2 norm of the residual and the inner product as the fitness, which can improve the efficiency of the sparse    www.nature.com/scientificreports/ decomposition algorithm. The experimental results showed that compared with other methods, the BatOMP algorithm is occupied with high efficiency, which can extract nonstationary signals form noise background without over constrained prior knowledges and avoid excessive decomposition. Testing results show that the proposed algorithm outperforms previous method in speeding up the convergence procedure and meanwhile ensuring high accuracy. Compared with the existing sparse decomposition algorithm, BatOMP requires much lower levels of hardware configuration. So, the new method will be helpful for to reducing data processing cost and enlarging the application fields.   Figure 11. A complete search process of the bat colony. Obviously, the bats gradually converge to the best solution from the original scattered position.