Underdetermined Blind Source Separation of Synchronous Orthogonal Frequency Hopping Signals Based on Single Source Points Detection

This paper considers the complex-valued mixing matrix estimation and direction-of-arrival (DOA) estimation of synchronous orthogonal frequency hopping (FH) signals in the underdetermined blind source separation (UBSS). A novel mixing matrix estimation algorithm is proposed by detecting single source points (SSPs) where only one source contributes its power. Firstly, the proposed algorithm distinguishes the SSPs by the comparison of the normalized coefficients of time frequency (TF) points, which is more effective than existing detection algorithms. Then, mixing matrix of FH signals can be estimated by the hierarchical clustering method. To sort synchronous orthogonal FH signals, a modified subspace projection method is presented to obtain the DOAs of FH. One superiority of this paper is that the estimation accuracy of the mixing matrix can be significantly improved by the proposed SSPs detection criteria. Another superiority of this paper is that synchronous orthogonal FH signals can be sorted in underdetermined condition. The experimental results demonstrate the efficiency of the two proposed algorithms.


Introduction
Frequency hopping (FH) signals have been widely used in wireless communications and radar systems due to their low probability of detection and interception [1][2][3][4][5]. To obtain useful information in communication reconnaissance field, it is an important task to sort multiple FH signals. Hop timing, frequency, and direction of arrival (DOA) are three important physical quantities for sorting FH signals. DOA is a vital parameter in radar and wireless communication applications [6][7][8][9][10][11][12][13][14]. Many methods for estimating FH parameters have been reported [15][16][17][18]. However, the methods in [15][16][17][18] can only estimate asynchronous FH signal parameters under overdetermined condition, in which the number of sensors is more than that of signals. To sort synchronous FH signals in underdetermined situation, in which the number of sensors is less than that of signals, Fu and Sha used underdetermined blind source separation (UBSS) method [19,20] based on time frequency (TF) analysis. In general, in UBSS situation, the major steps to sort synchronous FH signals are estimating mixing matrix and calculating DOAs [19,20]. Sha ignored the frequency effect on mixing matrix, so the results were not useful in practice. Fu modified time frequency ratio of mixtures algorithm to estimate frequencies and mixing matrix, and then DOA was obtained to sort synchronous FH signals. However, the mixing matrix estimation in Fu's method is not satisfactory, which leads to low estimation accuracy of DOA. Many other alternative methods are introduced, such as [21,22].
Many algorithms have been reported to deal with the UBSS problem; sparse component analysis (SCA) is an efficient method [21,22]. Based on SCA, UBSS can obtain favorable results by utilizing where n b is the th n signal amplitude, and ( ) n f t is the FH hopping instantaneous frequency.
Phase  ( ) n t is affected by frequency and modulation. T is the observation time. To avoid interference of irrelevant aspects, the amplitude and the phase can be ignored in this paper.
where ( ) m v t is the noise. The vector formulation of (3) is . Assume that the FH signals are not hopping in delay time, thus, according to (1) and (2), the transmission coefficient matrix ( ) t A is written as  The nth FH signal is written as where b n is the nth signal amplitude, and f n (t) is the FH hopping instantaneous frequency. Phase ϕ n (t) is affected by frequency and modulation. T is the observation time. To avoid interference of irrelevant aspects, the amplitude and the phase can be ignored in this paper. The delay time cannot be ignored when signals impinge on different antennas. The propagation time-delay from the mth antenna to the first antenna of the nth source can be formulated as follow according to Figure 1, where c denotes the speed of light, d denotes the element space and θ n denotes the nth direction of arrival (DOA). Therefore, the mth received signals can be written as where v m (t) is the noise. The vector formulation of (3) is where Assume that the FH signals are not hopping in delay time, thus, according to (1) and (2), the transmission coefficient matrix A(t) is written as where Ω m,n (t) = −2π f n (t)τ m,n . Obviously, according to (2), A(t) is a Vandermonde matrix. The data model in (4) is the same with blind source separation, so A(t) is the mixing matrix. The mixing matrix of FH signals changes with hopping frequencies and it is not easy to be estimated. In this paper, firstly the hop timings are obtained by other methods [40,41]. Then, the mixing matrix is constant in hop duration. The data model can be written in the kth hop duration where a n is the nth column of the mixing matrix. To identify FH signal slices, more parameters need to be estimated, such as DOAs and frequencies, which are included in the mixing matrix. Thus, the major task is to estimate the mixing matrix.

The UBSS Assumptions
Under UBSS condition, the number of sensors is less than the number of sources. Similar to other UBSS methods [19,20], the following two assumptions should be satisfied. First, any M × M sub-matrix of A is full rank so that all sources can be recovered. This assumption is easily satisfied due to the Vandermonde structure of mixing matrix. Second, suppose that only one source is active at some points in TF domain. FH signal is sparse in time-frequency domain. Single-source-points detection method is effective to deal with sparse signals. The method in [42] assumes that only one source is active at each time-frequency point. The algorithm in [27] requires that there exists at least one time-frequency point of single-source-occupancy for each source. The algorithms in References [35,36,43] assume that there exist some points in the TF domain where only a single source is active and such points must exist for each source. Therefore, the second assumption is reasonable. If the frequency hopping signal does not collide, this assumption is valid.

The Proposed Algorithms
In this section, firstly, the proposed SSPs detection algorithm is introduced to estimate mixing matrix. Secondly, the improved subspace projection algorithm is proposed to recover sources.

Single Source Detection
For simplicity, STFT is applied on both side of (6) without noise term. The mixtures in TF domain can be written as where X(t, f ) = [X 1 (t, f ), . . . , X M (t, f )] T are STFT coefficients of received signals, and S(t, f ) = [S 1 (t, f ), . . . , S N (t, f )] T are the STFT coefficients of sources.
Since the methods in [35,36] treat a large number of multi-source points as single-source points, the estimation results of mixing matrix are not satisfactory. Therefore, a new theorem is proposed to detect the single source points in TF domain. Theorem 1. At any TF point (t a , f a ) , if only one source s n (t) contributes its power, the following condition is required: where 0 < ε 1. ε is an empirical value and it depends on the signal ratio to noise (SNR). The effect of ε on estimation accuracy is discussed in Section 4.
Proof. The proof is given in Appendix A.

The Mixing Matrix Estimation
After detecting all SSPs, the mixing matrix can be estimated by clustering methods. Suppose that G n is a TF region which contains all SSPs of s n (t). If the second row of the mixing matrix can be obtained, the mixing matrix is obviously calculated based on the Vandermonde structure. According to (7), where a i,n represents the ith row of the nth column of mixing matrix. Because the first row elements of Vandermonde matrix are equal to one, the following result can be obtained To improve the mixing matrix estimation accuracy, based on (5) and (9), Ω 2,n which corresponds to the SSP (t a , f a ) ∈ G n , is represented as arctan Im{a 2,n } Re{a 2,n } + π if Re{a 2,n } < 0 and Im{a 2,n } Re{a 2,n } < 0 arctan Im{a 2,n } Re{a 2,n } others . (11) Cluster Ω 2,n into N different classes. The proposed algorithm utilizes the hierarchical clustering technique to estimate clustering centers Ω 2,n , n = 1, . . . , N . Although it may not be the best clustering algorithm, the clustering method will be studied in the future work.
Therefore, the mixing matrix estimation is calculated bŷ

FH Signals Recovery
An improved subspace projection method is proposed which is inspired by the method in [38]. The method in [38] assumed that there were no more than M active sources at each TF point. It calculated the minimum value of the product of the orthogonal projection matrix of mixing matrix and the TF domain mixtures to identify the sources. However, this method did not consider the effect of redundant signals. That is why the number of active sources at each TF point is always estimated as M. To improve the accuracy of the recovered FH signals, the proposed method in this paper sets a threshold to solve this problem. Suppose that there exist l (l ≤ M) FH signals at each TF point and l is a variable number. The corresponding mixing matrix column vectors are denoted by [a k 1 , . . . , a k l ], where, {k 1 , . . . , k l } ∈ {1, . . . , N}, k 1 < k 2 < · · · < k l . The l FH signals can be identified by where . . ., a k l . ε 0 is a small value and it is relative to redundant weak signals. The recovered FH sources in time domain can be estimated by inverse STFT.
The frequencies f n , n = 1, · · · , N can be calculated by fast Fourier transformation (FFT).

DOA Estimation and Splicing of FH Signals
For synchronous FH signals, hopping times of FH signals are the same. To splice FH signals, DOAs are required. According to the mixing matrix in (5), DOAs of FH signals in kth hop duration can be estimated byθ Sensors 2017, 17, 2074 6 of 20 DOAs do not change in a short time. If the DOAs of segments are approximately equal, these segments belong to the same FH signal. Therefore, comparing the DOAs between hop durations is an effective way to splice FH signals. The decision criterion is where DOA i and DOA j represent different segments, stands for the hypothesis of "not the same FH signals segments", H 0 stands for "the same FH signals segments" and r is a judgment threshold. The value of r will be described in Section 4.3.
In Summary, the proposed algorithm is given in Algorithm 1.

Algorithm 1 The proposed algorithm
Input: The received signals Step 1: Segment the received signals based on the known or estimated hop timings.
Step 2: In one hop duration, calculate the STFT detect signal-source-points.
Step 7: In another hop duration, repeat Step 2 to Step 6.

Simulations
In this section, numerical results are given to demonstrate the effectiveness of the proposed algorithms. In all experiments, M = 4. Suppose sample size L = 256 and N = 5. The synchronous FH signals parameters are shown in Table 1. In this section, the thresholds ε, ε 0 and r are set to 0.2, 0.1 and 2, respectively. The influence of thresholds on the estimation performance is discussed later. Hop timings can be estimated by the methods in [40,41]. This parameter is ignored in Table 1. In the following simulations, the SNR is written as Following [44], the performance of parameter estimation can be measured by Root mean square error (RMSE) which is defined as where η n is the nth original parameter, andη n is the nth estimation value.

Mixing Matrix Estimation
The performance of the mixing matrix estimation is demonstrated by mean square error (MSE) which is defined as A small MSE indicates superior quality. In addition, the columns of the estimated mixing matrix are normalized and the uncertainty of the columns is removed.
First, the underdetermined situation where M = 4 and N = 5 is considered. This experiment does 1000 Monte Carlo tests. Figure 2 illustrates the MSE performance of the proposed method comparing with Li's method in [35] and Nie's method in [36] with SNR increasing. It can be seen that the proposed method can obtain higher accuracy than other algorithms. Since the methods of Li and Nie treat a large number of multi-source points as single-source points, the estimation results of mixing matrix are not satisfactory. Figure 3a shows the scatter diagrams of real parts and imaginary parts of {a 2,n , n = 1, · · · , N} which are calculated by (10) when SNR = 15 dB. Figure 3a shows N clusters clearly. However, it is hard to identify clusters from Figure 3b,c, which are calculated by Nie's and Li's methods, respectively. Each cluster represents a 2,n , n = 1, . . . , N. Because the proposed method detects fewer error points than the other two methods, the proposed method works very well.  The imaginary part The imaginary part    Figure 4b,c, respectively, many disturbing points that lead to the bad results of the mixing matrix exist. {Ω 2,n , n = 1, · · · , N} calculated by (11) are shown in Figure 4a when SNR = 15 dB. The amplitudes represent the values of {Ω 2,n , n = 1, · · · , N}. Since the proposed method has eliminated many unreliable SSPs, the {Ω 2,n , n = 1, · · · , N} in Figure 4a clearly show N straight lines, which represent the clusters. Although, Nie's method and Li's method can also show straight lines in Figure 4b,c, respectively, many disturbing points that lead to the bad results of the mixing matrix exist. To analyze the complexity of these methods, CPU time is calculated. Our simulations are performed in MATLAB R2009a using an Inter Core I3-2100 +310 GHz processor. The operating system is Microsoft Windows XP and the memory is 2 GB. The average CPU time (s) for Monte Carlo tests against SNR levels is shown in Figure 5. The proposed method utilizes all mixtures in TF domain to detect SSPs, while Nie's method just uses two mixtures. That is why the proposed method cost more time than Nie's method. Because the numbers of SSPs are different in different To analyze the complexity of these methods, CPU time is calculated. Our simulations are performed in MATLAB R2009a using an Inter Core I3-2100 +310 GHz processor. The operating system is Microsoft Windows XP and the memory is 2 GB. The average CPU time (s) for Monte Carlo tests against SNR levels is shown in Figure 5. The proposed method utilizes all mixtures in TF domain to detect SSPs, while Nie's method just uses two mixtures. That is why the proposed method cost more time than Nie's method. Because the numbers of SSPs are different in different SNR, the cost time of proposed method is not the same. The cost time of the proposed method is increasing with the number of SSPs. To analyze the complexity of these methods, CPU time is calculated. Our simulations are performed in MATLAB R2009a using an Inter Core I3-2100 +310 GHz processor. The operating system is Microsoft Windows XP and the memory is 2 GB. The average CPU time (s) for Monte Carlo tests against SNR levels is shown in Figure 5. The proposed method utilizes all mixtures in TF domain to detect SSPs, while Nie's method just uses two mixtures. That is why the proposed method cost more time than Nie's method. Because the numbers of SSPs are different in different SNR, the cost time of proposed method is not the same. The cost time of the proposed method is increasing with the number of SSPs. is shown in Figure 6. It demonstrates that the performance of the proposed method is also Second, to analyze the source number effect on MSE performance, the underdetermined situation where M = 4 and N = 6 is considered. The comparison of the MSE for M = 4 and N = 6 is shown in Figure 6. It demonstrates that the performance of the proposed method is also better than the others when the source number increases to six. Comparing with Figure 2, it can be seen that the proposed method gets large MSE when N = 6. In other words, the performance becomes bad when the number of sources increases. Therefore, the number of FH signals cannot be much larger than that of antennas. better than the others when the source number increases to six. Comparing with Figure 2, it can be seen that the proposed method gets large MSE when  6 N . In other words, the performance becomes bad when the number of sources increases. Therefore, the number of FH signals cannot be much larger than that of antennas. Third, the influence of threshold  on the MSE of mixing matrix is discussed. When  is too small, not enough single-source-points exist to be used for estimating mixing matrix in low SNR. This will lead to wrong results. When  is too large, many more multi-source points exist to be used for estimating mixing matrix. This will also lead to wrong results. The effect of  is shown in Figure 7. Third, the influence of threshold ε on the MSE of mixing matrix is discussed. When ε is too small, not enough single-source-points exist to be used for estimating mixing matrix in low SNR.
This will lead to wrong results. When ε is too large, many more multi-source points exist to be used for estimating mixing matrix. This will also lead to wrong results. The effect of ε is shown in Figure 7. When ε = 0.1, the performance is bad in low SNR. When ε is greater than 0.3, the results are not credible. Through the above analysis, it is reasonable to set ε = 0.2. Third, the influence of threshold  on the MSE of mixing matrix is discussed. When  is too small, not enough single-source-points exist to be used for estimating mixing matrix in low SNR. This will lead to wrong results. When  is too large, many more multi-source points exist to be used for estimating mixing matrix. This will also lead to wrong results. The effect of  is shown in Figure 7.  The effect of sample size on the performance of mixing matrix is discussed in this part. Figure  8 shows that, when 256 L  , more accuracy can be obtained than However, if  is chosen too large, the accuracy is worse than 128 L  . Thus, in the case of large sample size, such as 256 L  , a small threshold  should be chosen, which leads to higher accuracy. The effect of sample size on the performance of mixing matrix is discussed in this part. Figure 8 shows that, when L = 256, more accuracy can be obtained than L = 128 if ε = 0.1. However, if ε is chosen too large, the accuracy is worse than L = 128. Thus, in the case of large sample size, such as L = 256, a small threshold ε should be chosen, which leads to higher accuracy.  Fourth, the influence of hop timing error on the MSE of mixing matrix is discussed. Figure 9 shows the impact of different hop timing errors on different sample sizes. When sample size is 128 L  , the performance of mixing matrix estimation becomes bad when hop timing error is greater than 30 samples. When sample size is 192 L  , the performance of mixing matrix estimation becomes bad when hop timing error is greater than 60 samples. When sample size is 256 L  , the performance of mixing matrix estimation is not bad even if hop timing error is greater than 100 samples. Thus, the greater the radio of the hop timing error to sample size is, the worse the result is. According to the existing hop timing estimation algorithms in [40,41], the hop timing can be estimated correctly in a certain range of error. Hop timing error is not very restricted to the proposed method. Fourth, the influence of hop timing error on the MSE of mixing matrix is discussed. Figure 9 shows the impact of different hop timing errors on different sample sizes. When sample size is L = 128, the performance of mixing matrix estimation becomes bad when hop timing error is greater than 30 samples. When sample size is L = 192, the performance of mixing matrix estimation becomes bad when hop timing error is greater than 60 samples. When sample size is L = 256, the performance of mixing matrix estimation is not bad even if hop timing error is greater than 100 samples. Thus, the greater the radio of the hop timing error to sample size is, the worse the result is. According to the existing hop timing estimation algorithms in [40,41], the hop timing can be estimated correctly in a certain range of error. Hop timing error is not very restricted to the proposed method.
Fourth, the influence of hop timing error on the MSE of mixing matrix is discussed. Figure 9 shows the impact of different hop timing errors on different sample sizes. When sample size is 128 L  , the performance of mixing matrix estimation becomes bad when hop timing error is greater than 30 samples. When sample size is 192 L  , the performance of mixing matrix estimation becomes bad when hop timing error is greater than 60 samples. When sample size is 256 L  , the performance of mixing matrix estimation is not bad even if hop timing error is greater than 100 samples. Thus, the greater the radio of the hop timing error to sample size is, the worse the result is. According to the existing hop timing estimation algorithms in [40,41], the hop timing can be estimated correctly in a certain range of error. Hop timing error is not very restricted to the proposed method.

FH Signals Recovery
The recovery performance is quantified by source to interference ratio (SIR), which is defined as

FH Signals Recovery
The recovery performance is quantified by source to interference ratio (SIR), which is defined as   First, Figure 10 compares the recovery performance between the improved subspace projection method and the method in [38]. These two methods utilize the mixing matrix that is estimated by the proposed SSPs detection method. It can be seen that the improved subspace projection method performs better than the method in [38]. Second, according to estimated mixing matrices that are obtained by the three methods, the improved subspace projection method is used to recover signals, and the performance is shown in Figure 11. It can be known that the higher the estimation accuracy of the mixing matrix is, the better (19) whereŝ n (t) is the nth recovered signal. A large SIR indicates superior quality.
First, Figure 10 compares the recovery performance between the improved subspace projection method and the method in [38]. These two methods utilize the mixing matrix that is estimated by the proposed SSPs detection method. It can be seen that the improved subspace projection method performs better than the method in [38].  First, Figure 10 compares the recovery performance between the improved subspace projection method and the method in [38]. These two methods utilize the mixing matrix that is estimated by the proposed SSPs detection method. It can be seen that the improved subspace projection method performs better than the method in [38]. Second, according to estimated mixing matrices that are obtained by the three methods, the improved subspace projection method is used to recover signals, and the performance is shown in Figure 11. It can be known that the higher the estimation accuracy of the mixing matrix is, the better the recovered performance is. This also proves that the proposed SSPs detection method is effective and can improve recovered signal performance. Second, according to estimated mixing matrices that are obtained by the three methods, the improved subspace projection method is used to recover signals, and the performance is shown in Figure 11. It can be known that the higher the estimation accuracy of the mixing matrix is, the better the recovered performance is. This also proves that the proposed SSPs detection method is effective and can improve recovered signal performance. Second, according to estimated mixing matrices that are obtained by the three methods, the improved subspace projection method is used to recover signals, and the performance is shown in Figure 11. It can be known that the higher the estimation accuracy of the mixing matrix is, the better the recovered performance is. This also proves that the proposed SSPs detection method is effective and can improve recovered signal performance. proposed method Li's method Nie's method Figure 11 The SIR comparison of different mixing matrix estimation methods. Figure 11. The SIR comparison of different mixing matrix estimation methods.
Third, frequencies can be estimated by FFT operator when the signals are obtained by the proposed recovery method in time domain. Figure 12 shows the comparison between the true and the estimated frequencies with SNR increasing. The square represents the estimated frequency and the straight line represents the true frequency. Figure 13 shows the RMSE of frequency estimation. Even when SNR is equal to 0 dB, the proposed methods can estimate the frequencies accurately. Third, frequencies can be estimated by FFT operator when the signals are obtained by the proposed recovery method in time domain. Figure 12 shows the comparison between the true and the estimated frequencies with SNR increasing. The square represents the estimated frequency and the straight line represents the true frequency. Figure 13 shows the RMSE of frequency estimation. Even when SNR is equal to 0 dB, the proposed methods can estimate the frequencies accurately.     Fourth, the influence of the threshold ε 0 on the SIR of the recovered FH signals is discussed. ε 0 effects the performance of waveform recovery. Figures 14 and 15 show that frequencies can be calculated by FFT even when using bad recovered waveforms. Thus, if the mixing matrix and frequencies are calculated, ε 0 has a minimal impact on the performance of DOA estimation.

DOA Estimation and Splicing FH Signals
The DOAs estimation can be obtained by (14) using the estimated mixing matrix and the frequency which are obtained by the proposed methods. Figure 16 shows the comparison of the true DOA and the estimated DOA with SNR increasing. The square represents the estimated value of DOA and the straight line means the true value of DOA. Figure 17 shows the RMSE of DOA estimation. It can be seen that the estimated DOAs are almost the same as the true ones when SNR is higher than 10 dB.  Figure 17 shows the RMSE of DOA estimation. It can be seen that the estimated DOAs are almost the same as the true ones when SNR is higher than 10 dB.     To sort FH signals, DOAs of different hop durations are estimated. Table 2 lists the RMSE performance of two hop durations when SNR is changing. Comparing the data in Table 2 with the data in Table 1, the DOA estimation is close to the true value when SNR is increasing. It can be seen that the RMSE of DOA estimation remains less than 1° when SNR is higher than 15 dB. In addition, the proposed algorithms can successfully sort synchronous orthogonal FH signals in the underdetermined case. The correlation of mixing matrix columns is related to the product of DOA and frequency. Under the assumption that any M M  sub-matrix of mixing matrix has full rank, the mixing matrix can be estimated even if DOA of FH signals are the same. Therefore, r can be selected based on the RMES of DOA estimation. r should be two times larger than that of RMSE of DOA estimation. According to Figure 17, the threshold r in (15) can be set to 2° so that different hops can be spliced. To sort FH signals, DOAs of different hop durations are estimated. Table 2 lists the RMSE performance of two hop durations when SNR is changing. Comparing the data in Table 2 with the data in Table 1, the DOA estimation is close to the true value when SNR is increasing. It can be seen that the RMSE of DOA estimation remains less than 1 • when SNR is higher than 15 dB. In addition, the proposed algorithms can successfully sort synchronous orthogonal FH signals in the underdetermined case. The correlation of mixing matrix columns is related to the product of DOA and frequency. Under the assumption that any M × M sub-matrix of mixing matrix has full rank, the mixing matrix can be estimated even if DOA of FH signals are the same. Therefore, r can be selected based on the RMES of DOA estimation. r should be two times larger than that of RMSE of DOA estimation. According to Figure 17, the threshold r in (15) can be set to 2 • so that different hops can be spliced.

Conclusions
This paper proposes a novel method to detect SSPs in TF domain and also an improved subspace projection method is proposed to recover signals. The synchronous FH sorting problem is solved in UBSS situation. The proposed SSPs detection method improves the complex-valued mixing matrix estimation accuracy. In this paper, DOAs and frequencies are also estimated by the improved subspace projection method. The RMSE of DOA estimation based on the proposed methods is less than 1 • when SNR is higher than 15 dB in the underdetermined case. The task of sorting FH signals has been completed by the proposed methods. We will try to improve DOA estimation accuracy in low SNR in the future work.
Utilizing the properties Re{a m,i } 2 + Im{a m,i } 2 = 1 and Re a m,j 2 + Im a m,j 2 = 1, (A15) is equivalent to the following Re{a m,i } = Re a m,j Im{a m,i } = Im a m,j . (A16) The condition in (A16) is hash in practice. Therefore, MSPs cannot satisfy (A10). The proof is completed.