A Parameter Estimation Algorithm for Multiple Frequency-Hopping Signals Based on Sparse Bayesian Method

Parameter estimation and network sorting for noncooperative wideband frequency-hopping (FH) signals have been essential and challenging tasks, especially in the case with little or even no prior information at all. In this paper, we propose a nearly blind estimation approach to estimate signal parameters based on sparse Bayesian reconstruction. Taking the sparsity in the spatial frequency domain of multiple FH signals into account, we propose a sparse Bayesian algorithm to estimate the spatial frequency parameters. As a result, the frequency and direction of arrival (DOA) parameters can be obtained. In order to improve the accuracy of the estimation parameters, we employmorphological filtermethods to further clean the data poisoned by the noise.Moreover, our method is applicable to the wideband signal models with little prior information.We also conduct extensive numerical simulations to verify the performance of our method. Notably, the proposed method works well even in low signal-to-noise ratio (SNR) environment.


Introduction
Frequency-hopping (FH) communication is widely used in the military communications, due to its high confidentiality, antijamming capacity, low susceptibility, and strong networking capabilities [1][2][3].However, these features also make it difficult to estimate the parameters of multiple FH signals, especially when there is no prior information about the hopping patterns.In general, the parameter estimation of FH signals often involves sorting and identifying of possible spatial FH signals belonging to different FH networks [2,4,5].Parameter estimation and network sorting are two fundamental but not well-solved problems in the field of FH signals reconnaissance.
Among various methods of parameter estimation, estimation methods based on maximum likelihood (ML) are practically intractable with few or no prior information, and thus several ML based methods have been proposed [6,7].Recently, another two techniques have attracted a lot of attentions, that is, time-frequency analysis and sparse reconstruction, since they can estimate the parameters by recovering clear time-frequency spectrum.Methods based on time-frequency analysis are attractive due to its simplicity and fast computation, but their performance degrades in low signal-to-noise ratio (SNR) scenarios [8,9].By contrast, sparse reconstruction based methods achieve good timefrequency spectrum and accurate parameter estimation in low SNR conditions [10][11][12].However, all these aforementioned methods suffer from degrading performance and mismatching of signals and parameters with multiple FH signals.
In order to achieve accurate estimation and network sorting results simultaneously, direction of arrival (DOA) information and blind separation (BS) is utilized with the assumption of multiple array elements.In many cases, FH signals are actually wideband, while many existing methods of DOA estimation consider the FH signals to be narrowband signals [13,14].To solve this issue, an expectation maximization (EM) approach is proposed for blind hop timing DOA and frequency estimation in the situation of possible bandwidth mismatch in [15]; the success of this approach mainly depends on the selection of the initial value.Unfortunately, it is difficult to obtain a reliable initial time-frequency spectrum using short-time Fourier transform (STFT) especially in low SNR environments.Paper [16,17] proposes a novel method for parameter estimation and With the assumption of wideband multiple FH signals, we propose a new method based on sparse reconstruction and morphological filtering to achieve better parameter estimation and network sorting results.We exploit the sparsity of the spatial frequencies of multiple FH signals to build redundant sparse dictionaries, based on which spatial frequencies are estimated for the reconstruction of FH signals and the DOA estimation can also be achieved.To improve the estimation accuracy in low SNR conditions, morphological filtering is used to achieve more accurate time-frequency and timespace spectrum.Therefore, parameters of FH signals can be estimated and then networks can be sorted with DOA.Moreover, the proposed method can be used for different multiple FH signals, different parameters, and different SNR s() = [ 0 (), . . .,  퐾−1 ()] conditions.
The rest of the paper is organized as follows.In signal model, the FH signal model and the sparse signal model are formulated.In proposed method, parameter estimation and network sorting algorithms are formulated mainly including sparse Bayesian reconstruction method and morphological filtering.In Algorithm Procedure, the procedure of the algorithm is summarized step by step.Several simulations and conclusions are given in Simulation and Analysis and Conclusion, respectively.

Signal Model
where  푘 denotes the amplitude of the th signal,  푘 () is the frequency of the th FH signal at the moment , and  푘 is the initial phase of the th FH signal.
The model after discrete sampling can be expressed as follows: where  푠 is the sampling frequency.Assume elements of the array are isotropic, ignoring the effect of mismatch and mutual of antenna.Assume the incidence angle of the th signal is  푘 , and its steering vector at th discrete moment can be described as Both frequency and incidence angle can be deduced from (3).Define  푘,푛 =  푘,푛  sin( 푘 )/ as "spatial frequency", where  푘 1 ,푛 ̸ =  푘 2 ,푛 , when  1 ,  2 ∈ {0, 1, . . .,  − 1}, and  1 ̸ =  2 , so (3) can be rewritten as Therefore, the signal vector at the th sampling moment can be expressed as follows: where whose spacing of division is Δ =  푙+1 −  푙 .The quantization error of each element in the spatial frequency corresponding to each signal is less than or equal to Δ/2, which makes the spatial frequencies at any moment constitute a small subset.In this way, (5) can be extended to the output model of redundant arrays as follows: where A(Φ)'s columns are steering vectors corresponding to each element in the sets Φ. s() are the extended version of vector s() with frequencies from Φ, which also extended from Φ 푛 .[s()] 푙 = [s()] 푘 ̸ = 0 if and only if  푙 =  푘,푛 .A(Φ) can be renamed as A for simplicity.
snapshots of received data are grouped as an operational frame.Data in two consecutive frames are not overlapped.The th frame can be expressed as where X 푔 = [x(), x( + 1), . . ., x( +  − 1)] 푇 and S 푔 = [s(), s( + 1), . . ., s( +  − 1)] and  denotes snapshots number.The position of nonzeros in s() is related to Φ 푛 .Φ usually have small sampling interval, whose number of samples satisfies  ≫  > .This means S is with great sparsity.In order to estimate parameters from (5), S 푔 should be bound to sparsity while obtaining best fitting of X 푔 to A S 푔 .Therefore, the optimization target function can be written as follows: where  is the penalty factor and  > 0. By the optimization above, the sparse solution S 푔 can be obtained, whose positions and amplitude of nonzero rows indicate the spatial frequencies of signals and signal amplitudes, respectively.

Proposed Method
The proposed method based on sparse Bayesian reconstruction is presented as the following steps.Firstly, signal frequency and DOA are roughly estimated by spatial frequencies obtained from reconstruction processing.Then, the estimated frequency and DOA are processed by morphological filtering, which produces precise parameter estimation and network sorting.

Spatial Frequency Estimation.
According to sparse Bayesian theory and ( 5), the probability density function of amplitudes from output data can be expressed as follows: where  2 푔 I 푀 denotes the variance of the th frame.To facilitate notation, intermediate variable  푔 = [ 1 (), . . .,  퐿 ()] 푇 is introduced to denote the power spectrum of incidence signal, and where Ω = diag ( 푔 ) and s() of different moments are mutually independent, and N(0, Ω) denotes Gaussian distribution with zero-mean variance of Ω.According to Bayesian Probabilistic Theory, the posterior probability density function (PDF) of S 푔 with respect to X 푔 can be expressed as where M S  and Σ S  denote the first and second moments of S 푔 .Given  푔 and  2 푔 , M S  and Σ S  can be estimated as From ( 12) and ( 13), it can be seen that the precision greatly depends on the prior distribution of  푔 and estimation of  2 푔 .Precise results can be obtained only when  푔 and  2 푔 can correspond to the energy of incidence signal.
To obtain  푔 more precisely, the observation data is used for optimization.The likelihood function of X 푔 with respect to  푔 is Taking logarithm of ( 14) and omitting constant terms, the target function with respect to  푔 is where R푔 = (1/)X 푔 X 푔 퐻 is the estimated array output covariance matrix.By minimizing (15), γ푔 , the power spectrum of incidence signal on spatial frequency sets Φ, can be obtained.
In practice, it is difficult to optimize directly on (15), so expectation maximization (EM) algorithm is used to estimate  푔 .Every iteration of EM includes an E-step and an Mstep.In E-step, the first-order moment  M S  and second-order moment  Σ S  are estimated by maximizing (11).In the Mstep, the  + 1 iteration of  푔 is estimated by calculating where ) 푙• denotes the th column of the matrix M S  after the th iteration and (  Σ ) 푙,푙 denotes the element on the th row and th column of the matrix Σ S after the th iteration, which can be deduced from ( 12) and (13).
To improve the speed of iteration and avoid exception caused by many zeros in  푔 , this iteration can be modified as where  is a small positive value.The estimated power 푔,푙 can be obtained after th iteration.The number of sources is known to be .The first  peaks ordered by amplitude are reserved for further processing.Every peak of the spectrum corresponding to every incidence signal usually contains multiple consecutive nonzero amplitudes because of quantization error of Φ.The spatial frequency can be roughly estimated via the linear interpolation of the main peak and the second largest amplitude beside the peak, as follows: where  (푞) 푘 () and where ), the notation (•) † denotes pseudo-inverse operation.When  푔 and  2 푔 meet the requirement for convergence, the iteration is terminated, and the power spectrum is estimated to be   푔 , which is taken into (12) to estimate the first-order moment  M S  .15) can be modified as follows in order to estimate spatial frequency precisely:

Precise Estimation of Spatial
where  푘 is the power of the th signal.The target function for signal power and spatial frequency is It can be inferred from ( 푘 ,  푘 )/ 푘 = 0, Taking (22) into ( 푘 ,  푘 )/ 푘 = 0, (23) can be deduced, where d( 푘 ) = a( 푘 )/ 푘 : The precise estimation of spatial frequency is as follows: In practical applications, the peak position can be obtained by searching in the smaller range [ 푘 (),  푘 ()], so that the spatial frequency can be precisely estimated.The computational complexity of the process above is moderate because it is conducted consecutively for  signals; the search range of which is only Δ.

DOA Estimation and Network
where angle( * ) means calculating the phase of a complex number.Therefore, the incidence angle can be estimated by  푘,푛 =  푘,푛  sin( 푘 )/ with signal frequencies: The incidence angle and frequency of  signals can be obtained from every frame.Considering the spatial incidence angle of the signal in the processing time is basically invariant, those angles are divided into  different clusters with known number of sources as .The frequencies corresponding to different spatial incidence angles are also categorized, so that different signals of different networks are sorted.

Parameter Estimation Based on Morphological Filtering.
In the frequency estimation of FH signals, one frame data may include information about two FH frequencies of the same signal; therefore, the estimation results are often in low credulity, since they suffered from great fluctuation in low SNR condition.
Many image filtering algorithms have been proposed to exclude these data points.Among these methods, median filter, wiener filter, and morphological filter are conventional methods to deal with binary image.Median filter is often used to suppress salt-and-pepper noise and it is a nonlinear filter.But the image we are processing has the feature that there is only one nonzero value in every time index; the median filter method is inefficient with the image of this type.Wiener filter is also a nonlinear filter.Under the criterion of minimum mean-square error, wiener filter is the best filter.Since the wiener filter needs iteration, it has a high computational complexity.morphological filtering method is used to improve the time-frequency spectrum since it has lower computational complexity and it is easier to be realized.
The obtained binary image is going through erosion and dilation operation.The scatter points of estimates can be eliminated by erosion, while the boundary can be restricted.The void in the signal can be filled by dilation, but the boundary would also be expanded.The erosion and dilation operations using structural element  on the image (, ) are expressed as where  푏 is the domain of .
There is only one nonzero value in each column of the image, and the probability of presence of consecutive points with greater fluctuation is small, so that the structural element is as [ 1 1].By erosion and dilation, outliers can be excluded while normal values are unaffected, generating a clarified time-frequency matrix   * 푘 (, ).The matrix is then turned back to the original estimates by coordinate transformation, with eroded scatter points modified by neighboring points.This transformation can be expressed as follows: where The modified time-frequency spectrum is turned into binary form  * 푘 (, ) and analyzed for component connectivity in order to get the hop moment, period, and frequency sets.The binarization is the same as (27).As the time-frequency points are consecutive by time, connected component detecting could be used to extract every hop of every signal.Eight connected components' detecting can get a better selection of the hop sequence than four connected components' detecting for its loose condition.The four connected components and the eight connected components could be seen in Figure 1.
The sets of points of each hop are  푘,0 , . . .,  푘,V  , . . .,  푘,푉  −1 , whose maximum and the minimum of time coordinates are  V  and  V  , indicating the moment of FH frequencies.The hop period estimate of the th signal is as follows: Therefore, the frequency estimate of V 푘 hop of th signal is The DOA estimate of th signal is

Algorithm Procedure and Computational Complexity
The procedure of joint parameter estimation and network sorting of multiple FH signals is summarized as follows. (
The computational complexity of the proposed method mainly depends on the heuristic and search algorithm.From the procedure of the algorithm, we can find that the computational complexity of one iteration mainly focuses on Steps 2 and 3. Therefore, it is in the order of ( 2 ), while the computational complexity of the algorithm in [11] is in the order of ( 2  2 ).It could be easily seen that the computational complexity of the proposed method is lower than the algorithm in [11].

Simulation and Analysis
In this section, the proposed method is analyzed by simulation to verify factors contributing to its performance.The simulated time-frequency spectrum and space-time pattern with SNR = 0 dB are depicted in Figures 2 and 3. Parameter estimates are listed in Tables 1 and 2.
As seen in Figure 2, the frequency estimates fluctuate greatly around the real value, corresponding to greater variance due to fewer number of snapshots in the sparse Bayesian method.By contrast, the estimates after morphological filtering are more accurate as shown in Figure 3.
As seen from Tables 1 and 2, the estimates of parameters of two sources are accurate.In the SNR of 0 dB, the error margin of DOA estimate is within 1 ∘ and the error margin of frequency sets is within 50 kHz, which indicates the proposed method is available for practical application.Figures 4, 5, and 6 show the estimation bias under frequency, DOA, and period in different SNR conditions of the proposed method and the particle filtering method in [14].The bias of signal frequencies, DOAs, and hop period are defined as follows: where  (푖) (푘,V  ) and   (푖) (푘,V  ) denote the frequency and estimate of the hop of the th signal in the th experiment;  (푖)  푘 and   (푖)  푘 denote the DOA and estimate of the th signal in the th experiment. (푖)  푘 and   (푖) 푘 denote the period and estimate of th signal in the th experiment. denotes the total number of experiments.
From Figures 4, 5, and 6, it can be seen that when SNR is greater than 5 dB, the bias of frequency and DOA of these two methods are similar, while when SNR is smaller than 5 dB, the proposed method shows much better performance.And the bias of period estimation of these two methods is nearly the same in all SNRs.

2. 1 .
FH Signal Model.Assume that  FH signals impinge onto a uniform linear array (ULA) with interantenna spacing ; the DOA is Θ = [ 0 ,  2 , . . .,  퐾−1 ], and the bandwidth of receiver is  =  max −  min , where  < ,  < /2 max , and  max and  min are the maximum and minimum frequency of the signals, respectively.The th FH signal can be expressed as follows: ) denote the spatial frequency of the largest and the second largest value of the th peak in the spectrum   (푞) 푔 , whose power is   ,2 , respectively.Based on the rough estimation above, the spatial frequency estimations  Φ )] in the th iteration can be obtained.The unbiased estimation of corresponding noise is

Figure 1 :
Figure 1: Four connected components and eight connected components.

( 4 )
), and Step 6 to obtain the time-frequency spectrum and space-time spectrum.Put space-time spectrums into different categories by different number of signals.Frequency and DOA of signal are separated based on the result of categorization.Parameter Estimation of FH Signals Based on Morphological Filtering

Experiment 1 .
Consider an ULA with eight elements and interelement spacing of 2.5 m.The bandwidth of receiver is [30 MHz, 60 MHz].Two far-field FH signals with incidence angle of  1 = −38 ∘ and  2 = 44 ∘ are assumed.The hop periods of two signals are 100 us and 50 us, respectively.The frequency set of the first signal is [42 MHz, 35 MHz, 53.5 MHz, 48.7 MHz], with [47 MHz, 53 MHz, 34.3 MHz, 42.5 MHz, 49 MHz, 56 MHz, 39 MHz] being the counterpart for the second signal.The total time length of signal procession is 300 us.Two signals are sampled with the rate of 150 MHz, and every operational frame has 100 snapshots.

Figure 2 :
Figure 2: Time-frequency spectrum and space-time spectrum after reconstruction.

2
Mathematical Problems in Engineering network sorting of wideband FH signals, which can only be used for multiple orthogonal signals but fails in asynchronous signals.

(𝑛)
, . . ., a 휙 , , . . ., a 휙 −1, ] is the array manifold array of the th sampling moment and s Sorting.The spatial frequencies of  signals can be obtained by algorithms listed in Sections 3.1 and 3.2. amplitudes  s 1 ,  s 2 , . . .,  s 퐾 can be obtained from  peak positions from the estimated firstorder moment  M S  , where  s 1 ,  s 2 , . . .,  s 퐾 are 1 ×  signal vectors.The number of samples is small, so phase-difference method is used to estimate frequency for  signals, with the following formula: −3is met, and output   푔 and  M S  .

Table 1 :
Estimate of DOA and period.Experiment 2. Consider an ULA with eight elements and interelement spacing of 2.5 m.The bandwidth of receiver is [30 MHz, 60 MHz].FH signals are generated with incidence angles randomly selected in (−90 ∘ , 90 ∘ ) and frequencies randomly selected in [30 MHz, 60 MHz].The frequency interval between adjacent two points is greater than 1 MHz.Assume the two frequencies meet | 1,푛 −  2,푛 | > 0.05 at any moment.The periods of two signals are 100 us and 50 us, respectively.The total time length of signal procession is 300 us.To get a convincing result, the signal number  is set to be 2 and 3, respectively.And signals are sampled with the rate of 150 MHz, and every frame has 100 snapshots.In every SNR condition from −5 dB to 10 dB by 1 dB, 100 Monte-Carlo simulations are performed.