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

Authors: Kun-feng Zhang; Ying Guo; Zisen Qi
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

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 employ morphological filter methods 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.


---

## Body

## 1. 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–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 time-frequency spectrum and accurate parameter estimation in low SNR conditions [10–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 narrow-band 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 network sorting of wideband FH signals, which can only be used for multiple orthogonal signals but fails in asynchronous signals.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 time-space 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 SNRst=s0t,…,sK-1t 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.


## 2. Signal Model

### 2.1. FH Signal Model

Assume thatK FH signals impinge onto a uniform linear array (ULA) with interantenna spacing d; the DOA is Θ=[θ0,θ2,…,θK-1], and the bandwidth of receiver is B=fmax-fmin, where K<M, d<c/2fmax, and fmax and fmin are the maximum and minimum frequency of the signals, respectively. The kth FH signal can be expressed as follows:(1)skt=bkej2πfktt+φk,where bk denotes the amplitude of the kth signal, fk(t) is the frequency of the kth FH signal at the moment t, and φk is the initial phase of the kth FH signal.The model after discrete sampling can be expressed as follows:(2)skn=bkej2πfk,nn/fs+φk,where fs 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 thekth signal is θk, and its steering vector at nth discrete moment can be described as(3)afk,n,θk=1,ej2πfk,ndsin⁡θk/c,…,ejM-12πfk,ndsin⁡θk/cT.Both frequency and incidence angle can be deduced from (3). Define ϕk,n=fk,ndsin⁡(θk)/c as “spatial frequency”, where ϕk1,n≠ϕk2,n, when k1,k2∈0,1,…,K-1, and k1≠k2, so (3) can be rewritten as(4)aϕn,k=1,ej2πϕk,n,…,ejM-12πϕk,nT.Therefore, the signal vector at thenth sampling moment can be expressed as follows:(5)xn=AΦnsn+nn,where A(Φn)=[aϕ0,n,…,aϕk,n,…,aϕK-1,n] is the array manifold array of the nth sampling moment and s(n)=[s0(n),…,sk(n),…,sK-1(n)] denotes incoming spatial signals with incidence angle vector Θ=[θ0,…,θk,…,θK-1] and spatial frequency vector Φn=[ϕ0,n,…,ϕk,n,…,ϕK-1,n]. The steering vector A(Φn) changes with the signal frequency of FH signal. n(n) denotes the Gaussian White Noise vector with zero mean and variance σ2IM, where IM denotes identity matrix.


### 2.2. Spatial Frequency Estimation Based on Sparsity

In this paper, the interelement spacing meetd<c/2fmax, and the incidence angle of signal is -90°<θk<90°, so that the range of spatial frequency ϕk,n is (−0.5, 0.5). If the sets of FH frequencies are finite, the possible combination of spatial frequencies are also finite, which satisfy the sparse condition in the spatial frequency domain. Spatial frequencies can be divided evenly into a complete discrete sets Φ=[ϕ1,…,ϕL], whose spacing of division is Δϕ=ϕl+1-ϕl. 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:(6)xn=AΦs-n+nn,where A(Φ)’s columns are steering vectors corresponding to each element in the sets Φ. s-(n) are the extended version of vector s(n) with frequencies from Φ, which also extended from Φn. [s-(n)]l=[s(n)]k≠0 if and only if ϕl=ϕk,n. A(Φ) can be renamed as A- for simplicity.N snapshots of received data are grouped as an operational frame. Data in two consecutive frames are not overlapped. The gth frame can be expressed as(7)Xg=A-S-g+N-,where Xg=[x(gN),x(gN+1),…,x(gN+N-1)]T and S-g=[s-(gN),s-(gN+1),…,s-(gN+N-1)] and N denotes snapshots number. The position of nonzeros in s-(n) is related to Φn.Φ usually have small sampling interval, whose number of samples satisfies L≫M>K. This means S- is with great sparsity. In order to estimate parameters from (5), S-g should be bound to sparsity while obtaining best fitting of Xg to A-S-g. Therefore, the optimization target function can be written as follows: (8)LS-∣Xg=Xg-A-S-g22+λS-gp,2,where λ is the penalty factor and λ>0.By the optimization above, the sparse solutionS-g can be obtained, whose positions and amplitude of nonzero rows indicate the spatial frequencies of signals and signal amplitudes, respectively.


## 2.1. FH Signal Model

Assume thatK FH signals impinge onto a uniform linear array (ULA) with interantenna spacing d; the DOA is Θ=[θ0,θ2,…,θK-1], and the bandwidth of receiver is B=fmax-fmin, where K<M, d<c/2fmax, and fmax and fmin are the maximum and minimum frequency of the signals, respectively. The kth FH signal can be expressed as follows:(1)skt=bkej2πfktt+φk,where bk denotes the amplitude of the kth signal, fk(t) is the frequency of the kth FH signal at the moment t, and φk is the initial phase of the kth FH signal.The model after discrete sampling can be expressed as follows:(2)skn=bkej2πfk,nn/fs+φk,where fs 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 thekth signal is θk, and its steering vector at nth discrete moment can be described as(3)afk,n,θk=1,ej2πfk,ndsin⁡θk/c,…,ejM-12πfk,ndsin⁡θk/cT.Both frequency and incidence angle can be deduced from (3). Define ϕk,n=fk,ndsin⁡(θk)/c as “spatial frequency”, where ϕk1,n≠ϕk2,n, when k1,k2∈0,1,…,K-1, and k1≠k2, so (3) can be rewritten as(4)aϕn,k=1,ej2πϕk,n,…,ejM-12πϕk,nT.Therefore, the signal vector at thenth sampling moment can be expressed as follows:(5)xn=AΦnsn+nn,where A(Φn)=[aϕ0,n,…,aϕk,n,…,aϕK-1,n] is the array manifold array of the nth sampling moment and s(n)=[s0(n),…,sk(n),…,sK-1(n)] denotes incoming spatial signals with incidence angle vector Θ=[θ0,…,θk,…,θK-1] and spatial frequency vector Φn=[ϕ0,n,…,ϕk,n,…,ϕK-1,n]. The steering vector A(Φn) changes with the signal frequency of FH signal. n(n) denotes the Gaussian White Noise vector with zero mean and variance σ2IM, where IM denotes identity matrix.


## 2.2. Spatial Frequency Estimation Based on Sparsity

In this paper, the interelement spacing meetd<c/2fmax, and the incidence angle of signal is -90°<θk<90°, so that the range of spatial frequency ϕk,n is (−0.5, 0.5). If the sets of FH frequencies are finite, the possible combination of spatial frequencies are also finite, which satisfy the sparse condition in the spatial frequency domain. Spatial frequencies can be divided evenly into a complete discrete sets Φ=[ϕ1,…,ϕL], whose spacing of division is Δϕ=ϕl+1-ϕl. 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:(6)xn=AΦs-n+nn,where A(Φ)’s columns are steering vectors corresponding to each element in the sets Φ. s-(n) are the extended version of vector s(n) with frequencies from Φ, which also extended from Φn. [s-(n)]l=[s(n)]k≠0 if and only if ϕl=ϕk,n. A(Φ) can be renamed as A- for simplicity.N snapshots of received data are grouped as an operational frame. Data in two consecutive frames are not overlapped. The gth frame can be expressed as(7)Xg=A-S-g+N-,where Xg=[x(gN),x(gN+1),…,x(gN+N-1)]T and S-g=[s-(gN),s-(gN+1),…,s-(gN+N-1)] and N denotes snapshots number. The position of nonzeros in s-(n) is related to Φn.Φ usually have small sampling interval, whose number of samples satisfies L≫M>K. This means S- is with great sparsity. In order to estimate parameters from (5), S-g should be bound to sparsity while obtaining best fitting of Xg to A-S-g. Therefore, the optimization target function can be written as follows: (8)LS-∣Xg=Xg-A-S-g22+λS-gp,2,where λ is the penalty factor and λ>0.By the optimization above, the sparse solutionS-g can be obtained, whose positions and amplitude of nonzero rows indicate the spatial frequencies of signals and signal amplitudes, respectively.


## 3. 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.


### 3.1. Spatial Frequency Estimation

According to sparse Bayesian theory and (5), the probability density function of amplitudes from output data can be expressed as follows:(9)pXg∣S-g=πσg2IMexp-σg-2Xg-A-S-g22,where σg2IM denotes the variance of the gth frame.To facilitate notation, intermediate variableγg=[γ1(g),…,γL(g)]T is introduced to denote the power spectrum of incidence signal, and (10)s-n~N0,Ω,n=1,…,N,where Ω=diag⁡(γg) and s-(n) 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-g with respect to Xg can be expressed as(11)pS-g∣Xg;γg,σg2=pXg∣S-g;σg2pS-g∣γg∫pXg∣S-g;σg2pS-g∣γgdS-g=πΣS-g-Nexp-trS-g-MS-gHΣS-g-1S-g-MS-g,where MS-g and ΣS-g denote the first and second moments of S-g. Given γg and σg2, MS-g and ΣS-g can be estimated as(12)M⌢S-g=ΩA-Hσg2I+A-ΩA-H-1Xp,(13)Σ⌢S-g=Ω-ΩA-Hσg2I+A-ΩA-H-1A-Ω.From (12) and (13), it can be seen that the precision greatly depends on the prior distribution of γg and estimation of σg2. Precise results can be obtained only when γg and σg2 can correspond to the energy of incidence signal.To obtainγg more precisely, the observation data is used for optimization. The likelihood function of Xg with respect to γg is(14)pS-gXg;γg,σg2=∫pXgS-g;σg2pS-gγgdS-g=πΣS-g-Nexp-trXgHσg2I+A-ΩA-H-1Xg.Taking logarithm of (14) and omitting constant terms, the target function with respect to γg is(15)κγg,σg2=lnσg2I+A-ΩA-H+trσg2I+A-ΩA-H-1R^g,where R^g=(1/N)XgXgH is the estimated array output covariance matrix. By minimizing (15), γ^g, 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 γg. Every iteration of EM includes an E-step and an M-step. In E-step, the first-order moment M⌢S-g and second-order moment Σ⌢S-g are estimated by maximizing (11). In the M-step, the q+1 iteration of γg is estimated by calculating ∂κ(γg,σg2)/∂γg=0 given σ2:(16)γ⌢g,lq+1=M⌢S-gql•22N+Σ⌢S-gql,l,where M⌢S-gql∙ denotes the lth column of the matrix MS-g after the qth iteration and Σ⌢S-gql,l denotes the element on the lth row and lth column of the matrix ΣS- after the qth iteration, which can be deduced from (12) and (13).To improve the speed of iteration and avoid exception caused by many zeros inγg, this iteration can be modified as(17)γ⌢g,lq+1=M⌢S-gql•22N1-Σ⌢S-gql,l/γ⌢g,lq+ξ,where ξ is a small positive value.The estimated powerγ⌢g,l(q) can be obtained after qth iteration. The number of sources is known to be K. The first K 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:(18)ϕ⌢kqg=γ⌢g,k,1qγ⌢g,k,1q+γ⌢g,k,2qϕ_kqg+γ⌢g,k,2qγ⌢g,k,1q+γ⌢g,k,2qϕ-kqg,k=1,…,K,where ϕ_k(q)(g) and ϕ-k(q)(g) denote the spatial frequency of the largest and the second largest value of the kth peak in the spectrum γ⌢g(q), whose power is γ⌢g,k,1(q) and γ⌢g,k,2(q), respectively.Based on the rough estimation above, the spatial frequency estimationsΦ⌢(q)=ϕ⌢1qg,…,ϕ⌢Kqg in the qth iteration can be obtained. The unbiased estimation of corresponding noise is(19)σ⌢g2Φ⌢q=PAΦ⌢q⊥Xg22NM-K,where PA(Φ⌢(q))⊥=I-A(Φ⌢(q))A†(F⌢(q)), the notation (•)† denotes pseudo-inverse operation. When γg and σg2 meet the requirement for convergence, the iteration is terminated, and the power spectrum is estimated to be γ⌢g, which is taken into (12) to estimate the first-order moment M⌢S-g.


### 3.2. Precise Estimation of Spatial Frequency

The reconstruction in Section3.1 would lead to estimation error in the range of [-Δϕ/2,Δϕ/2]. To reduce the quantization error, a precise estimator for spatial frequency is given in this section, based on the results of previous reconstruction.By the algorithm in Section4(1), the noise variance σ⌢g2, the range of spatial frequencies φ_1g,φ-1g,…,φ_Kg,φ-Kg, and the covariance matrix of signal waveform Ω⌢ can all be obtained, where ϕ_k(g) and ϕ-k(g) denote spatial frequency for the largest and the second largest peaks in the kth peak of the spectrum γ⌢g. Let Φ-k=ϕ1,…,ϕ_kg-Δϕ,ϕ-kg+Δϕ,…,ϕL and A-k=A(Φ-k). Ω⌢k can be obtained by crossing out two rows and two columns corresponding to ϕ_k(g) and ϕ-k(g). Σ-k=σ⌢g2I+A-kΩ-kA-kH. Equation (15) can be modified as follows in order to estimate spatial frequency precisely:(20)κηk,ϕk=ln⁡Σ-k+ηkaϕkaHϕk+trΣ-k+ηkaϕkaHϕk-1R^g,where ηk is the power of the kth signal.The target function for signal power and spatial frequency is(21)η⌢k,ϕ⌢k=argminη⌢k,ϕ⌢k⁡κηk,ϕk.It can be inferred from∂κ(ηk,ϕk)/∂ηk=0,(22)η⌢k=aHϕkΣ-k-1R⌢-Σ-kΣ-k-1aϕkaHϕkΣ-k-1aϕk2.Taking (22) into ∂κ(ηk,ϕk)/∂ϕk=0, (23) can be deduced, where d(ϕk)=∂a(ϕk)/∂ϕk:(23)ReaHϕkΣ-k-1aϕkaHϕkΣ-k-1R⌢g-R⌢gΣ-k-1aϕkaHϕkΣ-k-1dϕk=0.The precise estimation of spatial frequency is as follows:(24)ϕ⌢kg=argmaxϕk∈ϕ_kg,ϕ-kgReaHϕkΣ-k-1aϕkaHϕkΣ-k-1R⌢-R⌢Σ-k-1aϕkaHϕkΣ-k-1dϕk-1,k=1,…,K.In practical applications, the peak position can be obtained by searching in the smaller rangeϕ_k(g),ϕ-k(g), so that the spatial frequency can be precisely estimated. The computational complexity of the process above is moderate because it is conducted consecutively for K signals; the search range of which is only Δϕ.


### 3.3. DOA Estimation and Network Sorting

The spatial frequencies ofK signals can be obtained by algorithms listed in Sections 3.1 and 3.2.  K amplitudes s⌢1,s⌢2,…,s⌢K can be obtained from K peak positions from the estimated first-order moment M⌢S-g, where s⌢1,s⌢2,…,s⌢K are 1×N signal vectors. The number of samples is small, so phase-difference method is used to estimate frequency for K signals, with the following formula:(25)f⌢kg=12π∑n=0N-2Wnangle⁡s⌢k∗ns⌢kn+1,k=0,1,…,K-1,Wn=3N2N2-11-n-N/2-1N/22,n=0,1,…,N-1,where angle⁡(∗) means calculating the phase of a complex number.Therefore, the incidence angle can be estimated byϕk,n=fk,ndsin⁡(θk)/c with signal frequencies:(26)θ⌢kg=arcsin⁡cϕ⌢kgf⌢kgd.The incidence angle and frequency ofK 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 K different clusters with known number of sources as K. The frequencies corresponding to different spatial incidence angles are also categorized, so that different signals of different networks are sorted.


### 3.4. 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.ConsiderK frequency estimates f⌢kg, k=0,…,K-1; g=0,…,G-1, which is converted to binary image data after quantization by the following formulas:(27)f⌢kg,n=1,fmin+nΔf≤f⌢kg<fmin+n+1Δf0,else,n=0,…,BΔf-1,where Δf and B denote interval and bandwidth.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 elementb on the image f(x,y) are expressed as(28)fΘb=min⁡fx-x′,y-y′-bx′,y′∣x′,y′∈Db,f⊕b=max⁡fx-x′,y-y′+bx′,y′∣x′,y′∈Db,where Db is the domain of b.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 as11. By erosion and dilation, outliers can be excluded while normal values are unaffected, generating a clarified time-frequency matrix f⌢k∗(g,n). 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:(29)f⌢k∗g=f⌢kg,g∈G1f⌢kgm,gm=argmings∈G1gs-g,g∈G0,where(30)G0=g:∑n=0B/Δf-1f⌢k∗g,n=0,g=0,1,…,G-1,G1=g:∑n=0B/Δf-1f⌢k∗g,n=1,g=0,1,…,G-1.The modified time-frequency spectrum is turned into binary formf-k∗(g,n) 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.Figure 1
Four connected components and eight connected components.
(a)
Four connected components
(b)
Eight connected componentsThe sets of points of each hop areGk,0,…,Gk,vk,…,Gk,Vk-1, whose maximum and the minimum of time coordinates are g-vk and g_vk, indicating the moment of FH frequencies. The hop period estimate of the kth signal is as follows:(31)T⌢k=NVk-2fs∑vk=1Vk-2g-vk-gvk-1.Therefore, the frequency estimate ofvk hop of kth signal is(32)f⌢k,vk=1g-vk-g_vk+1∑g=g_vkg-vkf⌢k∗g.The DOA estimate ofkth signal is(33)θ⌢k=1G∑g=0G-1θ⌢kg.


## 3.1. Spatial Frequency Estimation

According to sparse Bayesian theory and (5), the probability density function of amplitudes from output data can be expressed as follows:(9)pXg∣S-g=πσg2IMexp-σg-2Xg-A-S-g22,where σg2IM denotes the variance of the gth frame.To facilitate notation, intermediate variableγg=[γ1(g),…,γL(g)]T is introduced to denote the power spectrum of incidence signal, and (10)s-n~N0,Ω,n=1,…,N,where Ω=diag⁡(γg) and s-(n) 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-g with respect to Xg can be expressed as(11)pS-g∣Xg;γg,σg2=pXg∣S-g;σg2pS-g∣γg∫pXg∣S-g;σg2pS-g∣γgdS-g=πΣS-g-Nexp-trS-g-MS-gHΣS-g-1S-g-MS-g,where MS-g and ΣS-g denote the first and second moments of S-g. Given γg and σg2, MS-g and ΣS-g can be estimated as(12)M⌢S-g=ΩA-Hσg2I+A-ΩA-H-1Xp,(13)Σ⌢S-g=Ω-ΩA-Hσg2I+A-ΩA-H-1A-Ω.From (12) and (13), it can be seen that the precision greatly depends on the prior distribution of γg and estimation of σg2. Precise results can be obtained only when γg and σg2 can correspond to the energy of incidence signal.To obtainγg more precisely, the observation data is used for optimization. The likelihood function of Xg with respect to γg is(14)pS-gXg;γg,σg2=∫pXgS-g;σg2pS-gγgdS-g=πΣS-g-Nexp-trXgHσg2I+A-ΩA-H-1Xg.Taking logarithm of (14) and omitting constant terms, the target function with respect to γg is(15)κγg,σg2=lnσg2I+A-ΩA-H+trσg2I+A-ΩA-H-1R^g,where R^g=(1/N)XgXgH is the estimated array output covariance matrix. By minimizing (15), γ^g, 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 γg. Every iteration of EM includes an E-step and an M-step. In E-step, the first-order moment M⌢S-g and second-order moment Σ⌢S-g are estimated by maximizing (11). In the M-step, the q+1 iteration of γg is estimated by calculating ∂κ(γg,σg2)/∂γg=0 given σ2:(16)γ⌢g,lq+1=M⌢S-gql•22N+Σ⌢S-gql,l,where M⌢S-gql∙ denotes the lth column of the matrix MS-g after the qth iteration and Σ⌢S-gql,l denotes the element on the lth row and lth column of the matrix ΣS- after the qth iteration, which can be deduced from (12) and (13).To improve the speed of iteration and avoid exception caused by many zeros inγg, this iteration can be modified as(17)γ⌢g,lq+1=M⌢S-gql•22N1-Σ⌢S-gql,l/γ⌢g,lq+ξ,where ξ is a small positive value.The estimated powerγ⌢g,l(q) can be obtained after qth iteration. The number of sources is known to be K. The first K 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:(18)ϕ⌢kqg=γ⌢g,k,1qγ⌢g,k,1q+γ⌢g,k,2qϕ_kqg+γ⌢g,k,2qγ⌢g,k,1q+γ⌢g,k,2qϕ-kqg,k=1,…,K,where ϕ_k(q)(g) and ϕ-k(q)(g) denote the spatial frequency of the largest and the second largest value of the kth peak in the spectrum γ⌢g(q), whose power is γ⌢g,k,1(q) and γ⌢g,k,2(q), respectively.Based on the rough estimation above, the spatial frequency estimationsΦ⌢(q)=ϕ⌢1qg,…,ϕ⌢Kqg in the qth iteration can be obtained. The unbiased estimation of corresponding noise is(19)σ⌢g2Φ⌢q=PAΦ⌢q⊥Xg22NM-K,where PA(Φ⌢(q))⊥=I-A(Φ⌢(q))A†(F⌢(q)), the notation (•)† denotes pseudo-inverse operation. When γg and σg2 meet the requirement for convergence, the iteration is terminated, and the power spectrum is estimated to be γ⌢g, which is taken into (12) to estimate the first-order moment M⌢S-g.


## 3.2. Precise Estimation of Spatial Frequency

The reconstruction in Section3.1 would lead to estimation error in the range of [-Δϕ/2,Δϕ/2]. To reduce the quantization error, a precise estimator for spatial frequency is given in this section, based on the results of previous reconstruction.By the algorithm in Section4(1), the noise variance σ⌢g2, the range of spatial frequencies φ_1g,φ-1g,…,φ_Kg,φ-Kg, and the covariance matrix of signal waveform Ω⌢ can all be obtained, where ϕ_k(g) and ϕ-k(g) denote spatial frequency for the largest and the second largest peaks in the kth peak of the spectrum γ⌢g. Let Φ-k=ϕ1,…,ϕ_kg-Δϕ,ϕ-kg+Δϕ,…,ϕL and A-k=A(Φ-k). Ω⌢k can be obtained by crossing out two rows and two columns corresponding to ϕ_k(g) and ϕ-k(g). Σ-k=σ⌢g2I+A-kΩ-kA-kH. Equation (15) can be modified as follows in order to estimate spatial frequency precisely:(20)κηk,ϕk=ln⁡Σ-k+ηkaϕkaHϕk+trΣ-k+ηkaϕkaHϕk-1R^g,where ηk is the power of the kth signal.The target function for signal power and spatial frequency is(21)η⌢k,ϕ⌢k=argminη⌢k,ϕ⌢k⁡κηk,ϕk.It can be inferred from∂κ(ηk,ϕk)/∂ηk=0,(22)η⌢k=aHϕkΣ-k-1R⌢-Σ-kΣ-k-1aϕkaHϕkΣ-k-1aϕk2.Taking (22) into ∂κ(ηk,ϕk)/∂ϕk=0, (23) can be deduced, where d(ϕk)=∂a(ϕk)/∂ϕk:(23)ReaHϕkΣ-k-1aϕkaHϕkΣ-k-1R⌢g-R⌢gΣ-k-1aϕkaHϕkΣ-k-1dϕk=0.The precise estimation of spatial frequency is as follows:(24)ϕ⌢kg=argmaxϕk∈ϕ_kg,ϕ-kgReaHϕkΣ-k-1aϕkaHϕkΣ-k-1R⌢-R⌢Σ-k-1aϕkaHϕkΣ-k-1dϕk-1,k=1,…,K.In practical applications, the peak position can be obtained by searching in the smaller rangeϕ_k(g),ϕ-k(g), so that the spatial frequency can be precisely estimated. The computational complexity of the process above is moderate because it is conducted consecutively for K signals; the search range of which is only Δϕ.


## 3.3. DOA Estimation and Network Sorting

The spatial frequencies ofK signals can be obtained by algorithms listed in Sections 3.1 and 3.2.  K amplitudes s⌢1,s⌢2,…,s⌢K can be obtained from K peak positions from the estimated first-order moment M⌢S-g, where s⌢1,s⌢2,…,s⌢K are 1×N signal vectors. The number of samples is small, so phase-difference method is used to estimate frequency for K signals, with the following formula:(25)f⌢kg=12π∑n=0N-2Wnangle⁡s⌢k∗ns⌢kn+1,k=0,1,…,K-1,Wn=3N2N2-11-n-N/2-1N/22,n=0,1,…,N-1,where angle⁡(∗) means calculating the phase of a complex number.Therefore, the incidence angle can be estimated byϕk,n=fk,ndsin⁡(θk)/c with signal frequencies:(26)θ⌢kg=arcsin⁡cϕ⌢kgf⌢kgd.The incidence angle and frequency ofK 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 K different clusters with known number of sources as K. The frequencies corresponding to different spatial incidence angles are also categorized, so that different signals of different networks are sorted.


## 3.4. 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.ConsiderK frequency estimates f⌢kg, k=0,…,K-1; g=0,…,G-1, which is converted to binary image data after quantization by the following formulas:(27)f⌢kg,n=1,fmin+nΔf≤f⌢kg<fmin+n+1Δf0,else,n=0,…,BΔf-1,where Δf and B denote interval and bandwidth.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 elementb on the image f(x,y) are expressed as(28)fΘb=min⁡fx-x′,y-y′-bx′,y′∣x′,y′∈Db,f⊕b=max⁡fx-x′,y-y′+bx′,y′∣x′,y′∈Db,where Db is the domain of b.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 as11. By erosion and dilation, outliers can be excluded while normal values are unaffected, generating a clarified time-frequency matrix f⌢k∗(g,n). 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:(29)f⌢k∗g=f⌢kg,g∈G1f⌢kgm,gm=argmings∈G1gs-g,g∈G0,where(30)G0=g:∑n=0B/Δf-1f⌢k∗g,n=0,g=0,1,…,G-1,G1=g:∑n=0B/Δf-1f⌢k∗g,n=1,g=0,1,…,G-1.The modified time-frequency spectrum is turned into binary formf-k∗(g,n) 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.Figure 1
Four connected components and eight connected components.
(a)
Four connected components
(b)
Eight connected componentsThe sets of points of each hop areGk,0,…,Gk,vk,…,Gk,Vk-1, whose maximum and the minimum of time coordinates are g-vk and g_vk, indicating the moment of FH frequencies. The hop period estimate of the kth signal is as follows:(31)T⌢k=NVk-2fs∑vk=1Vk-2g-vk-gvk-1.Therefore, the frequency estimate ofvk hop of kth signal is(32)f⌢k,vk=1g-vk-g_vk+1∑g=g_vkg-vkf⌢k∗g.The DOA estimate ofkth signal is(33)θ⌢k=1G∑g=0G-1θ⌢kg.


## 4. Algorithm Procedure and Computational Complexity

The procedure of joint parameter estimation and network sorting of multiple FH signals is summarized as follows.(1) Spatial Frequency Estimation Based on Sparse Bayesian LearningStep 1.
Initializeγ⌢g(0) and (σ⌢g2)0, where γ⌢g can be solved by the least-square solution of signal amplitudes. M⌢S-g(0)=A-H(A-A-H)-1Xg(σ⌢g2)0, γ⌢g,l(0)=M⌢S-g0l∙22/N, (σ⌢g2)0=0.1×XgF2/MN.Step 2.
Updateγ⌢g by (12), (13), and (17).Step 3.
Updateσ⌢g2 by (19).Step 4.
Repeat Steps2 and 3 until the convergence condition γ⌢g(q+1)-γ⌢g(q)2/γ⌢g(q)2<10-3 is met, and output γ⌢g and M⌢S-g.(2) Precise Estimation of Spatial FrequencyStep 5.
Conduct precise estimation of spatial frequency by (24), producing ϕ⌢1(g),…,ϕ⌢K(g).(3) DOA Estimation and Network SortingStep 6.
Estimate the instantaneous frequency and DOA by (25) and (26).Step 7.
Update frame for operation and repeat (1), (2), and Step6 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.(4) Parameter Estimation of FH Signals Based on Morphological FilteringStep 8.
The time-frequency spectrum generated by (27) is binarized and morphologically filtered and modified by (28).Step 9.
Detect the connected components of modified binary spectrum and extract the time-frequency point setGk,0,…,Gk,vk,…,Gk,Vk-1 of each hop. Estimate hop period, frequency, and DOA by (31), (32), and (33).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 Steps2 and 3. Therefore, it is in the order of O(ML2), while the computational complexity of the algorithm in [11] is in the order of O(M2L2). It could be easily seen that the computational complexity of the proposed method is lower than the algorithm in [11].


## 5. Simulation and Analysis

In this section, the proposed method is analyzed by simulation to verify factors contributing to its performance.Experiment 1.
Consider an ULA with eight elements and interelement spacing of 2.5 m. The bandwidth of receiver is[30MHz,60MHz]. 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.The simulated time-frequency spectrum and space-time pattern with SNR = 0 dB are depicted in Figures2 and 3. Parameter estimates are listed in Tables 1 and 2.Table 1
Estimate of DOA and period.
DOA 1/°
DOA 2/°
Period 1/us
Period 2/us
Real value
38
−44
100
50
Estimate
38.93
−43.01
101
50.6Table 2
Estimate of frequency sets.
Frequency set 1/MHz
Frequency set 2/MHz
Real value
42
35
53.5
48.7
47
53
34.3
42.5
49
56
39
Estimate
42.02
35.03
53.51
48.67
47.02
52.96
34.33
42.46
49.04
55.97
38.95Figure 2
Time-frequency spectrum and space-time spectrum after reconstruction.Figure 3
Time-frequency spectrum and space-time spectrum after morphological filtering.As seen in Figure2, 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 Tables1 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.Experiment 2.
Consider an ULA with eight elements and interelement spacing of 2.5 m. The bandwidth of receiver is[30MHz,60MHz]. FH signals are generated with incidence angles randomly selected in (-90°,90°) and frequencies randomly selected in [30MHz,60MHz]. The frequency interval between adjacent two points is greater than 1 MHz. Assume the two frequencies meet ϕ1,n-ϕ2,n>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 K 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.Figures4, 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:(34)ef=1I∑kVk∑k,ifk,vki-f⌢k,vki,eθ=1KI∑k,iθki-θ⌢ki,eT=1KI∑k,iTki-T⌢ki,where f(k,vk)(i) and f⌢(k,vk)(i) denote the frequency and estimate of the hop of the kth signal in the ith experiment; θk(i) and θ⌢k(i) denote the DOA and estimate of the kth signal in the ith experiment. Tk(i) and T⌢k(i) denote the period and estimate of kth signal in the ith experiment. I denotes the total number of experiments.Figure 4
Bias of frequency estimation in different SNRs.Figure 5
Bias of DOA estimation in different SNRs.Figure 6
Bias of period estimation in different SNRs.From Figures4, 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.


## 6. Conclusions

The reconnaissance and processing of multiple FH signals are challenging tasks in modern warfare. By sparse Bayesian reconstruction method, the parameter estimation and network sorting of multiple wideband FH signals are successfully achieved and precise estimation of parameters in low SNR conditions is conducted by the proposed method. The procedure of the algorithm has been deduced in theoretical analysis. And the excellent performance of the algorithm is verified by simulation.



---

*Source: 6129120-2017-09-06.xml*