Single-Channel Blind Source Separation using Adaptive Mode Separation-Based Wavelet Transform and Density-Based Clustering with Sparse Reconstruction

In this paper, the problem of single-channel blind source separation (SCBSS) is addressed using a novel approach that combines the adaptive mode separation-based wavelet transform (AMSWT) and the density-based clustering with sparse reconstruction. The proposed method is performed in the time–frequency domain and in a reverberant environment. First, using the Fourier transform, the amplitude spectrum of the observed mixture signal is obtained. Then, using variational scaling and wavelet functions, the AMSWT is used to adaptively extract spectral intrinsic components (SICs). To obtain a better time–frequency resolution, the AMSWT is applied to each mode. Thus, the SCBSS problem is transformed into a non-underdetermined. Then, for each frequency bin, the density-based clustering, reformulated to the eigenvector clustering problem, is performed to estimate the mixing matrix. Finally, sparse reconstruction is used to reconstruct the sources. The proposed method is evaluated using objective measures of separation quality. A computational complexity evaluation based on time consumption is also performed. Simulation results show that the proposed method is very effective for solving the SCBSS problem and provides better separation performances than the reference methods. However, the proposed method is computationally expensive.

In the literature, the BSS methods are classified as being linear and nonlinear, instantaneous and convolutive, and overcomplete and underdetermined. The convolutive mixture model of BSS is an effective way to represent the speech signal mixing mechanism in a reverberant environment [7,30]. The BSS problem can be formulated either in the time domain or in the frequency domain. The BSS can be also treated in the joint time-frequency (TF) domain where computationally efficient BSS algorithms are available.
In most situations and for many practical uses, only one-channel recording of mixture signals is available. This particular instance of the underdetermined source separation problem called single-channel source separation (SCSS) has been the subject of many studies. To address the single-channel audio source separation problem, numerous strategies have been proposed in the literature [12]. In [14], the authors attempted to combine the maximum-likelihood estimation and nonnegative matrix factorization (NMF) based on the Itakura-Saito divergence measurement. The shorttime Fourier transform (STFT) representation of the observed single-channel signal has been subjected to the NMF-based approach in [39]. The method requires the use of extra training data. A combination of the empirical mode decomposition (EMD) and independent component analysis (ICA), as well as wavelet transformations, has been suggested in [31]. However, the wavelet transform needs some specified basis functions to represent a signal, and there is no rigorous mathematical theory underpinning the EMD or its improved algorithms [20]. The bark scale aligned wavelet packet decomposition has been introduced in [26] where the separation step has been performed using the Gaussian mixture model (GMM), which was employed before the Fourier transform. In [45], the authors proposed the variational mode decomposition (VMD) method to solve the SCBSS problem. The separation process is performed using joint approximate diagonalization based on the fourth-order cumulant matrices. In [36], the authors presented a novel method for SCBSS in a noisy environment. The method is based on selecting the TF units of signal presence and computing the mixture spectral amplitude. The separation process is performed using TF masking. In [25], an adaptive signal separation has been proposed. The method uses a time-varying parameter that adapts locally to instantaneous frequencies and a linear chirp (linear frequency modulation) to model the signal components. The single-channel technique has been explored for muscle artifact removal from multichannel EEG [6].
The classic TF representation is computed using the STFT. The STFT does not reflect the time-varying information. Moreover, it yields a time-frequency representation with only uniform time and frequency resolution. A new adaptive mode separation-based wavelet transform (AMSWT) has been proposed in [24] based on [10,16]. The AMSWT method involves solving a recursive optimization problem to adaptively extract spectral intrinsic components (SICs). The limited support of each spectral mode is implemented to establish the spectral boundaries for wavelet bank configuration. Then, the spectral boundaries of the created wavelet bank configuration are used to highlight the spectral information. The AMSWT strategy is fully adaptive in the sense that one does not require prior knowledge.
In [41], a new method to solve the underdetermined BSS problem for convolutive mixture has been proposed. The method operates in the time-frequency domain, and it combines the density-based grouping and sparse source reconstruction. The densitybased clustering is introduced to estimate the mixing matrix, which is converted to an eigenvector clustering issue. The rank-one structure of the local covariance matrices of the mixture TF vectors is first used to extract the eigenvectors. By combining weight clustering and density-based clustering, the eigenvectors are subsequently grouped and tweaked to provide an approximated mixing matrix. The source reconstruction is transformed into a l p -norm minimization using an iterative Lagrange multiplier method. The Lagrange multiplier used to solve optimization problems under constraints aims to enforce the constraint, while the quadratic penalty improves the convergence. In the iterative formula, both the primal and dual variables are updated iteratively.
In this paper, a new method to solve the SCBSS problem is proposed. The method combines the AMSWT [24] and density-based clustering with the sparse reconstruction method introduced in [41]. The method is performed in three stages: (i) The amplitude spectrum of the observed mixture signal is obtained using STFT. The convolution in the time domain can be approximated by a multiplication in the STFT domain. (ii) A better TF resolution is obtained using the variational scaling and wavelet functions, which are applied to the spectral intrinsic components (SICs) extracted adaptively using the AMSWT. By creating virtual multichannel signals of the TF representation, the underdetermined single-channel problem is transformed to a nonunderdetermined problem. (iii) For each TF representation and each frequency bin, the density-based clustering, which is converted to an eigenvector clustering problem, and the sparse reconstruction, which is converted to a minimization problem, are, respectively, performed to estimate the mixing matrix and sources.
The BSSeval toolbox [15] is used to evaluate the proposed method's performance. The evaluation is performed in terms of many criteria such as source-to-distortion ratio (SDR), source-to-artifact ratio (SAR) and source-to-interference ratio (SIR). The proposed method is compared to the variational mode decomposition (VMD) method [45], adaptive spectrum amplitude estimator and masking method [36] and the nonnegative tensor factorization of modulation spectrograms method [3].
The following sections make up the remaining content. The SCBSS problem formulation is presented in Sect. 2. The AMSWT method, the density-based clustering method and the source reconstruction are the main focus of Sect. 3. Simulation results are presented in Sect. 4. Finally, conclusions and discussions are given in Sect. 5.

Convolutive Mixture Model
., x M (t)] T be a vector of M observed sources obtained via the mixing of N independent sources s(t) [s 1 (t), .., s N (t)] T . The BSS problem aims to estimate the N sources from the M mixtures. The convolutive mixture occurs through the propagation of the sound through space and multiple paths caused by reflections from different objects, especially in rooms and closed environments. The convolutive mixture is modeled as follows: The matrix form is given as: where h ji denotes the impulse response from source i to sensor j, and H is an MxN matrix that contains the kth filter coefficients.
In most cases and for many practical purposes, only one-channel recording is accessible. Numerous studies have examined this instance known as single-channel source separation. In this case, M 1. The convolutive SCBSS in the time-frequency domain is described as follows: The conventional source separation techniques are ineffective in this scenario. The issue in SCBSS might be viewed as a single observation combined with numerous unidentified sources.

Single-Channel Blind Source Separation Method
The different steps of the proposed method for single-channel blind source separation are summarized by the flowchart shown in Fig. 1.
The spectrum of the observed signal is obtained by the STFT. The convolution in the time domain is transformed into a multiplication in the STFT domain. The AMSWT approach is used to obtain an optimal spectral mode separation. Thus, the SCBSS problem is transformed into a non-underdetermined problem by establishing virtual multichannel signals of the TF representation of the observed signals. Then, the M time-frequency representations of the mixture are divided into Q nonoverlapping blocks.  As a preprocessing step at the mixing matrix estimation stage, the TF representation of the observed signal is whitened for each frequency bin x d . The whitening process is performed using the eigenvector matrix U x and the eigenvalue matrix and it is expressed by The estimation of the mixing matrix is reformulated into an eigenvector clustering issue. The ambiguity of scaling is solved by rescaling the estimated mixing matrix by restricting the first row. The order of the reconstructed sources is aligned, by grouping the nearby source TF vectors based on their correlation, in terms of power ratio, to resolve the permutation ambiguity [10].
The post-processing stage involves de-whitening the predicted mixing matrix by H U x 1/2 x H. Then, the source reconstruction is reformulated into a sparse minimization problem, whose solution was achieved using an initialization-corrected iterative Lagrange multiplier approach.
Finally, the estimated sources are obtained in the TF domain, which are transformed into the time domain using the modified method proposed in [27].

Adaptive Mode Separation-Based Wavelet Transform
The STFT is a TF representation, which has an even bandwidth distribution across all frequency channels and suffers from the TF resolution limitation due to the fixed window size. The speech signal is substantially nonperiodic and nonstationary. Therefore, the use of the STFT will result in mistakes, particularly when complex transitory phenomena like voice mixing occur in the signal.
The AMSWT performs a time-frequency analysis using variational scaling and wavelet functions. The method is built on the alternating direction method of multiplier (ADMM) solver [19], which then defines a bank of variational scaling functions and wavelets depending on the spectral boundaries that have been defined. As a result, the approximation coefficients are obtained as the inner product of the analyzed signal x with variational scaling function. The inner product of the analyzed signal x with variational wavelets yields the detail coefficients, which are expressed as: and where x is the input signal.
In [24], under the amplitude-modulated and frequency-modulated (AM-FM) assumption, the intrinsic modes u(t) have distinguishable features in the frequency domain. Using the ADMM solver, the spectral modes can be adaptively obtained similarly to the intrinsic mode functions (IMFs) extraction, to estimate compact modes: where x(t) is the signal to be decomposed under the constraint that over all modes' summation should be equal to the input signal, δ(.) is the Dirac impulse and δ(t) + j π t * u k (t) denotes the original data and its Hilbert transform. The variables u k , ω k and k denote the modes, their central frequencies and the mode number, respectively.
The spectral segmentation boundary number can be determined empirically as follows: K min n ∈ Z + |n ≥ 2ρlnN (7) where N is the signal length and ρ is the scaling exponent determined by the detrended fluctuation analysis (DFA). According to [24], (6) is solved using a quadratic penalty term; the parameter λ denotes the Lagrangian multiplier used to render the problem unconstrained Therefore, u k is determined recursively as the data-fidelity constraint. The center frequencies ω n+1 k are updated as the center of gravity of the corresponding mode's power spectrum using the following equation As a result, rather than using a predefined wavelet bank, we create adaptive wavelet banks based on spectral modes and their corresponding center frequencies, which represent the intrinsic components.
Authors in [24] used the mode bandwidth and central frequencies to define the boundaries between each mode; however, in the literature, some authors just used the average of the two central frequencies as spectral boundary, which ignores the spectrum distribution.
Consider the kth mode with an average frequency ω k and a spectral bandwidth β k . Then, the boundary k between the kth mode and the (k + 1)th mode is given by where 0 0 and k π . The authors apply the same principle used in the production of both Littlewood-Paley and Meyer's wavelets [9] for variational scaling functions and wavelets based on spectral boundaries. ∅ k and ψ k are, respectively, defined by the following equation, with γ is the parameter that ensures no overlap between the two consecutive transitions. (12) and and β(ν) is an arbitrary function defined as follows: The adaptive mode separation-based wavelet transform algorithm is summarized as follows: Step 1: Time-frequency representation Input: Observed mixture.
• Using the Fourier transform, obtain the amplitude spectrum signal.
• Obtain the appropriate spectrum spectral modes (segments). Execute the first inner loop and the second inner loop to update u k according to (9), and update ω k according to (10), respectively • Compute the proper spectral boundaries using (11). Then, using (12) and (13), the bank of variational scaling functions and wavelets based on the spectral boundaries is defined. • Finally, using (4) and (5), respectively, apply variational scaling and wavelet functions to each mode to obtain the time-frequency distribution.
Output: time-frequency distribution of the observed mixture.

Density-Based Clustering
In [41], the authors introduced the eigenvector clustering as an alternative to estimate the mixing matrix. The eigenvector clustering is based on two factors, which are the local density ρ q and the minimum distance δ q that may be taken between the point q and any additional points with a higher density. They are given, respectively, by and δ q min where the region for each data point is defined by a cut-off distance τ c , and υ qk are the elements of the similarity matrix: From the eigenvectors A whose elements are a q , the similarity matrix V is generated as follows:υ qk a q − (a H q a k ) 2 F , where . F denotes Frobenius norm [29] expressed as follows: The eigenvector extraction is based on the local covariance matrix R X where h i is called the steering vector representing each direction of the mixing matrix. According to [41], there is at least one subblock indexed as q i for which the associated local covariance R X q i has roughly a rank-one structure. This condition is exploited in [16] where the authors applied the eigenvalue decomposition (EVD) to the local covariance matrix R X q , which results in the following equation: where U q and q denote the eigenvector matrix and eigenvalue matrix, respectively. The extracted vector denoted a q corresponds to the largest eigenvalue of q and also represents the first eigenvector in U q . To obtain the eigenvector matrix A [a 1 , . . . , a Q ], the eigenvector extraction is done subblock wisely.
According to [41], the global maximum in the density indexed as q * has a minimum distance δ q * defined as follows:δ q * max q,k 1,...,Q (υ qk ) if ρ q * max q 1,...,Q (ρ q )(20) The two components are multiplied together to provide the following score: To get γ q Q q 1 , the scores from (20) are applied to all of the subblocks. The obtained scores are then arranged in a decreasing order. As a result, the eigenvectors with the greatest N scores define the clusters, which are denoted by C [c 1 , ..., c N ].
As mentioned in [41], it would be difficult to cluster the eigenvectors using solely the density-based strategy described above. To address this issue, a weight clustering approach to further tune the projected clusters is used [42]. The procedure of the weighted eigenvector clustering can be summarized in three steps.
First, the eigenvector is weighted by a kernel function defined as follows: b qk e ω 2 where k 1, .., N and ω qk a q − a q H c k c k 2 F .
Then, the weighted covariance matrix is created as: Finally, the EVD is applied to the weighted covariance matrix R b k : As an updated of cluster c k where k 1, ..., N , the eigenvector that corresponds to the largest eigenvalue from the equation (24) is extracted.
The mixing matrix estimation algorithm is summarized as follows: Step 2 : Mixing matrix estimation Input : X which represents the TF representation of the observed signal whose elements x d .

For each block q ∈ Q do
• Calculate the local covariance matrix of R X q using R X • Construct the eigenvector matrix A using (19).

End
• Using the eigenvector matrix A, compute the similarity matrix defined by (17) For each block q ∈ Q do • Calculate the local density ρ q and the minimum distance δ q and the score γ q using (15), (16) and (21)

Source Reconstruction
In [41], the sparsity-based method has been introduced as an alternative to reconstruct the estimated source signal using a l p -norm-based minimization measurement. (The convergence is guaranteed for 0 < p < 1.) The method consists of converting the source reconstruction problem to a sparse reconstruction minimization problem. A designed iterative Lagrange multiplier approach with an appropriate initialization procedure is used to solve this minimization problem.
The source reconstruction is performed to find the sparsest term of s d . For this, the maximum posterior likelihood of s d is expressed as: where the complex-valued super-Gaussian distribution P s i,d is given by: where p and γ control the shape and variance of the probability function, denotes the gamma function and H represents the estimated mixing matrix. The problem returns to solve the equivalent optimization problem given as follows: The Lagrange multiplier method is used to solve the optimization problem. Hence, the problem is reformulated to an unconstrained optimization problem as follows: where α is the Lagrange multiplier. The implicit solution of the problem is given as follows: where −1 (s d ) The iterative scheme to obtain the solution s d is given as follows: The source reconstruction algorithm is summarized as follows: Step 3: Estimation of the TF representation of the sources Input: Time-frequency representation of the observed signal denoted X whose elements x d and estimated mixing matrix H For each frequency bin do • Initialize the sources as s is less than a given threshold. End Aware that s Output: time-frequency representation of the estimated sources.
For each frequency bin d, since the iterative method computes successive approximations to the solution of the problem, the stopping criterion minimizes the iterative absolute error. The tolerance or threshold of the stopping criteria is determined to guarantee the best algorithm performance without resulting in a high computing time.

Simulation Results
To investigate the effectiveness of the proposed method, numerical simulations have been performed in a reverberant environment. The TIMIT database [37] and NOIZEUS database [1] were used to build the speech dataset, which was chosen at random (available online). The sampling rate of the speech signals is f s 16kHz, and the speakers might be either female or male. Using the technique outlined in [17], the propagation environment is simulated as a reverberant room shown in Fig. 2.
The room impulse response from the source i to the sensor is illustrated in Fig. 3. By adjusting the reverberant time, a variety of convolutive mixed signals can be produced. It is crucial to evaluate the transmission duration of the signal decay to 60 dB to reflect the room reverberation.
As an illustration, let the three sources shown in Fig. 4a. The three sources are convolutedly mixed in the virtual room shown in Fig. 1 using the room impulse responses shown in Fig. 2. The observed single-channel signal is shown in Fig. 4b. Figure 4c shows a frame of 1024 sample length of the observed signal. The obtained modes are shown in Fig. 4d. For this example, the decomposition of the observed signal results  Fig. 4e. A comparison between the STFT of the estimated frame and the original frame of the observed signal is shown in Fig. 4f. The estimated sources are shown in Fig. 4g. A comparison between the TF representations of original sources and estimated sources is shown in Fig. 4h.
As can be seen, the estimated sources are highly similar to the original sources. The proposed method based on the AMSWT method and density-based clustering with sparse reconstruction provides an accurate estimate of the source signals and results in a spectral content located with high accuracy.
The BSSeval toolbox [15] is used to analyze the performance of the proposed approach. The estimated sources are expressed as s s target + e interf + e noise + e artif for the objective performance criteria measurement, where s target refers to the source signals, e interf stands for interference from other sources, e noise stands for distortion brought on by noise and e artif includes all other artifacts introduced by the separation algorithm.
The parameter p of the l p -norm-based minimization method can have a significant impact on source reconstruction performance [41]. Many tests have been performed using different values of p to assess its effect on the source-to-distortion ratio (SDR) using the given dataset. Table 1 displays the obtained SDRs for the parameter p varying from 0.1 to 0.9 by a step of 0.2. The duration of the signal presented as an example is 24 sec   As can be seen, the SDR marginally increases as p increases and reaches its maximum when p 0.7. The parameter p is set to 0.7 in the subsequent experiments. Changing the value of the parameter p to take advantage of the source sparsity proves that the sparse reconstruction based on l p -norm-based minimization method is very effective.
The estimated sources' performances are evaluated using the SDR, the source-toartifact ratio (SAR) and the source-to-interference ratio (SIR) criteria, and compared with the performances of the estimated sources obtained via the VMD method [45], adaptive spectrum amplitude estimator and masking method [36] and the nonnegative tensor factorization of modulation spectrograms method [3]. The Figure 5 shows the mean square error (MSE) between the original signal and the estimated sources obtained via the proposed method and reference methods. The comparison is performed for different reverberation conditions where the reverberation time is varied from 100 ms to 500 ms by a step of 50 ms. As observed, the proposed method provides the smallest MSE even in a highly reverberant environment. Figures 6, 7 and 8 show, respectively, the SDR, SAR and SIR obtained by the proposed method and reference methods for different reverberant times. As can be seen, the proposed method results in a better performance in terms of the three criteria compared to the VMD, adaptive spectrum amplitude estimator and masking and nonnegative tensor factorization of modulation spectrogram methods in a reverberant environment. The proposed method results in higher performance criteria even in a highly reverberant environment.
The proposed method has been compared to the reference methods in terms of time computing. In general, the computational complexity [5,8] is a measure of the execution time. The Fourier transform of a signal of length N has a computational complexity of O(NlogN ) [5,8]. Then, using variational scaling and wavelet functions, the AMSWT is introduced to adaptively extract spectral intrinsic components (SICs). The AMSWT method is built on the ADMM solver, with a computational complexity of O n 2 . The density-based clustering has a computational complexity of O n 3 . The     computational complexity required to compute the similarity matrix is O n 2 , and the sparse reconstruction used to reconstruct the estimated source has a computational complexity of O n 3 . The sparse methods are computationally expansive. The experiments have been carried out using a PC with a 2.4 GHz processor and 4 GB of RAM. The comparison has been performed for reverberation time conditions of 100 ms and 400 ms. The obtained results are shown in Fig. 9. As can be seen, the proposed method has a computational cost than the reference methods both in a weakly reverberant environment and in a highly reverberant environment. The SCBSS methods in the time-frequency domain are computationally expensive.

Conclusion
A new method to solve the SCBSS problem has been presented. The method combines the adaptive mode separation-based wavelet transform (AMSWT) with adaptive mode separation and the density-based clustering with sparse reconstruction. The SCBSS problem is transformed into a non-underdetermined. The method operates in the time-frequency domain and a reverberant environment. The proposed method has been tested on speech datasets constructed from TIMIT and NOIZEUS databases for various reverberation time conditions. Simulation experiments indicate that the proposed method results in the smallest MSE and the highest values of SIR, SAR and SDR compared to the reference methods. The simulation results demonstrate the effectiveness of the proposed method to solve the SCBSS problem even in a highly reverberant environment. In terms of computational complexity, the proposed method is expensive.