A Regularized Weighted Smoothed L0 Norm Minimization Method for Underdetermined Blind Source Separation

Compressed sensing (CS) theory has attracted widespread attention in recent years and has been widely used in signal and image processing, such as underdetermined blind source separation (UBSS), magnetic resonance imaging (MRI), etc. As the main link of CS, the goal of sparse signal reconstruction is how to recover accurately and effectively the original signal from an underdetermined linear system of equations (ULSE). For this problem, we propose a new algorithm called the weighted regularized smoothed L0-norm minimization algorithm (WReSL0). Under the framework of this algorithm, we have done three things: (1) proposed a new smoothed function called the compound inverse proportional function (CIPF); (2) proposed a new weighted function; and (3) a new regularization form is derived and constructed. In this algorithm, the weighted function and the new smoothed function are combined as the sparsity-promoting object, and a new regularization form is derived and constructed to enhance de-noising performance. Performance simulation experiments on both the real signal and real images show that the proposed WReSL0 algorithm outperforms other popular approaches, such as SL0, BPDN, NSL0, and Lp-RLSand achieves better performances when it is used for UBSS.


Introduction
The problem that UBSS [1,2] needs to address is how to separate multiple signals from a small number of sensors. The essence of this problem is to solve the optimal solution of the undetermined linear system of equations (ULSE). Fortunately, as a new undersampling technique, compressed sensing (CS) [3][4][5] is an effective way to solve ULSE, which makes it possible to apply CS to UBSS.
The model of CS is shown in Figure 1. According to this figure, it can be see that CS boils down to the form, where Φ = [φ φ φ 1 , φ φ φ 2 , ..., φ φ φ n ] ∈ R m×n is a sensing matrix with the condition of m n and φ φ φ i ∈ R m , i = 1, 2, ..., n, which can be further represented as Φ = ψ ψ ψϕ ϕ ϕ, while ψ ψ ψ is a random matrix, and ϕ ϕ ϕ is the sparse basis matrix. y ∈ R m is the vector of measurements. Moreover, b ∈ R m denotes the additive noise. To solve the ULSE in Equation (1), we try to recover the sparse signal x from the given {y, Φ} by CS. According to CS, this problem is transformed into solving the L 0 -norm minimization problem.
where denotes error. This rather wonderful attempt is actually supported by a brilliant theory [6]. Based on this theory, in the noiseless case, it is proven that the sparsest solution is indeed a real signal when x x x is sufficiently sparse and Φ satisfies the restricted isometry property (RIP) [7]: where K is the sparsity of signal x and δ K ∈ (0, 1) is a constant. In Equation (2), the L 0 -norm is nonsmooth, which leads an NP-hard problem. In practice, two alternative approaches are usually employed to solve the problem [8]: • Greedy search by using the known sparsity as a constraint; • The relaxation method for the P 0 .
For greedy search, the main methods are based on greedy matching pursuit (GMP) algorithms, such as orthogonal matching pursuit (OMP) [9,10], stage-wise orthogonal matching pursuit (StOMP) [11], regularized orthogonal matching pursuit (ROMP) [12], compressive sampling matching pursuit (CoSaMP) [13], generalized orthogonal matching pursuit (GOMP) [14,15], and subspace pursuit (SP) [16,17] algorithms. The objective function of these algorithms is given by: As shown in the above equation, the features of GMP algorithms can be concluded as: • Using sparsity as prior information; • Using the least squares error as the iterative criterion.
The advantage of GMP algorithms is that the computational complexity is low, but the reconstruction accuracy is not high in the noise case.
At present, the relaxation method for P 0 is widely used. The relaxation method is mainly divided into two categories: the constraint-type algorithm and the regularization method. The constraint-type algorithm can also be divided into L 1 -norm minimization methods and smoothed L 0 -norm minimization methods. The representative algorithm of the former is the BPalgorithm [18], and the latter is the smoothed L 0 -norm minimization (SL0) algorithm. For the SL0 algorithm, the objective function can be expressed as: where F σ (x) is a smoothed function, which approximates the L 0 -norm when σ → 0. Compared with L 1 or L p , a small σ is selected to make the function close to L 0 -norm [8]; therefore, F σ (x) are closer to the optimal solution. Based on the idea of approximation, Mohimani used a Gauss function to approximate the L 0 -norm [19], which is described as: According to the equation, we can know: when σ is a small enough positive value, the Gauss function is almost equal to the L 0 -norm. Furthermore, the Gauss function is differentiable and smoothed; hence, it can be optimized by optimization methods such as the gradient descent (GD) method. Zhao proposed another smoothed function: the hyperbolic tangent (tanh) [20], This smoothed function makes a closer approximation to the L 0 -norm than the Gauss function, as shown in [19], with the same σ; hence, it performs better in sparse signal recovery. Indeed, a large number of simulation experiments confirmed this view.
Another relaxation method is the regularization method. For CS, sparse signal recovery in the noise case is a very practical and unavoidable problem. Fortunately, the regularization method makes the solution of this problem possible [21,22]. The regularization method can be described as a "relaxation" approach that tries to solve the following unconstrained recovery problem: where λ > 0 is the parameter that balances the trade-off between the deviation term Φx − y 2 2 and the sparsity regularizer υ(x). The sparse prior information is enforced via the regularizer υ(x), and a proper υ(x) is crucial to the success of the sparse signal recovery task: it should favor sparse solutions and make sure the problem P υ can be solved efficiently in the meantime.
For regularization, various sparsity regularizers have been proposed as the relaxation of the L 0 -norm. The most popular algorithms are the convex L 1 -norm [22,23] and the nonconvex L p -norm to the p th power [24,25]. In the noiseless case, the L 1 -norm is equivalent to the L 0 -norm, and the L 1 -norm is the only norm with sparsity and convexity. Hence, it can be optimized by convex optimization methods. However, according to [8], in the noisy case, the L 1 -norm is not exactly equivalent to the L 0 -norm, so the effect of promoting sparsity is not obvious. Compared to the L 1 -norm, the nonconvex L p -norm to the p th power makes a closer approximation to the L 0 -norm; therefore, L p -norm minimization has a better sparse recovery performance [8].
In view of the above explanation, in this paper, a compound inverse proportional function (CIPF) function is proposed as a new smoothed function, and a new weighted function is proposed to promote sparsity. For the noise case, a new regularization form is derived and constructed to enhance de-noising performance. The experimental simulation verifies the superior performance of this algorithm in signal and image recovery, and it has achieved good results when applied to UBSS. This paper is organized as follows: Section 2 introduces the main work of this paper. The steps of the ReRSL0algorithm and the selection of related parameters are described in Section 3. Experimental results are presented in Section 4 to evaluate the performance of our approach. Section 5 verifies the effect of the proposed weighted regularized smoothed L 0 -norm minimization (WReSL0) algorithm in UBSS. Section 6 concludes this paper.

Main Work of This Paper
In this paper, based on the P F in Equation (9), we propose a new objective function, which is given by: According to this equation, We not only propose a smoothed function approximating the L 0 -norm, but also propose a weighted function to promote sparsity. This section focuses on the relevant contents of W = [w 1 , w 2 , ...w n ] T and H σ (x).

New Smoothed Function: CIPF
According to [26], some properties of the smoothed functions are summarized in the following: Property: Let f : R → [−∞, +∞] and, define f σ (r) ≈ f σ (r/σ) for any σ > 0. The function f has the property, if: (a) f is real analytic on (r 0 , ∞) for some r 0 ; It follows immediately from Property that { f σ (r)} converges to the L 0 -norm as σ → 0 + , i.e., Based on Property, this paper proposes a new smoothed function model called CIPF, which satisfies Property and better approximates the L 0 -norm. The smoothed function model is given as: In Equation (12), α denotes a regularization factor, which is a large constant. By experiments, the factor α is determined to be 10, which is a good result of the simulation. σ represents a smoothed factor, and when it is smaller, it will make the proposed model closer to the L 0 -norm. Obviously, where H σ (x) ≈ ||x|| 0 for small values of σ, and the approximation tends to equality when σ → 0. Figure 2 shows the effect of the CIPF model approximating the L0-norm. Obviously, the CIPF model makes a better approximation.
In conclusion, the merits of the CIPF model can be summarized as follows: • It closely approximates the L 0 -norm; • It is simpler in form than that in the Gauss and tanh function models.
These merits make it possible to reduce the computational complexity on the premise of ensuring the accuracy of sparse signal reconstruction, which is of practical significance for sparse signal reconstruction.

New Weighted Function
Candès et al. [27] proposed the weighted L 1 -norm minimization method, which employs the weighted norm to enhance the sparsity of the solution. They provided an analytical result of the improvement in the sparsity recovery by incorporating the weighted function with the objective function. Pant et al. [28] applied another weighted smoothed L 0 -norm minimization method, which uses a similar weighted function to promote sparsity. The weighted function can be summarized as follows: From the two weighted functions, we can find a phenomenon: a large signal entry x i is weighted with a small w i ; on the contrary, a small signal entry x i is weighted with a large value w i . By analysis, the large w i forces the solution x to concentrate on the indices where w i is small, and by construction, these correspond precisely to the indices where x is nonzero.
Combined with the above idea, we propose a new weighted function, which is given by: As for Candès et al., when the signal entry is zero or close to zero, the value of w i will be very large, which is not suitable for computation by a computer. Although Pant et al. noticed the problem and improved the weighted function to avoid it, the constant ζ depends on experience. Actually, the proposed weighted function can avoid the two problems. Moreover our weighted function can be satisfied with the phenomenon. When the small signal entry x i can be weighted with a large w i and a large signal entry x i can be weighted with a small w i , this can make the large signal entry and small signal entry closer. In this way, the direction of optimization can be kept as consistent as possible, and the optimization process tends to be more optimal. Therefore, the proposed weighted function can have a better effect.

WReSL0 Algorithm and Its Steps
Here, in order to analyze the problem more clearly, we rewrite Equation (10) as follows: σ . Therefore, we can obtain the gradient of CIPF, which is written as: According to Equation (15), as in [28], we can obtain: Solving the problem of ULSE is to solve the optimization problem in Equation (10). As for this problem, there are many methods, such as split Bregman methods [29][30][31], FISTA [32], alternating direction methods [33], gradient descent (GD) [34], etc. In order to reduce the computational complexity, this paper adopts the GD method to optimize the proposed objective function.
Given σ, a small target value σ min , and a sufficiently large initial value σ max , after referring to the annealing mechanism in simulated annealing [35], this paper proposes a monotonically-decreasing sequence {σ t |t = 2, 3, ..., T}, which is generated as: where γ = , θ is a constant that is larger than one, and T is the maximum number of iterations. Using such a monotonically-decreasing sequence can avoid the case of too small of a σ leading to the local optimum.
Similar to SL0, WReSL0 also consists of two nested iterations: the external loop, which begins with a sufficiently large value of σ, i.e, σ max , responsible for the gradually decreasing strategy in Equation (17), and the internal loop, which for each value of σ, finds the maximizer of According to the GD algorithm, the internal loop consists of the gradient descent step, which is given by: where d d d = g g g and µ denotes a step size factor. This part is similar to SL0, followed by solving the problem: arg min where x * denotes the optimal solution. By regularization, this form can be converted to another form as follows, arg min where λ is the regularization parameter, which is adapted to balance the fit of the solution to the data y and the approximation of the solution to the maximizer of H σ (x). Weighted least squares (WLS) can be used to solve this problem, and the solution is: By calculation, Equation (21) is equivalent to: where I n and I m are both identity matrices of size n × n and m × m, respectively. Therefore, we can obtain: According to the above analysis and derivation, we can get: The initial value of the internal loop is the maximizer of H σ (x) obtained for σ max . To increase the speed, the internal loop is repeated a fixed and small number of times (L). In other words, we do not wait for the GD method to converge in the internal loop.
According to the explanation above, we can conclude the steps of the proposed WReSL0 algorithm, which are given in Table 1. As for σ, it can be shown that function H σ (x) remains convex in the region where the largest magnitude of the component of x is less than σ. As the algorithm starts at the original value y, the above choice of σ 1 ensures that the optimization starts in a convex region. This greatly facilitates the convergence of the WReSL0 algorithm. Table 1. Weighted regularized smoothed L 0 -norm minimization (WReSL0) algorithm using the GD method.
, and T is the maximum number of iterations.

Selection of Parameters
The selection of parameters µ and σ will affect the performance of the WReSL0 algorithm; thus, this paper discusses the selection of these two above parameters in this section.

Selection of Parameter µ
According to the algorithm, each iteration consists of a descent step followed by a projection step. If for some values of i, we have |x i | σ, then the algorithm does not change the value of x i in that descent step; however, it might be changed in the projection step. If we are looking for a suitably large µ, a suitable choice is to make the algorithm force all those values of x satisfying |x i | σ toward zero. Therefore, we can get: and: Combining Equations (24) and (25), we can further obtain: By calculation, we can obtain: According to the above derivation, we have come to the conclusion that µ ≈ σ 2 2α . Therefore, we can set µ = σ 2 2α .

Selection of Parameter σ
According to Equation (17), the descending sequence of σ is generated by (17)). Parameter σ min and parameter σ max should be appropriately selected. The selection of σ min and σ max is discussed below. For the initial value of σ, i.e., σ max , here, letx = max{|x|}; suppose there is a constant b, in order to make the algorithm converge quickly; let parameter σ max satisfy: From the equation, we can see that constant b satisfies 1−b b ≥ 0; thus 0 < b ≤ 1, and here, we define constant b as 0.5. Hence, σ max = √ α max{|x|}. For the final value σ min , when σ min → 0, H σ min (x) → ||x|| 0 . That is, the smaller σ min , the more H σ min (x) can reflect the sparsity of signal x, but at the same time, it is also more sensitive to noise; therefore, the value σ min should not be too small. Combining [19], we choose σ min = 0.01.

Performance Simulation and Analysis
The numerical simulation platform is MATLAB 2017b, which is installed on a computer with a Windows 10, 64-bit operating system. The CPU of the simulation computer is the Intel (R) Core (TM) i5-3230M, and the frequency is 2.6 GHz. In this section, the performance of the WReSL0 algorithm is verified by signal and image recovery in the noise case.
Here, some state-of-the-art algorithms are selected for comparison. The parameters are selected to obtain the best performance for each algorithm: for the BPDNalgorithm [36], the regularization parameter λ = σ N 2log(n); for the SL0 algorithm [19], the initial value of smoothed factor δ max = 2 max{|x|}, the final value of smoothed factor δ min = 0.01, scale factor is set as step size L = 5, and the attenuation factor ρ = 0.8; for the NSL0algorithm [20], the initial value of smoothed factor δ max = 4 max{x}, the final value of smoothed factor δ min = 0.01, the step size L = 10, and the attenuation factor ρ = 0.8; for L p -RLSalgorithm [24], the number of iterations T = 80, the norm initial value p 1 = 1, the norm final value p T = 0.1, the initial value of regularization factor 1 = 1, the final value of regularization factor T = 0.01, and the algorithm termination threshold E t = 10 −25 ; for the WReSL0 algorithm, the initial value of smoothed factor σ max = √ c max{|x|}, the final value of smoothed factor σ min = 0.01, the iterations T = 30, the step size L = 5, and the regularization parameter λ = 0.1. All experiments are based on 100 trials.

Signal Recovery Performance in the Noise Case
In this part, we discuss signal recovery performance in the noise case. We add noise b to the measurement vector y; moreover, b = δ N Ω, Ω is randomly formed and follows the Gaussian distribution of N (0, 1). For signal recovery under noise conditions, we evaluate the performance of algorithms by the normalized mean squared error (NMSE) and the CPU running time (CRT). NMSE is defined as ||x − x|| 2 /||x|| 2 . CRT is measured with tic and toc. In order to analyze the de-noising performance of the WReSL0 algorithm in context closer to the real situation, we constructed a certain signal as an experimental object in the experiments in this section. The signal is given by: where α 1 = 0.2, α 2 = 0.1, β 1 = 0.3, and β 2 = 0.4. f 1 = 50 Hz; f 2 = 100 Hz; f 3 = 200 Hz; and f 4 = 300 Hz. Here, t is a sequence with t = [1, 2, 3, ..., n], and T s is sampling interval with the value of 1 f s . f s is the sampling frequency with the value of 800 Hz. The object that needs to be reconstructed can be expressed as: where x x x ∈ R n is a sparse signal in the frequency domain, and it is the Fourier transform expression of X X X , y y y ∈ R m . Here, let n = 128, m = 64. Moreover, Φ can be represented as Φ = ψ ψ ψϕ ϕ ϕ; here, ψ ψ ψ is a randn matrix generated by a Gaussian distribution, and ϕ ϕ ϕ is a sparse basis matrix generated by Fourier transform. Here, ϕ ϕ ϕ can be given by Fourier I n×n , and I n×n is a unit matrix. This target signal X X X is sparse in Fourier space; hence, the signal X X X can be recovered from given {y y y, Φ} by CS recovery methods. Figure 3 shows the signal recovery effect. Obviously, BPDN and SL0 do not perform well, while NSL0, L p -RLS and the proposed WReSL0 perform quite well. This verifies that the regularization mechanism has a good de-noising effect. Figure 4 shows the frequency spectrum of the recovered signal by the selected algorithms. The spectrum of the signal recovered by our proposed WReSL0 algorithm is almost the same as the original signal, while other algorithms fail to achieve this effect.     Table 2 shows the CRT of all algorithms. The n changes according to a given sequence [170, 220, 270, 320, 370, 420, 470, 520]. From the table, for any n, SL0 has the shortest computation time, followed by WReSL0, NSL0, and L p -RLS, and BPDN has the longest computation time. The BPDN algorithm is generally implemented by the quadratic programming method, and the computational complexity of this method is very high, thus resulting in a large increase in the overall computation time of the algorithm. Furthermore, in L p -RLS, the iterative process adopts the conjugate gradient method with high complexity, while NSL0 and WReSL0 do not. Compared with NSL0, WReSL0 is more prominent in the decrease of computation time. The performance of each algorithm under different noise intensities is shown in Figure 5. When δ N = 0, SL0 outperforms other algorithms, but with the increase of δ N , the effect of SL0 becomes worse and worse. This result further illustrates that the traditional constrained sparse recovery algorithm does not have the performance of anti-noising. For BPDN, NSL0, L p -RLS, and WReSL0, they all applied the regularization mechanism, and they are indeed superior to SL0 in the noise case. Therefore, the proposed WReSL0 in this paper has the best de-noising performance.

Image Recovery Performance in the Noise Case
Real images are considered to be approximately sparse under some proper basis, such as the DCT basis, DWT basis, etc. Here, we choose the DWT basis to recover these images. We compare the recovery performances based on the four real images in Figure 6: boat, Barbara, peppers, and Lena. The size of these images is 256 × 256; the compression ratio (CR; defined as m/n) is 0.5; and the noise δ N equals 0.01. We still choose SL0, BPDN, NSL0, and L p -RLS to make comparisons. For image recovery, the object of image processing is given by: Here, Y Y Y, X X X, B B B are matrices, and among these, Y Y Y, B B B ∈ R m×n , X X X ∈ R n×n . In order to meet the basic requirements of CS, we perform the following processing: To perform image recovery, we valuate it by the peak signal to noise ratio (PSNR) and the structural similarity index (SSIM). PSNR is defined as: PSNR = 10 log(255 2 /MSE) (33) where MSE = ||x − x|| 2 2 , and SSIM is defined as: Among these, µ p is the mean of image p, µ q is the mean of image q, σ p is the variance of image p, σ q is the variance of image q, and σ pq is the covariance between image p and image q. Parameters c 1 = z 1 L and c 2 = z 2 L, for which z 1 = 0.01, z 2 = 0.03, and L is the dynamic range of pixel values. The range of SSIM is [−1, 1], and when these two images are the same, SSIM equals one.   Figure 7 shows the recovery effect of boat and Barbara with noise intensity δ N = 0.01. For boat and Barbara, the recovered images by SL0 and BPDN have obvious water ripples, while recovered images by other algorithms have no such water ripples. Similarly, for peppers and Lena, the recovered images by SL0 and BPDN are blurred compared with the recovered images by other algorithms. The NSL0, L p -RLS, and WReSL0 algorithms are also effective at noisy image recovery. For the NSL0, L p -RLS, and WReSL0 algorithms, their recovery effects are very similar. In order to further analyze the advantages and disadvantages of the algorithms, we analyze the PSNR and SSIM of the images recovered by these algorithms, and the results are shown in Tables 3 and 4. By observation and analysis, L p -RLS performs better than NSL0, and at the same time, WReSL0 outperforms L p -RLS. Hence, the WReSL0 proposed by this paper is superior to the other selected algorithms in image processing.

Application in Underdetermined Blind Source Separation
The problem of UBSS stems from cocktail reception, which is shown in Figure 8. Suppose the source signal matrix S(t) = [s 1 (t), s 2 (t), ..., s m (t)] T , the mixed matrix (Sensors) A is m × n (m n) matrix, the Gaussian noise G(t) = [g 1 (t), g 2 (t), ..., g m (t)] T is generated by Gaussian distribution, and the observed mixed signal matrix X(t) = [x 1 (t), x 2 (t), ..., x n (t)] T ; therefore, the general mathematical models of UBSS can be summarized as: Source Signal S1 Source Signal S2

Source Signal S3
Sensor A1 Sensor A2 Recovered Signal S1 Recovered Signal S2 Recovered Signal S3 Figure 8. Schematic diagram of cocktail reception signal mixing.
In fact, each signal has L data collected; therefore, X ∈ R m×L , A ∈ R m×n (m n), S(t) ∈ R n×L , and G ∈ R m×L , and G can be represented as δ N W (W obeys N (0, 1)). The purpose of UBSS is to use the mixed signal matrix x(t) to estimate the sof the source signal matrix s(t). In fact, this is the process of solving the underdetermined linear system of equations (ULSE). For this problem, we can use the two-step method to solve it, which is shown in Figure 9.
Clustering CS Figure 9. Schematic diagram of two-step method for UBSS.

X(t) A(t) S(t)
From Figure 9, firstly, we get the mixed matrix by the clustering method and then use CS technology to separate the signal, so as to restore the original signal.

Solving the Mixed Matrix by the Potential Function Method
In this section, we choose the potential function method to solve the mixed matrix A. To verify the performance of the proposed WReSL0 algorithm better, we choose four simulated signals and four real images to organize experiments in this section.
Suppose there are four source signals, which are: where f 1 = 310 Hz, f 2 = 210 Hz, f 3 = 110 Hz, and f 4 = 10 Hz. The length of each source signal s i (i = 1, 2, 3, 4) is 1024, and the sample frequency is 1024 Hz. These four signals are shown in Figure 10.
The four source images are the classic standard test images: boat, Barbara, peppers, and Lena, which are in Figure 6.
Suppose there are two sensors that receive signals and another two sensors that receive images. Mixed matrices A and B are set as: By this mixed matrix and added Gaussian noise (δ N = 0.1), we can get the two mixed signals, which are shown in Figure 11, and the two mixed images, which are shown in Figure 12. Then, we can get the estimated mixed matrixÂ andB by clustering by the potential function method [37]. As shown in Figure 13, the potential function method can cluster well. By clustering, we get the estimated values of A and B, as follows:Â By calculation, the error of solving the mixed matrix is 64%. This error range is much smaller than the classical k-means and fuzzy c-means, thus laying a foundation for the reconstruction of compressed sensing.

Using CS to Separate Source Signals
The next problem is to get S(t) from known A(t) andX(t). Here, we solve this problem by CS. The solution process is similar to the image reconstruction process. The difference is that the sparse basis used here is the Fourier basis. Then, we apply the proposed RWeSL0 algorithm to this process. First, we transform the obtained x(t) into column vectors: Then, we use the Fourier (for the sparse signal) or DWT (for the image) basis for sparse representation and extend the matrix and the valuated mixed matrix to obtain the sensing matrix.
For this equation, ⊗ denotes the Kronecker product sign, Fourier(·) represents the Fourier transform, and DWT represents the discrete wavelet transform. Therefore, the CS-UBSS model can be described as:

=B(t)ΨΘ(t) + G(t)
= ΦΘ(t) + G(t) (41) where Θ is the Fourier transform or DWT of S(t), so Θ is a sparse signal. As for UBSS in the images, firstly, each image matrix needs to be transformed into a row vector, then the four row vectors form a matrix S(t). At the same time, the sparse basis in Equation (40) needs to be replaced by DWT. Then, we can recover the source signal by CS. In summary, the above can be described as the flowchart in Figure 14.

The Effect of the WReSL0 Algorithm Applied to UBSS
In this section, we evaluate the effect of the WReSL0 algorithm applied to UBSS by the separation of signals and spectrum analysis.
The effect of the separation of signals is shown in Figure 15: the source signals are well separated, and the separation signals and the original signals are very similar. Figure 16 displays the error between the original source signal and the recovered source signal. It indicates that the error between the original source signal and the recovered source signal is fairly small, and the WReSL0 algorithm can better deal with the problem of UBSS. In addition, We get the time-frequency diagram of the restored signal by short-time Fourier transform. Figure 17 is the time-frequency diagram. From this figure, we find that each signal has the same frequency as the original signal, and it also validates the rationality of the proposed algorithm for UBSS.

Performance Comparisons of the Selected Algorithms
Here, we use the SL0, NSL0, and L p -RLS algorithms and the classical shortest path method (SPM) [38] to make a comparison in different noise cases. In order to analyze the situation of signal recovery clearly, we apply average SNR (ASNR) (for the signal) and average peak SNR (APSNR) (for the image) to evaluate. Let the original source signal be s i and the recovered source signal beŝ i , so ANSR is defined as: and PSNR is defined as: where M and N are the width and height of the image.
The ASNR comparisons are shown in Table 5. From the table, we can see that ASNR attenuates sharply when δ N increases from 0.15-0.2. The reason is that the error of the valuated mixed matrixÂ increases obviously, which leads those CS recovery algorithms to perform poorly. In fact, from this table, our proposed RWeSL0 algorithm performs well when δ N is less than 0.15, and when δ N is greater than 0.15, the L p -RLS algorithm performs best, followed by our proposed RWeSL0 algorithm.
The APSNR comparisons are shown in Table 6. In this table, It is clear that APSNR is not high, and it drops greatly when δ N increases from 0.15-0.2. From Figure 18, we can see that these separated images seem to be enveloped in mist, which leads to a low APSNR. Therefore, we will try our best to improve this problem in the future.
In summary, the CS technique can be used in UBSS and performs well especially for the signal recovery. Our proposed WReSL0 algorithm can perform well in UBSS for the signal recoverywhen the noise is small; and regarding image recovery, we will develop this in the future.

Conclusions
In this paper, we propose the WReSL0 algorithm to recover the sparse signal from given {y, Φ} in the noise case. The WReSL0 algorithm is constructed under the GD method, in which the update process of x in the inner loop adopts the regularization mechanism to enhance the de-noising performance. As a key part of the WReSL0 algorithm, a weighted smoothed function W T H σ (x) is proposed to promote sparsity and provide the guarantee of robust and accurate signal recovery. Furthermore, We deduced the value of µ and the initial value σ max to ensure the optimization performance of the algorithm. Performance simulation experiments on both real signals and real images show that the proposed WReSL0 algorithm performs better than the L 1 or L p regularization methods and the classical L 0 regularization methods. Finally, we apply the proposed WReSL0 algorithm to solve the problem of UBSS and also make comparisons with the classical SPM, SL0, NSL0, and Lp-RLS algorithms. Experiments show that this algorithm has some advanced performance. In addition, we would also like to apply the the proposed algorithm to other CS applications such as the RPCA [39], SAR imaging [40], and other de-noising methods [41].
Author Contributions: All authors have made great contributions to the work. L.W., X.Y., H.Y., and J.X. conceived of and designed the experiments; X.Y. and H.Y. performed the experiments and analyzed the data; X.Y. gave insightful suggestions for the work; X.Y. and H.Y. wrote the paper.
Funding: This research received no external funding.