Array model
The fractional Fourier transform of LFM signal
DOA estimation by FrFT-MSPICE
DOA estimation for wideband LFM signals with a few snapshots
Deliang Liu1,
Zhuo Li2, 3Email author,
Xiwei Guo1 and
Shen Zhao1
Received: 14 November 2016
Published: 7 February 2017
The wideband linear frequency modulation (LFM) signals are widely used in information systems. The conventional direction-of arrival (DOA) estimation algorithms of LFM signals rely on a large number of snapshots, some of which are not reliable in numerous practical applications such as underwater array processing. To solve the above problem, we present a modified sparse iterative covariance (MSPICE)-based estimation method in fractional Fourier transform (FrFT) domain to estimate the DOA of wideband LFM signals. First, we extend the original SPICE algorithm in FrFT domain with a specific transform order for wideband LFM signals. Then, we utilize the energy centrobaric modification method to make the original SPICE more accurate without adding more computational complexity. The simulation results demonstrate the effectiveness of the proposed method.
Direction of arrival
Fractional Fourier transform
Iterative approach
Wideband LFM
DOA estimation arises an important problem in a number of fields including sonar, radar, and wireless communication in recent years [1–4]. The classical delay-and-sum (DAS) method suffers from low resolution and high sidelobe levels, while the well-known multiple signal classification (MUSIC) method, which requires a large amount of snapshots, provides super-resolution DOA estimation for narrowband uncorrelated signals [5]. Furthermore, the wideband LFM signal, which is an important class of non-stationary signals, has been widely used in the abovementioned information systems. So many DOA estimation approaches have been proposed for them, such as MUSIC method in FrFT domain [6], Estimating signal parameters via rotational invariance techniques (ESPRIT) in FrFT domain [7] and ambiguity-function-based techniques [8], etc. However, none of these methods is able to provide high-angular resolution depending on very low snapshots which is the case when the environment being sensed by the array is stationary for a short duration of time.
IAA, SLIM, and SPICE algorithms [9] proposed by Stoica are sparsity-based techniques that can estimate the DOA depending on a few snapshots with high-resolution and low sidelobes. But they do not apply to wideband LFM signals, because the steering vector of the received signal model is time-variant. So some extensions to this method have been proposed to deal with the wideband LFM signals, such as spatial resampling (SR) method [10], the accuracy of which is not high due to using an approximate model for wideband LFM signals.
In this paper, a FrFT-MSPICE method is proposed to estimate the DOA for wideband LFM signals using a few snapshots based on the modified SPICE algorithm and FrFT. First, wideband LFM signal was represented as a sine wave with a single frequency in FrFT domain, so that the time-invariant steering vector can be obtained for the original SPICE algorithm. Then, we improve the original SPICE for higher accuracy and utilize it to estimate the DOA values in FrFT domain. The accuracy of the SPICE depends on the number of scanning points in the region, but increasing this number for higher accuracy will also dramatically increase the computational complexity. So, we utilize the energy centrobaric modification method [11] to improve the accuracy of SPICE without adding too much computational cost. Finally, the simulation results confirm the effectiveness of our algorithm.
2 Array model
Let us consider an active radar system as shown in Fig. 1. The linear array of the radar has M sensors uniformly placed along the x axis. The transmitter is located at the origin point. The distance between two adjacent sensors is d. T 1, T 2, …, T K are K Far-field targets at θ, where θ = [θ 1, θ 2,…, θ K]. Here K is usually unknown, so it is considered to be the amount of potential targets (scanning points) in the region, and it is much larger than the amount of actual ones. Only a few signal power estimates of the potential targets will be non-zero, so sparsity-based algorithm can be used in array processing applications.
The geometry of the array system
The transmitter emits a LFM signal that can be expressed as
$$ x(t)= a \exp \left( j2\pi {f}_0 t+ j\pi \mu {t}^2\right) $$
where a is the signal amplitude, f 0 is the centre frequency, and μ is the chirp rate.
The received signal at the sensor m can be expressed as the sum of K delayed versions of x(t), given by
$$ {x}^m(t)={\displaystyle \sum_{k=1}^K{\rho}_k{x}_k\left( t-{\tau}_k^m\right)}+{e}^m(t) $$
where e m (t) for m = 1, 2, …, M is the additive Gauss white noise at the sensor m, ρ k is the backscattering coefficient of target k for k = 1,2,…,K. \( {\tau}_k^m \) is the time delay of the kth signal traveling to the sensor m relative to the reference sensor (the first sensor), which can be expressed as
$$ {\tau}_k^m=\left( m-1\right) d \sin {\theta}_k/ c $$
where c is the wave speed.
The received signals at the sensors can be written in matrix form as
$$ \mathrm{x}(t)=\mathrm{a}\mathrm{s}(t)+\mathrm{e}(t) $$
where s(t) = [ s 1 (t), s 2 (t), … , s K (t)] T is the waveform vector, a = [a(θ 1), a(θ 2),…, a(θ K )] T is the steering vector with
$$ \mathrm{a}\left({\theta}_k\right)={\left[ \exp \left(- j2\pi {f}_0{\tau}_1+ j\pi \mu {\tau}_1^2\right) \exp \left(- j2\pi \mu {\tau}_1 t\right),\dots, \exp \left(- j2\pi {f}_0{\tau}_M+ j\pi \mu {\tau}_M^2\right) \exp \left(- j2\pi \mu {\tau}_M t\right)\right]}^T $$
If the bandwidth is small compared to the carrier frequency, a(θ k ) can be considered as time-invariant (i.e., the term exp(-j2πμτ m t) in (5) can be neglected), but for wideband LFM, it can not. So, the steering vector a(θ k ) of wideband LFM depends on the time t, and the SPICE algorithm for narrowband signals cannot be applied directly to wideband signals. To solve this problem, we derive FrFT method.
3 The fractional Fourier transform of LFM signal
The FrFT of signal x(t) is represented as [7, 12, 13]
$$ X\left(\alpha, u\right)={F}^p\left[ x(t)\right]={\displaystyle {\int}_{-\infty}^{+\infty } x(t){K}_{\alpha}}\left( t, u\right) d t $$
$$ {K}_{\alpha}\left( t, u\right)=\left\{\begin{array}{c}\hfill \sqrt{1- j \cot \alpha} \exp \left\{ j\pi \left[{t}^2 \cot \alpha -2 tu \csc \alpha +{u}^2 \cot \alpha \right]\right\},\alpha \ne n\pi \hfill \\ {}\hfill \delta \left( t- u\right),\kern15.25em \alpha =2 n\pi \hfill \\ {}\hfill \delta \left( t+ u\right),\kern13.5em \alpha =\left(2 n\pm 1\right)\pi \hfill \end{array}\right. $$
with p is the transform order, F p is the FrFT operator, K α (t,u) is the kernel function, α is the rotation angle, α = pπ/2.
As a generalization of the standard Fourier transform, the FrFT can be regarded as a counterclockwise rotation of the signal coordinates around the origin in the time-frequency plane, and the rotation angle is α. When α =2nπ + π/2, FrFT is equal to Fourier transform.
The FrFT of x(t) in (1) about angle α can be represented as
$$ X\left(\alpha, u\right)= a\sqrt{\frac{1+ j \tan \alpha}{1+\mu \tan \alpha}}\times \exp \left[ j\pi \frac{u^2\left(\mu - \tan \alpha \right)+2 u{f}_0 \sec \alpha -{f}_0^2 \tan \alpha}{1+\mu \tan \alpha}\right] $$
When α = α d = arctan(μs c 2), we can get
$$ X\left({\alpha}_d, u\right)= C \exp \left[ j2\pi {f}_0 \cos {\alpha}_d u\right] $$
where s c is the scale factor to normalize signals. In this paper, we use Ozaktas's fast sampling-type discrete FrFT method [14] to compute the digital values for FrFT, so, here, \( {s}_c=\sqrt{N}/{f}_s \) with N denotes the number of snapshots, f s is the sampling frequency. \( C= a \cos {\alpha}_d\sqrt{1+ j \tan {\alpha}_d} \exp \left[- j\pi {f}_{{}_0}^2 \sin {\alpha}_d \cos {\alpha}_d\right] \), which is a constant. Therefore, after the FrFT, the LFM signal becomes a sine wave with a single frequency f = f 0 cos α d .
Simulation is run to prove this dechirp property of FrFT with f 0 = 2.4 MHz, μ 0 = 8 × 1013Hz/s, f s = 600 MHz, N = 301. The result is shown in Fig. 2. The LFM signal in time domain (Fig. 2a) is changed into a sine signal (Fig. 2b) in FrFT domain with rotation angle α d .
The dechirp property of FrFT. a Time domain. b FrFT domain
According to the time shift characteristic of FrFT, we can get
$$ {F}^p\left[ x\left( t-\tau \right)\right]= \exp \left( j\pi {\tau}^2 \sin \alpha \cos \alpha - j2\pi u\tau \sin \alpha \right)\times X\left({\alpha}_d, u-\tau \cos \alpha \right) $$
At the sensor m, the received signal in FrFT domain which is reflected by target k can be expressed as
$$ \begin{array}{l}{X}_k^m\left({\alpha}_d, u\right)={F}^p\left[ x\left( t-{\tau}_k^m\right)\right]={\rho}_k aC \exp \left[ j2\pi \left({f}_0 \cos {\alpha}_d-{\tau}_k^m \sin {\alpha}_d\right) u\right]\\ {}\kern8.25em \times \exp \left\{ j\pi \left[{\left({\tau}_k^m\right)}^2 \sin {\alpha}_d \cos {\alpha}_d-2 f{ \cos}^2{\alpha}_d{\tau}_k^m\right]\right\}\end{array} $$
Because \( {\tau}_k^m \sin {\alpha}_d \) is very small, so we can get
$$ {X}_k^m\left({\alpha}_d, u\right)\approx {\rho}_k a C \exp \left( j2\pi {f}_0 \cos {\alpha}_d u\right)\times \exp \left\{ j\pi \left[{\left({\tau}_k^m\right)}^2 \sin {\alpha}_d \cos {\alpha}_d-2 f{ \cos}^2{\alpha}_d{\tau}_k^m\right]\right\} $$
Then (10) can be reformulated as
$$ {X}_k^m\left({\alpha}_d, u\right)={A}_k^m{S}_k(u) $$
$$ {A}_k^m= \exp \left\{ j\pi \left[{\left({\tau}_k^m\right)}^2 \sin {\alpha}_d \cos {\alpha}_d-2 f{ \cos}^2{\alpha}_d{\tau}_k^m\right]\right\} $$
$$ {S}_k(u)={\rho}_k a C \exp \left( j2\pi {f}_0 \cos {\alpha}_d u\right) $$
Put (3) into (13), we can get
$$ \begin{array}{l}{A}_k^m\left({\theta}_k\right)={A}_k^m= \exp \left\{ j\pi \right[{\left( m-1\right)}^2{d}^2{ \sin}^2{\theta}_k/{c}^2 \sin {\alpha}_d \cos {\alpha}_d\\ {}\kern5.75em -2 f{ \cos}^2{\alpha}_d\left( m-1\right) d \sin {\theta}_k/ c\left]\right\}\end{array} $$
Take FrFT on (2), we can get
$$ {X}^m\left({\alpha}_d, u\right)={\displaystyle \sum_{k=1}^K{X}_k^m\left({\alpha}_d, u\right)}+{E}^m\left({\alpha}_d, u\right) $$
Therefore, the FrFT of (4) can be given as
$$ \mathrm{X}=\mathsf{A}\mathsf{S}+\mathsf{E} $$
$$ S= diag\left\{{S}_1(u),{S}_2(u),\dots, {S}_K(u)\right\}, X={\left[{X}^1,{X}^2,\dots, {X}^M\right]}^T, A={\left[{A}_1,{A}_2,\dots, {A}_K\right]}^T $$
\( {X}^m={\left[{X}_1^m\left({\alpha}_d, u\right),{X}_2^m\left({\alpha}_d, u\right),\dots, {X}_K^m\left({\alpha}_d, u\right)\right]}^T,{A}_k={\left[{A}_k^1\left({\theta}_k\right),{A}_k^2\left({\theta}_k\right),\dots, {A}_k^M\left({\theta}_k\right)\right]}^T \)After FrFT, the steering vector will not depend on u, so we can use the SPICE algorithm to estimate the DOA values θ k .
4 DOA estimation by FrFT-MSPICE
The discrete form of (17) can be given as
$$ X(n)=\mathrm{A}\mathrm{S}(n)+ E(n) $$
where n = 1, 2,…, N is the number of snapshots.
Let P be a K × K diagonal matrix, whose diagonal contains the power at each angle on the scanning grid. The initial estimates \( {\widehat{P}}_k^{(0)} \) can be obtained using the SFLS method [15]
$$ {\widehat{P}}_k^{(0)}=\frac{{\displaystyle \sum_{n=1}^N{\left|{A_k}^H X(n)\right|}^2}}{{\left({A_k}^H{A}_k\right)}^2 N} $$
The noise covariance matrix R of X(n) can be given as
$$ {R}^{(i)}= E\left(\mathrm{X}{\mathrm{X}}^H\right)=\mathrm{A}{\widehat{\mathrm{P}}}^{\left( i-1\right)}{A}^H $$
Then in the i + 1 iteration, the power at each angle on the scanning grid can be updated as
$$ {p}_k^{\left( i+1\right)}=\frac{\left|{\beta}_k^{(i)}\right|}{w_k^{1/2}{\displaystyle \sum_{l=1}^K{w}_l^{1/2}}\left|{\beta}_l^{(i)}\right|} $$
$$ {w}_k=\frac{\left|\right|{a}_k\left|\right|{}^2}{tr(R)},\kern2em \left|{\beta}_k\right|{}^2={\left[\mathrm{QX}\left({X}^H X\right){X}^H{Q}^H\right]}_{k k}/{N}^2, Q=\mathrm{P}\mathrm{A}{\mathrm{R}}^{-1} $$
Then iterate (20)–(21) until \( {\widehat{P}}_k \) will not change obviously.
Finally, we search the positions of the spectral peak of \( {\widehat{P}}_k \) which are the final DOA estimates. The details of SPICE algorithm can be found in [9].
It is obvious that the accuracy of the DOA estimation depends very much on the angular scanning grid size. The higher accuracy we need, the smaller the grid size is, which means K is larger. But increasing K will dramatically increase the computational complexity of the algorithm. In order to make the DOA estimation more accurate without increasing too much complexity, we utilize the energy centrobaric modification method which is commonly used for modification of the frequency estimation in discrete spectrum.
Actually, because the existence of noise, when SPICE scans angles around the real DOA value, the energy estimation (\( {\widehat{P}}_k \)) will not be zero, which means there is "energy leakage" around the real DOA in the angular spectrum.
Suppose that there is a peak at k = k* with the energy P k* in the angular spectrum. Then, we search for P k * ', which is the largest P k around k*. Finally, we utilize these two spectral lines to estimate θ. But if P k * + 1 and P k * − 1 are nearly equal, we consider θ = k*, and it does not need modification.
$$ \widehat{\theta}=\left\{\begin{array}{c}\hfill \frac{k*{P}_{k*}+ k*\hbox{'}{P}_{k*\hbox{'}}}{P_{k*}+{P}_{k*\hbox{'}}}\kern2em \left|{P}_{k*+1}-{P}_{k*-1}\right|>\eta \hfill \\ {}\hfill \kern2.75em k*\kern4.5em \left|{P}_{k*+1}-{P}_{k*-1}\right|\le \eta \hfill \end{array}\right. $$
where P k * ' = max(P k * + 1, P k * − 1), η is a specified tolerance.
The whole iteration procedure can be summarized in Table 1.
Iteration procedue of MSPICE
5 Simulation results
We evaluate the performance of the proposed FrFT-MSPICE algorithm and compare it with other alternative methods.
5.1 Angular spectrum comparison
We investigate the performance of our algorithm about angular resolution, sidelobe level, and accuracy in the angular spectrum. We also compare it with other two algorithms: one is FrFT-MUSIC algorithm [6] which utilizes MUSIC algorithm combined with FrFT. But in order to apply for active radar, the rotation angle α has been changed from -arccotμ to arctanμ; the other is SR-IAA algorithm [10] which uses spatial resampling method combined with IAA.
The LFM signal with center frequency of 12 MHz and a bandwidth of 6 MHz is used in the simulation. The array contains 32 sensors which are uniformly spaced at half wavelength. The sampling frequency is 60 MHz and the number of snapshots is 2. We consider 3 targets at 70, 72, and 110°. The noise is assumed to be Gaussian random processes with zero mean and the SNR is 20 dB. The scanning grid for the three algorithms is uniform in the range from 1 to 180°, with 1° increment between adjacent grid points, so K = 180. The number of iterations is 10. One hundred independent runs are simulated to obtain the angular spectrum. Furthermore, each power value is normalized as 10 log 10(P k /P max), where P max is the maximum value of P k , k = 1,2,…,180.
We can see from Fig. 3a that the FrFT-MUSIC method can exactly estimate the DOA of the target at 110°, but it suffers from lower angular resolution than other two with low snapshots so that it can not separate the two targets at 70 and 72° and consider them as one target at 71°. From Fig. 3b, the DOA estimation of SR-IAA is 62, 65, and 117°. We can see that although the SR-IAA method has a higher resolution and lower sidelobe level than FrFT-MUSIC under low snapshots conditions and it can separate the two targets at 70 and 72°, but the DOA estimation errors are larger than other two algorithms because it uses an approximate model for wideband signals which omits the second order item of the original LFM model. Besides, we also found that this algorithm would be more accurate when the signal bandwidth is much smaller. From Fig. 3c, our FrFT-MSPICE algorithm produces the angular spectrum with lower sidelobe levels, higher angular resolution and higher accuracy in low snapshots conditions.
Comparison of the angular spectrum. a FrFT-MUSIC. b SR-IAA. c FrFT-CSPICE
5.2 Accuracy comparison between MSPICE and SPICE
In this section, we compare the accuracy of the proposed MSPICE with the original SPICE to demonstrate the effectiveness of using the modification method on SPICE. We use narrowband signal here because FrFT and SR may cause additional errors. Figure 4 shows the angle root mean squared error (RMSE) of the two algorithms for varying SNR from 0 to 30 dB. The scanning grid is uniform in the range from 1 to 180° with 1° increment. The target is placed at 60.0, 60.1,…, 60.9° respectively, and each RMSE in Fig. 4 is calculated as the sum of these ten individual RMSE. The results of 100 Monte Carlo trials are shown in the plot. We can observe that our modified SPICE algorithm has a higher accuracy than the original SPICE using the same scanning grid size without increasing much computational complexity.
RMSE comparison between CSPICE and SPICE
This paper has presented a FrFT-MSPICE method for the DOA estimation of wideband LFM signal. We extend the SPICE algorithm in FrFT domain so that the DOA of wideband LFM signals can be estimated with a few snapshots. The proposed method has high angular resolution and low sidelobe levels. We also utilize the energy centrobaric modification method in order to increase the accuracy of the SPICE algorithm without imposing too much additional computational burden. The simulation results have demonstrated the effectiveness of the proposed method.
This research was supported the by National Natural Science Foundation of China under Grant nos. 61601494, 61602346, 61501493.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Shijiazhuang Mechanical Engineering College, Shijiazhuang, China, 050003
Tianjin Key Laboratory of Wireless Mobile Communications and Power Transmission, Tianjin Normal University, Tianjin, China, 300387
College of Electronic and Communication Engineering, Tianjin Normal University, Tianjin, China, 300387
Z Li, K Liu, Y Zhao et al., MaPIT: An enhanced pending interest table for NDN with mapping bloom filter. IEEE Commun. Lett. 18(11), 1915–1918 (2011)View ArticleGoogle Scholar
Z Li, Y Chen, H Shi et al., NDFB-GSM-R: a novel high-speed railway communication system via named data networking. EURASIP J. Wirel. Commun. Netw. 48, 1–5 (2016)View ArticleGoogle Scholar
X Liu, Z Li, P Yang et al., Information-centric mobile ad hoc networks and content routing: a survey. Ad Hoc Netw. (2016). http://dx.doi.org/10.1016/j.adhoc.2016.04.005
Z Li, L Song, H Shi et al., Approaching the capacity of k-user MIMO interference channel with interference counteraction scheme. Ad Hoc Netw. (2016). http://dx.doi.org/10.1016/j.adhoc.2016.02.009
T Yardibi, J Li, P Stoica, M Xue, Source localization and sensing: a nonparametric iterative adaptive approach based on weighted least squares. IEEE Trans. Aerosp. Electron. Syst. 46(1), 425–443 (2010)View ArticleGoogle Scholar
F Ye, Y Li, S Yang, 2-D DOA estimation of LFM signal based on FRFT, in Proc. 1st ISSCAA (Harbin, China, 2006), pp. 1080–1083Google Scholar
P Luo, K Liu, W Shi, G Yan, 2-D DOA estimation of wideband LFM signals for arbitrary planar array. in Proc. 10th ICSP, 2010, pp. 307–310Google Scholar
N Ma, JT Goh, Ambiguity-function-based techniques to estimate DOA of broadband chirp signals. IEEE Trans. Signal Process. 54(5), 1826–1839 (2006)View ArticleGoogle Scholar
P Stoica, P Babu, J Li, New method of sparse parameter estimation in separable models and its use for spectral analysis of irregularly sampled data. IEEE Trans. Signal Process. 59(1), 35–46 (2011)MathSciNetView ArticleGoogle Scholar
Z Chen, J Li, P Stoica, KW Lo, Iterative adaptive approach for wide-band active sonar array processing. in Proc. IEEE OCEANS (Sydney, 2010), pp. 1–10Google Scholar
C Zheng, K Ding, Z Yang, Noise influence on frequency estimation accuracy from energy centrobaric correction method for discrete spectrum. in Proc. ICIA '09, 2009, pp. 1477–1481Google Scholar
R Tao, F Zhang, Y Wang, Research progress on discretization of fractional Fourier transform. Science in China Series F: Information Sciences 51(7), 859–880 (2008)MathSciNetView ArticleMATHGoogle Scholar
MF Erden, MA Kutay, HM Ozaktas, Repeated filtering in consecutive fractional Fourier domains and its application to signal restoration. IEEE Trans. Signal Process. 47(5), 1458–1462 (1999)View ArticleGoogle Scholar
HM Ozaktas, O Arikan, MA Kutay et al., Digital computation of the fractional Fourier transform. IEEE Trans. Signal Process. 44(9), 2141–2150 (1996)View ArticleGoogle Scholar
JD Scargle, Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data. Astrophys. J. 263, 835–853 (1982)View ArticleGoogle Scholar