A Novel Underdetermined Blind Source Separation Algorithm of Frequency-Hopping Signals via Time-Frequency Analysis

To address the significant performance degradation of conventional underdetermined blind source separation algorithms for frequency-hopping (FH) signals under time-frequency (TF) overlapping conditions, this brief presents a novel three-stage scheme based on the TF distribution of FH signals. In the first stage, key parameters of the FH signal are estimated using a TF binary graph. In the second step, the initial mixing matrix is estimated for non-overlapping and overlapping carrier frequencies employing density peaks clustering and tensor decomposition methods, respectively. In the third step, the final mixing matrix, directions of arrival (DOA), and source signals are estimated using the expectation-maximization algorithm within the nonnegative matrix factorization model. Finally, different segments of the FH signals are spliced together based on the DOAs of different source signals. Comprehensive experimental results demonstrate the superior performance of the proposed algorithm compared to state-of-the-art algorithms. Even at a signal-to-noise ratio of 5 dB, the correlation coefficient of the estimated source signals can still reach 0.91.


I. INTRODUCTION
E LECTRONIC warfare has become one of the primary forms of modern warfare [1].Frequency-hopping (FH) communication is a communication method in which the carrier frequencies transmitted by two cooperating parties are randomly hopped within a common frequency set with the same hopping pattern [2].Due to advancements in antenna systems [3], [4] and its benefits, such as low interception probability and strong anti-jamming capabilities, FH communication is extensively employed in military communication systems [5], [6].
The separation and subsequent processing of intercepted mixed FH signals enable the acquisition of valuable battlefield information from opponents.Effective estimation of  [7], [8], [9] to resolve this issue.UBSS primarily comprises two steps: mixing matrix estimation (MME) and source signal recovery (SSR).
The MME algorithms based on single source point detection (SSPD) in UBSS have been successively proposed in [10], [11], [12].These algorithms rely on the assumption that source signals are sparse in the time-frequency (TF) domain.The core procedure involves detecting the TF points dominated by single source signal.However, when multiple source signals substantially overlap in the TF domain, the SSPD algorithm struggles to effectively detect the SSPs for these overlapped source signals, leading to significant estimation errors in both the quantity of source signals and steering vectors.In contrast to the SSPD algorithm, which leverages the sparsity of non-stationary signals, the MME algorithm based on tensor decomposition (TD) [13], [14], [15] primarily utilizes the independence of source signals.Nevertheless, ensuring the uniqueness of tensor decomposition imposes strict restrictions on the number of sources and mixed signals, which considerably limits the algorithm's practical applications.
In [16], SSR is transformed into an l p norm minimization problem, with the solution obtained via the Lagrange multiplier algorithm.When the number of source signals is equal to or less than that of mixed signals, the subspace projection algorithm [14] effectively performs the SSR task.However, the performance is not ideal under lower signalto-noise ratio (SNR) conditions.In [12], the k-sparse component analysis algorithm is proposed to identify the source signal subspace using singular value decomposition (SVD).Nonnegative matrix factorization (NMF) algorithms [17], [18] are frequently employed to address underdetermined speech separation problems but require sufficient prior information.
Influenced by the aforementioned investigations, we propose a three-stage scheme for addressing the UBSS of FH signals, leveraging the TF distribution of FH signals.Our algorithm highlights the following features: (1) It thoroughly utilizes the FH signal parameter information contained within the mixed signals, focusing only on the effective frequency bins in the last two stages, thereby enhancing robustness to noise and saving calculation time; (2) When estimating the mixing matrix, the density peak clustering (DPC) algorithm clusters only the TF points near a single carrier frequency at each time, estimating a unique cluster center that also serves as the steering vector, thus avoiding potential errors caused by overall clustering; (3) As we construct the tensor based on mixed signals in average overlapping frequency bin, the uniqueness condition of tensor decomposition is relaxed to some extent; (4) By inputting the effective frequency bins near carrier frequencies into the NMF model, the initial estimated mixing matrix provides proper initialization for the expectation-maximization (EM) algorithm, ultimately achieving the tasks of final MME, direction of arrival (DOA) estimation, and SSR. Notation: and Im(•) denote the operation of inverse, conjugate transpose, transpose, conjugate, pseudo-inverse, Khatri-Rao product, outer product, AND, taking modulus, taking phase angle, arccosine, diagnonalization, trace, taking determinant, maximum, minimum, mathematical expectation, taking real part and taking imaginary part, respectively.The function ξ(•) is defined as follow: when the input is greater than 0, the output is 1; otherwise, it is 0.

II. PROPOSED METHOD
This brief considers a uniform linear array (ULA) consisting of M sensors, with each sensor capable of receiving N FH signals, where N > M. The key feature of FH signals is that the carrier frequency remains constant during the same hop duration and subsequently hops to a new frequency at the following hop timing.As illustrated in Fig. 1, the proposed algorithm consists of three stages.
The mixed signal received by the sensor [15] can be written as: where τ m,i denotes the time delay for the ith FH signal to reach the mth sensor, and n m (t) is the Gaussian white noise of the mth sensor.τ m,i is mainly related to DOA θ i , speed of light c l , and sensor space d.TF distribution is applied to the mixed signal, we can write that where X(t, f ), S(t, f ), and N(t, f ) are the TF coefficients of mixed signal vector, FH signal vector and noise vector at TF point (t, f ), respectively.A = [a 1 , a 2 , . . ., a N ] ∈ C M×N is a vandermonde matrix at moment t, and the steering vector can be written as ] T .The mixing matrix A remains constant within each hop duration and changes as the carrier frequencies hop.Therefore, it is necessary to estimate the hop timings and subsequently divide the mixed signal into different hop durations.
In this brief, we assume that the carrier frequencies of most FH signals are non-overlapped and that the mixing matrix is row full rank.

A. FH Signal Parameters Estimation
In this brief, a single-sensor FH signal parameter estimation algorithm was employed.To achieve high TF resolution, the modified B-distribution (MBD) [19] was utilized to obtain the TF energy distribution ρ(t, f) ∈ R T×F of a single mixed signal.The TF binary graph will facilitate the extraction of the TF ridges from ρ(t, f) more effectively, as demonstrated: where α 1 is an energy threshold serving to select high-energy TF points.The step of taking the derivative in the frequency direction is intended to locate local energy peaks.Since each FH signal is not entirely overlapped, the frequency distance between nonzero TF points should be greater than α 2 , which is the carrier frequency interval.The values within each frequency bin in B(t, f) are accumulated, resulting in the formation of an effective frequency set F e , which is composed of frequency bins situated within the α 2 range near the local peak frequency bins.
The nonzero points in F e form the TF ridges, and subsequent operations are performed on these points.The number of nonzero points and their frequencies at each moment are recorded.The most frequently occurring number of nonzero points is identified as the estimated number of source signals N. Moments with a nonzero point count of N are extracted from F e , and their frequencies are subjected to clustering.The number of clusters and the clustering centers are the number of hop durations K and the carrier frequencies f , respectively.The moments when the clustering categories change are considered as starting points for forward searching.If a moment t satisfies it is considered to belong to the kth hop duration.By taking the average of the last moment of the kth hop duration and the initial moment of the (k + 1)th hop duration, the kth estimated hopping time tk can be obtained.

B. Initial Mixing Matrix Estimation
We have estimated hop timings t, carrier frequencies f , the number of FH signals N and hop durations K.These are the key elements for solving the UBSS in this brief.
Given that within a hop duration, the time resolution is not a concern, the short-time Fourier transform (STFT) is utilized in the MME and SSR algorithms due to its corresponding inverse transformation and lower computational complexity.
The kth split mixed signal x k is converted to the TF domain to obtain X k (t, f) ∈ C T k ×F k .In TF analysis, the utilize of window functions inevitably leads to the presence of main lobes and sidelobes in the resulting spectral estimates.In this brief, the STFT employs the Hamming window with a main lobe bandwidth of η ≈ 2f s γ , where f s and γ represent the sampling frequency and the length of the window, respectively.When the distance between two carrier frequencies is less than 2η, TF overlapping occurs.Based on this characteristic, carrier frequencies are classified into overlapped and non-overlapped carrier frequencies.Our algorithm merely calculates the frequency bins near the non-overlapped carrier frequency or the overlapped frequency bin each time.
We assume that, in the vicinity of the i 1 th non-overlapping carrier frequency, the TF points within [ fk (i 1 ) − η : fk (i 1 ) + η] frequency bins are dominated by a single source signal.For these SSPs, we can derive the following equation: where X k o (t, f ) and X k o+1 (t, f ) are the STFT coefficients of the oth and (o + 1)th mixed signals at the TF point (t, f ).θ i 1 and a i 1 (2) denote the DOA of the i 1 th FH signal and the second element of the i i th steering vector, respectively.By using (4), we can construct the real-valued ratio matrix The conventional SSPD algorithm estimates multiple cluster centers from all SSPs through clustering, while the DPC algorithm in this brief only needs to estimate a steering vector from the corresponding V at each time.The core of the DPC algorithm is to find those points that have a higher local density μ p and a relatively larger distance ε p from higher local density points as the cluster centers.These two concepts [20] can be explained by the following equations: where D pq denotes the distance between p and q.D α 3 is a distance threshold, which represents the minimum value greater than α 3 of all the distances of the ratio matrix V.This algorithm effectively separates TF point (t c , f c ) with the highest μ p ε p value.The ratio of this TF point ) is the estimated steering vector element âi 1 (2), and we can easily get the estimated steering vector âi 1 = [1, âi 1 (2), . . ., âi 1 (2) M−1 ].The DPC algorithm in this brief directly estimates the points with the most significant influence in the set as the output, which reduces the harmful impact of noise.
Due to the poor performance of SSPD algorithm in overlapped carrier frequencies, this brief proposes the TF domain tensor decomposition method to complete the estimation of overlapped steering vectors.
The presence of noise is temporarily ignored.We assume that there are Ṅ FH signals are overlapped, the average frequency fk a (i 2 ) bin of these overlapping carrier frequencies is applied to construct a three-dimensional tensor T ∈ C M×M×L .It is formed by stacking L mixed signal covariance matrices.The lth mixed signal covariance matrix can be defined as: where H ] ∈ C Ṅ× Ṅ are the overlapped source signal vector, overlapped mixing matrix and lth source signal covariance matrix, respectively.And T i,j,l = κ l (i, j).The Canonical Polyadic Decomposition (CPD) of T is given by where Q ∈ C L× Ṅ is composed of diagonal elements of L source signal covariance matrices, and represents the ith column of the matrix Q.The unfolding of , it can be expressed as: 1,:,: , T T 2,:,: , . . ., T T M,:,: .(8) The condition that the CPD of tensor is unique can be expressed as: where Ȧ represent the submatrix of Ȧ without the last row.
The singular value decomposition of T is denoted as T = U V H , where U ∈ C M 2 × Ṅ , ∈ C Ṅ× Ṅ , and V ∈ C L× Ṅ .There is a nonsingular matrix G ∈ C Ṅ× Ṅ that makes the following equations hold: where E and 0 refer to the identity matrix and zero matrix, respectively.Ȧ is the submatrix of Ȧ without the first row.If (9) holds, U 1 is a column full rank matrix, and the result of CP decomposition is unique.By combining the steering vectors estimated by DPC and CPD, the initial estimated mixing matrix Âk of the kth hop duration will be obtained.

C. Final Mixing Matrix Estimation, DOA Estimation, and Source Signal Recovery
In the NMF model, the TF energy distribution of each source signal is considered as the result of the sum of G components, each component is the product of two nonnegative matrices.The model [17] can be described as follows: where ) where x ∈ C T k ×F k ×M×M is the mixture covariance tensor.The EM algorithm is divided into two steps: conditional expectations of natural statistics and parameters update Ḃy iteratively updating the parameters, the objective function can be minimized.To overcome the EM algorithm's sensitivity to initialization, the initial mixing matrix Â is utilized as the initialization.Since the main information of FH signals is concentrated within the main lobe bandwidth of the carrier frequencies, only these frequency bins are used as input for the algorithm.Upon completing the final estimation of the mixing Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
When the mixing matrix and carrier frequencies are known, the following equation can be applied in DOA estimation.
The DOA of each source signal remains relatively constant within a short time frame.Consequently, source signals from adjacent hop durations with similar DOAs can be considered as belonging to the same FH signal.We utilize this criterion to sort the separated FH signals, effectively addressing the permutation ambiguity issue.

III. EXPERIMENTS AND ANALYSIS A. Implementation Issues
Each FH signal has a sampling frequency f s of 128 MHz, a sampling time of 0.6 ms, utilizes binary phase-shift keying (BPSK) as the modulation mode, and the carrier frequency interval is 125 kHz.Based on the fact that the carrier frequencies of the FH signal are randomly generated in the simulation, non-overlapping or overlapping conditions may randomly occur in the multiple Monte Carlo trials.The frequency and time resolution of the MBD are 31250 Hz and 256 samples, respectively.α 2 , γ , the number of fast Fourier transform (FFT) points, L, the iterations of the EM algorithm, and G are set to 125 kHz, 512, 512, 140, 100, and 2, respectively.All the simulations were run with AMD Ryzen 5 2600X Six-Core Processor (3.60 GHz, 16 GB memory).
To quantitatively evaluate the algorithms in this brief, the root mean square error (RMSE), the normalized mean square error (NMSE), and the correlation coefficient ψ are used as evaluation criteria [14].RMSE t , RMSE f , and RMSE θ represent the RMSE of hopping timings, carrier frequencies, and DOAs, respectively.

B. FH Signal Parameters Estimation Performance
In this subsection, we have verified the performance of the FH signal parameters estimation algorithm by setting different numbers of hop durations and source signals.A total of 1000 Monte Carlo trials were performed at each SNR (5 dB, 25 dB).FH signals with different numbers of hop durations and source signals were randomly generated.Furthermore, the estimated carrier frequency error RMSE f , estimated hop timing error RMSE t , and estimation success rate were measured.
The parameter α 1 serves as an energy threshold.When its value is too small, many noise TF points will be misidentified as nonzero TF points and erroneously counted in the TF binary graph.Conversely, when its value is too large, many TF ridge points will be recognized as zero TF points, leading to a decline in the algorithm's performance.The estimation success rate is a priority metric, as the success of the parameter estimation algorithm directly impacts the subsequent UBSS results.Fig. 2 clearly illustrates the effect of α 1 .When the value of α 1 is set between 0.45 and 0.55, the estimation success rate and estimation accuracy of carrier frequency under various conditions are relatively stable.The former are all above 99%, while the latter's RMSE f is also close to 0. Although the estimation error of hop timings fluctuates within this range, under the maximum error condition (M=3, N=5, K=3, 5 dB),  RMSE t is 39.6 samples, which is still less than the time resolution of the TF distribution.This shows that the algorithm is able to accurately estimate the parameters of the FH signals under different SNRs, which provides a solid foundation for the UBSS.

C. Underdetermined Blind Source Separation Performance
In this subsection, the performance of the proposed algorithm for MME, DOA estimation, and SSR is verified through comparisons with state-of-the-art (SOAT) algorithms.
The value of α 3 determines the average proportion of neighboring points in the vicinity of each point to the total number of points.In previous DPC algorithms, this value typically ranges from 0.01 to 0.02.However, in this brief, each DPC only needs to separate the single dominant point from a category of points, so the value of α 3 will increase.As can be seen from Fig. 3, the performance of the MME algorithm is optimal under various conditions when the value of α 3 ranges from 0.3 to 0.4.
Since the algorithm in [14], [15] has a strict limitation on the relationship between the number of mixed and source signals, all simulations in this subsection are performed under the condition (M=3, N=4, K=2).There are 1000 Monte Carlo trials for each SNR ([5 dB:40 dB], in 5 dB steps).
In Fig. 4(a) and Fig. 4(b), the accuracy of our algorithm outperforms the comparison algorithm in both the MME and the DOA estimation at each SNR.In comparison to the SSPD algorithm in [12]: firstly, our algorithm identifies overlapped frequency bins in advance and overcomes the dependence of the SSPD algorithm on sparsity through TF CPD; Secondly, the clustering algorithm in [12] clusters SSPs of all categories simultaneously.In contrast, we use the DPC algorithm to cluster TF points of a category each time to obtain the unique cluster center, which avoids unnecessary estimation errors.
In Fig. 4(c), the SSR performance of our algorithm has obvious advantages over the comparison algorithms at lower SNRs.At an SNR of 5 dB, ψ of our algorithm still reaches Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.0.91 compared to the highest value of 0.48 in the comparison algorithm.In the process of recovering source signals, the comparison algorithms need to traverse all TF points.In contrast, our algorithm merely calculates the frequency bins adjacent to the carrier frequencies based on the TF distribution of FH signals, which diminishes the negative effect of noise.Therefore, our algorithm can achieve such superior performance at relatively lower SNRs.However, when the SNR increases to the range of 35 dB to 40 dB, the performance of the SSR algorithm proposed in this brief does not have an advantage compared to the comparison algorithm.Since we only process the frequency bins within the main lobe bandwidth containing the primary information of the source signals during the SSR procedure, this approach can indeed improve performance at lower SNRs.However, it also sacrifices some source signal information within the sidelobes, and this drawback becomes more evident at higher SNRs.
In Table I, it can be observed that the proposed algorithm has the shortest simulation time compared to other algorithms.The primary reason for this is that our algorithm focuses on calculating frequency bins near the carrier frequencies, which saves a significant amount of time.

IV. CONCLUSION
In this brief, we have proposed a three-stage scheme for solving the UBSS of FH signals via TF analysis.Firstly, we employ the TF binary graph to estimate FH signal parameters, which are provided as key information for the next two stages.Secondly, we exploit a novel algorithm integrating DPC and CPD to address the mixing matrix estimation problem under the condition of TF overlapping.Thirdly, the EM algorithm with proper initialization can efficiently recover source signals from the effective frequency bins.In cases of lower SNR at

A
Novel Underdetermined Blind Source Separation Algorithm of Frequency-Hopping Signals via Time-Frequency Analysis Yifan Wang , Yibing Li , Qian Sun , and Yingsong Li , Senior Member, IEEE Abstract-To address the significant performance degradation of conventional underdetermined blind source separation algorithms for frequency-hopping (FH) signals under time-frequency (TF) overlapping conditions, this brief presents a novel threestage scheme based on the TF distribution of FH signals.In the first stage, key parameters of the FH signal are estimated using a TF binary graph.In the second step, the initial mixing matrix is estimated for non-overlapping and overlapping carrier frequencies employing density peaks clustering and tensor decomposition methods, respectively.In the third step, the final mixing matrix, directions of arrival (DOA), and source signals are estimated using the expectation-maximization algorithm within the nonnegative matrix factorization model.Finally, different segments of the FH signals are spliced together based on the DOAs of different source signals.Comprehensive experimental results demonstrate the superior performance of the proposed algorithm compared to state-of-the-art algorithms.Even at a signal-to-noise ratio of 5 dB, the correlation coefficient of the estimated source signals can still reach 0.91.Index Terms-Underdetermined blind source separation, frequency-hopping signal, time-frequency analysis, nonnegative matrix factorization.

Fig. 2 .
Fig. 2. The effect of α 1 settings on the performance of the parameters estimation.(a) RMSE f .(b) RMSE t .(c) Estimation success rate.

Fig. 3 .
Fig. 3.The effect of α 3 settings on the performance of the MME algorithm.(a) The mixing matirx performance.(b) The DOA estimation performance.

Fig. 4 .
Fig. 4. Underdetermined blind source separation performance.(a) The comparison of mixing matrix estimation performance.(b) The comparison of DOA estimation performance.(c) The comparison of source signal recovery performance.
GN×T k are left factor matrix and right factor matrix of |S(t, f)| 2 , respectively.The noise covariance tensor v ∈ R F k ×M×M is formed by stacking the noise variance matrices of each frequency bin.The objective function can be expressed as: the gth component of the ith source signal.W i ∈ R F k ×G and H i ∈ R G×T k are left factor matrix and right factor matrix of |S i (t, f)| 2 , respectively.W, H, A, v are parameters of EM algorithm, whereW = [W 1 , W 2 , . . ., W N ] ∈ R F k ×GN and H = [H 1 , H 2 , . . ., H N ] ∈ R

TABLE I THE
SIMULATION TIME COMPARISON OF DIFFERENT ALGORITHM