An efficient joint estimation of wideband polynomial-phase signal parameters and direction-of-arrival in sensor array

We consider the joint estimation of the direction-of-arrival (DOA) and parameters of wideband polynomial-phase signals (PPSs) in sensor array. Unlike concurrent methods that require multidimensional searches, the proposed method requires 1D searches for all the parameters of interest. In this way, we can efficiently estimate the considered parameters in applications where large antenna arrays, containing tens or hundreds of sensors, are used. As special cases, we consider in detail the estimation of the second- and third-order PPSs. The former are estimated using the high-order ambiguity function (HAF), while the latter are estimated using the cubic phase function (CPF), known to outperform the HAF in terms of both accuracy and signal-to-noise ratio (SNR) threshold. In both cases, the estimation of the highest order parameter reaches the Cramér-Rao lower bound (CRLB), while the DOA estimation is above the CRLB for around 1 dB (second-order PPS) and around 6 dB (third-order PPS).


Introduction
An important application of polynomial-phase signal (PPS) estimation is related to the underwater monitoring of vessels and marine fauna [1][2][3], where large hydrophone arrays, containing tens or hundreds of sensors, are a common tool [4][5][6]. Numerous publications address the problem of estimating the PPS parameters along with the direction-of-arrival (DOA) [7][8][9][10][11][12]. Gershman et al. [9], consider the joint estimation of DOA and PPS parameters using a technique called the polynomial-phase beamformer. Although this approach is more efficient than the maximum likelihood (ML) technique, it still requires a search over a multidimensional parameter space. For example, the second-order PPS estimator, also referred to as the chirp beamformer, requires a search over a 3D parameter space. As a search tool, genetic algorithms (GAs) are used in [9]. In [10], monocomponent PPSs are considered and the high-order ambiguity function (HAF) is used to obtain the coarse estimates, which are in turn refined using the extended Kalman filter (EKF). The EKF can also be used with multicomponent PPSs [11]. An important recent advance has been proposed in [12], where multiple PPSs impinging on non-calibrated arrays, common for practical trials, are considered. The effects of non-calibration are avoided and signal components are separated using a blind separation technique [7,12,13]. In [12], the precise estimation is performed by means of the phase unwrapping for each component separately under the assumption that the considered signal is narrowband.
The signal estimation in underwater environment attracts a considerable attention in the recent past. In [2], a new approach for channel estimation in underwater acoustic communications, based on the Mellin transform, is introduced. Time-frequency analysis can be successfully used for underwater dispersive channel estimation [3]. An algorithm for azimuth/elevation direction-finding that enlarges the array aperture without introducing additional sensors and nonuniform interelement spacing, and that avoids the directioncosine ambiguity is proposed in [4].
In this article, we propose a method for the joint PPS and DOA estimation that efficiently reduces the search space to 1D. As such, the method is suitable for application in underwater acoustics where hydrophone arrays comprising numerous sensors are found. By calculating the high-order instantaneous moment (HIM) of the considered signal, we transform it to a form from which the highest PPS parameter and the DOA can be estimated very accurately and efficiently. The proposed method is described in detail for the cases of second-and thirdorder PPSs. We consider monocomponent PPSs impinging on a uniform linear array (ULA) of omnidirectional sensors, assuming that advanced techniques can be used for components separation and array calibration [7,12,13]. Furthermore, we do not introduce a restriction that the considered signal is narrowband as in [12] since such a restriction is hard to be satisfied in underwater systems [2]. The main tools used in the proposed technique are the HAF [14] and the cubic phase function (CPF) [15].
The rest of the article is organized as follows. The considered signal model, along with the polynomialphase beamformer, is described in Section 2. The proposed estimation method is described in Section 3, where the special cases of the second-and third-order PPSs are considered in detail. Simulation results and concluding remarks are given in Sections 4 and 5, respectively.

Array signal model
The model of a PPS impinging on a ULA with M omnidirectional sensors can be described as where a(θ, n) is the M × 1 array steering vector, x(n) the PPS, v(n) the M × 1 vector of complex Gaussian zero-mean noise, and N the number of samples. Without loss of generality, we will assume that N is odd. The Kth order PPS is defined as where A is the amplitude, a k is the kth order phase parameter and Δ is the sampling interval. The array steering vector can be modeled as [10] where the instantaneous frequency (IF) of the signal is given by and ω(n) is assumed to be constant during the time necessary for the wave to travel across the array aperture. The parameter ψ is related to the signal's DOA, θ, as follows: where d is the spacing between two adjacent ULA sensors and c is the propagation speed. Our aim is to estimate the vector of unknown parameters V = [θ, a 1 , ..., a K ] from y(n).
Gershman et al. [9], propose an estimator, known as the polynomial-phase beamformer, where the vector V is estimated by maximizing the following function: where (·) H represents the Hermitian operator. With (6), the dimensionality of the problem is reduced by one compared to the ML approach, but it still remains quite high. The applicability of the polynomial-phase beamformer is limited to lower PPS orders. To that end, in [9], a chirp is considered as the transmitted waveform, and the obtained estimator is referred to as the chirp beamformer. Such an estimator entails a 3D search and the GA is applied for the 3D search optimization. However, a proper GA setup requires a great number of performed tests and selections. Poorly chosen GA setup could lead to local optima that are far away from the true parameter values [16]. Furthermore, with the increase of the PPS order, the number of local optima in the optimization function increases, which in turn increases the probability that the GA lands on a local optimum.

DOA and PPS estimation algorithm
The Kth order PPS arriving on the mth ULA sensor can be rewritten as follows: In order to estimate the DOA and PPS parameters of x m (n), we propose to calculate the Kth higher-order instantaneous moment (HIM) of x m (n) as [14] with the first two HIM orders given by where τ is the time lag parameter. Note that the HIM is implemented through recursive auto-correlations. The Kth order HIM of x m (n) equals where m = 0, . . . , M − 1, |n| ≤ N−1 2 − (K − 1)τ ,and HIM K [x m (n), τ] represents the product of two complex sinusoids, one with index n and the other with m. The frequency of the former is ω n = 2 K-1 K!τ K-1 Δ K a K , whereas the frequency of the latter is ω m = 2 K-1 K!τ K-1 Δ K-1 a K ψ. The estimation of the highest order parameter a K and DOA θ therefore boils down to sinusoid frequency estimation. If we denote the frequency estimations of ω n and ω m asω n andω m , respectively, a K and θ can be estimated aŝ The frequency estimationsω n andω m can be obtained using the periodogram maximization procedure [17]. The discrete Fourier transform (DFT) of the HIM is referred to as the HAF. The estimation of a K thus requires the calculation of the DFT of HIM K [x m (n), τ] with respect to index n and the DFT maximization. Similarly, the θ estimation requires the calculation of the DFT of HIM K [x m (n), τ] with respect to index m and the DFT maximization. The DFT maximization requires 1D search [17].
The a K and θ estimates can be improved by averaging results over m and n, respectively. Since for lower signal-to-noise ratios (SNRs) the estimation can be plagued by outliers, we propose to perform an a-trimmed averaging instead of standard averaging. The a-trimmed averaging does not take into consideration a percentage of extreme estimates, which most probably correspond to outliers. The a-trimmed average of an N-element array X is defined as [18] Trim where Trim a [·] is the a-trimmed average operator, a the percentage of discarded elements, X s represents the array X sorted in ascending/descending order, and ⌊·⌋ and ⌈·⌉ represent the round down and round up operators, respectively.
Lower order PPS parameters can be obtained from the dechirped signalŝ by repeating the procedure defined by (10) and (11).
In the sequel, we will explain in detail the estimation of DOA and parameters of second-and third-order PPSs.

Estimation algorithm for K = 2
Let us consider the second-order PPS, i.e., case K = 2. Now we have The second-order HIM equals and parameters a 2 and θ can be estimated aŝ In (15) and (16), DFT n [·] and DFT m [·] represent the DFT operators with respect to n and m, respectively.
The parameter a 1 can be obtained by maximizing the DFT of with respect to index n and averaging the estimates obtained for all sensor indices m = 0,..., M -1.
In the considered K = 2 case, instead of performing a 3D search as proposed in [9], we are able to estimate all the parameters by performing three 1D searches per sensor. The estimates are improved by averaging results obtained for all sensors, or by using the a-trimmed average (13).
The estimation algorithm is given below, followed by the calculation complexity analysis. for Estimate a m 2 from the DFT of HIM 2 [x m (n), τ] calculated with respect to n, i.e., end for Estimate a 2 using the a-trimmed average operator as follows: end for Estimate a 1 aŝ In the algorithm's calculation complexity analysis, we will assume that the N-samples DFT calculation requires N log 2 N complex additions and multiplications, and that the sorting of an N-samples long real valued sequence requires N log 2 N comparison/exchange operations [20,Section 5.2.2]. In addition, 1D searches and scaling operations will not be accounted for in the analysis due to their low complexity.  N N θ N a 1 N a 2 ) operations, where N θ , N a 1 , and N a 2 represent the number of points in the θ, a 1 , and a 2 grids, respectively, used in the maximization procedure. Clearly, the proposed approach offers a significant computational cost reduction with respect to the chirp beamformer.

Estimation algorithm for K = 3
When the PPS order is K = 3, we have where A 0 = a 0 + a 1 mψ, For K = 3, the HAF is not the optimal solution for the PPS estimation since the A 3 estimation requires the calculation of HIM 3 [x m (n), τ], which incorporates two auto-correlations. Each auto-correlation increases the SNR threshold a by about 6 dB [21]. Therefore, we will use the CPF which offers lower SNR threshold and more precise estimation [15]. The CPF is defined as In noise-free case, the CPF is maximized at = φ m (n) = 2A 2 + 6A 3 (n ) = 2a 2 + 6a 3 (n ) + 6a 3 ψm, (19) where φ m (n) represents the second-order phase derivative of x m (n). Parameters A 2 and A 3 are estimated by locating maxima of the CPF calculated at two time instants and solving a set of two linear equations [15]. Therefore, in order to estimate A 3 , we need to perform one auto-correlation less compared to the HAF. After estimating A 3 , the DOA θ is estimated using (12), which entails the calculation of HIM 3 [x m (n), τ] according to (10). The parameter a 2 can then be estimated from A 2 (see (17)). The estimates of a 3 , a 2 , and θ can be improved by averaging over all sensors (a 3 and a 2 ) and time instants (θ).
Parameters A 0 and A 1 are estimated from the dechirped signalŝ The estimates of A 0 , A 1 , A 2 , and A 3 can be refined using the method proposed in [22] and in turn used to refine a 3 , θ, a 2 , a 1 , and a 0 , respectively. The refinement method is outlined in Appendix 1.
Again, all the parameters are estimated via 1D searches, as opposed to the ML approach and polynomial-phase beamformer that require 5D and 4D searches, respectively.
The estimation algorithm follows, along with the calculation complexity summary.
for m = 0 to M -1 Estimate  N a 2 , and N a 3 represent the number of points in the θ, a 1 , a 2 , and a 3 grids, respectively, used in the maximization procedure.

Estimation for higher and unknown orders
Underwater acoustic signals can be modeled by PPSs of order higher than three [23]. The proposed algorithm for joint estimation of the PPS parameters and DOA, presented in Section 3, works with an arbitrary PPS order. The highest order parameter a K and DOA are estimated using (11) and (12), respectively, while lower order PPS parameters are estimated from the dechirped signal (14). Keep in mind, however, that the SNR threshold in the HAF-based approach increases with the PPS order [21]. Then, if the underlying application requires estimation at lower SNR values, we could use the product HAF (PHAF) [19] or the hybrid CPF-HAF approach [24], instead of the HAF. These approaches are characterized by lower SNR threshold for higher PPS orders.
When the PPS order is unknown, we could use the strategy of increasing the HIM order until the DC component is obtained [19]. When the HIM order and that of the PPS coincide, a complex sinusoid with frequency proportional to the highest order parameter is obtained [14,19]. If, however, the HIM order exceeds the PPS order, a DC component is obtained. Another approach for determining unknown PPS order is presented in [25].

Simulation results
In our examples, we evaluate the proposed estimation method on the Kth order PPS x (n) = A exp (j(a 0 + a 1 (n ) + · · · + a K (n ) K )), where (nΔ) [-1, 1], for K = 2 and 3. The signal's DOA is θ = π/6, d = 1.5 m and c = 1500 m/s. The method's performance is evaluated by means of the root meansquare error (RMSE),calculated as where a is the true parameter's value,â l is its estimate in the lth Monte-Carlo simulation, and N MC is the number of Monte-Carlo simulations. Herein, N MC = 500. The SNR is defined as SNR = 10log 10 (A 2 /s 2 ). In this section, in addition to DOA, we will present results only for the estimation of parameter a 2 for K = 2, and parameters a 3 and a 2 for K = 3, since the estimation accuracy of a 1 is highly influenced by their accuracy. The results for a 1 will be reported briefly.
Example 1. In the first example, we consider the second-order PPS x(n) = A exp(j(19π + 5π(nΔ) + 11π(nΔ) 2 )) in three different scenarios. In all scenarios, the RMSE is calculated for a 2 and θ. In the first scenario, we calculate the RMSE versus the SNR that is varied from -6 to 14 dB in steps of 1 dB. The signal length is fixed to N = 257, i.e., Δ = 1/128, and the number of sensors is M = 100. The RMSE curves are depicted in the left two subplots (top and bottom) in Figure 1. Along with the RMSE curves, the CLRB curves are given (for the derivation of the CRLB see Appendix 2). Above the SNR threshold, which is around 2 dB, the a 2 estimation reaches the CRLB, whereas the θ RMSE is larger than the CRLB by about 1 dB.
In the second scenario, the RMSE is calculated for fixed N = 257 and SNR = 10 dB, and variable number of sensors M, that takes values from 10 to 190 in increments of 10. The RMSE curves are depicted in the middle two subplots in Figure 1. In the third scenario, for fixed SNR = 10 dB and M = 100, we varied the number of samples N from 101 to 1001 with a step of 100. The RMSE curves are given in the right two subplots in Figure 1. In both the second and third scenarios, the a 2 estimation reaches the CRLB, while the θ RMSEs are above the corresponding CRLBs by about 1 dB.
As for the estimation accuracy of a 1 , we obtained the difference between the a 1 RMSE and corresponding CRLBs of around 1 dB in all the considered scenarios.
The proposed approach is characterized by a very accurate estimation of parameter a 2 , while the RMSE of the DOA estimation is larger than the CRLB for a couple of dBs in all the considered scenarios. Interestingly, in the chirp beamformer, the discrepancy between the RMSEs of the estimated parameters and corresponding CRLBs increase as the SNR increases [9,Figures 4 and 5]. In addition, in the DOA estimation, the chirp beamformer's RMSE is closer to the corresponding CRLB with smaller number of sensors [9, Figure 7]. In our approach, however, the RMSE curves follow the corresponding CRLB curves with the increase of SNR and the number of sensors.
Example 2. Here, we consider the third-order PPS x (n) = A exp(j(19π + 5π(nΔ) + 11π(nΔ) 2 + 7π(nΔ) 3 )). The same scenarios and setups as in the previous example are considered and the obtained RMSEs ofâ 2 ,â 3 andθ are shown in Figure 2. The parameters are estimated using the algorithm described in Section 3.2. In the refinement method, the moving average (MA) filter length is M = 5, as suggested in [22]. The a-trimmed averaging is performed with a = 10%.
Again, the left three subplots in Figure 2 depict the RMSE versus SNR curves. Now, the performance threshold is around 0 dB. The a 3 estimation reaches the CRLB, whereas the a 2 and θ RMSEs are about 3 and 6 dB above the CRLBs, respectively.
In the case of a varying M (middle three subplots in Figure 2) and a varying N (right three subplots in Figure  2), the RMSEs ofâ 2 andθ are about 3 dB and 6 dB above the CRLBs, respectively.
The estimation accuracy of a 2 is noticeably lower than that of a 3 . This is due to the fact that the a 2 estimation is influenced by the estimation accuracy of three parameters, namely A 2 , a 3 , and θ, whose estimation RMSEs contribute to the estimation RMSE of a 2 .
The a 1 RMSE is around 2.5 dB above the corresponding CRLBs in all the considered scenarios.
Note also that the DOA estimation is worse than the PPS parameters estimation, which is also the case in the previous example. From (12), the DOA estimation accuracy is influenced by the estimation accuracy of a K and ω m , and non-linearity of the arcsine function. This accuracy loss is inherent to the proposed method, however it is justified by significant computational benefits.
In our approach, the DOA and PPS parameters are estimated from multiple sinusoids with frequencies proportional to the considered parameter. The obtained estimates are averaged using the α trimmed operator. Nevertheless, the proposed approach is by no means the general one. Multiple sinusoids can be combined in several ways to obtain an estimation more accurate than with a single sinusoid. For example, in order to increase the SNR, we could estimate the frequency from the sum of the spectra of the considered sinusoids. In order to improve the estimation accuracy, interpolation in frequency should be used. Another approach would be to multiply the spectra of considered sinusoids, an approach similar to the PHAF. Our initial simulations show that both of these approaches result in a lower SNR threshold than with the proposed approach, but with worse accuracy.

Conclusion
The calculation complexity reduction is a very important goal in applications where large antenna arrays are  used. The underwater monitoring of vessels and marine fauna is one such application. In this article, we proposed an efficient method for the estimation of DOA and parameters of PPS impinging on a ULA. As opposed to concurrent methods, it significantly reduces the computational complexity without a significant loss in accuracy. All the parameters are estimated through 1D searches. In addition, the proposed method can be used for parameter estimation of an arbitrary order PPS, which is very important in underwater environment where acoustic signals can be modeled by PPSs of order higher than three. The estimation of the second-and third-order PPSs are considered in detail. In both cases, the estimation of the highest order parameter achieves the CRLB, while the DOA estimation is above the CRLB for around 1 dB (second-order PPS) and around 6 dB (third-order PPS). Future research will consider estimation accuracy improvements, especially of the DOA, as well as providing more robust results in terms of the SNR threshold decreasing.

Appendix 1: PPS estimation refinement
Here, we revisit the PPS estimation refinement method proposed in [22]. The considered signal is where x(n) is the Kth order PPS and v(n) is zero-mean complex Gaussian noise. We will assume that we already have the coarse estimates of all the PPS parametersâ k , k = 1, 2, ..., K. The coarse estimate of a 0 is not needed at this point. It will be estimated along with the refinements of other parameters. The refinement algorithm is given below.