Performance Analysis of Fractional Fourier Domain Beam-Forming Methods for Sensor Arrays

ABSTRACT The problem of estimating the direction of arrival (DoA) of multiple far-field moving targets, producing broad-banded chirp signals, in both active and passive mode of operation is addressed. As the chirp signal, commonly used in both radar and sonar systems is better processed in fractional Fourier domain, the detection is done using fractional Fourier transform (FrFT) in the proposed work. Subsequently, the DoA estimation is achieved using both conventional and subspace methods. Even though FrFT beam-forming method has been studied in the past, its performance in respect of varying signal-to-noise ratios (SNRs) and multiple target detection has not been extensively analyzed yet. The results of simulation experiments reported in the present work show that FrFT enjoys better estimation accuracy over the conventional frequency domain beam-forming for chirp signals using Fourier transform (FT), in terms of computational efficiency, accuracy and resolution with low SNR, limited snapshots and sensors, and spatially coherent/multi-path signals. The performance metrics used for the study are (i) root mean square error, (ii) 3-dB beam width and (iii) CPU time for accuracy, resolution and computation time, respectively. It is seen from computer simulations that MUSIC outperforms other DoA estimation techniques from the performance curves for both active and passive systems. Graphical Abstract


Introduction
Sensor array signal processing has been an active area of research which deals with the problem of direction of arrival (DoA) estimation and beam-forming in a wide variety of applications such as radars [1,2], sonars [3], seismology [4], as well as in communication [5]. Supported by a rich collection of literature, the conventional approach to DoA estimation is the Delay Sum (DS) or Conventional Beam-Former (CBF) [6], which has been found to be very effective for narrow-band signals; but suffers from a poor detection accuracy and is not very effective in resolving multiple sources due to the presence of high side lobe levels. In order to improve the resolution, Capon et al. [7] developed the Minimum Variance Distortion-less Response (MVDR), by normalizing the response in all directions; but keeping the response constant in the direction of interest. But this algorithm is computationally complex due to the requirement of inverting the full-rank spatial covariance matrix, dimensioned by the size of the aperture. Subspace-based methods such as spectral Multiple Signal Classification (MUSIC) [8], spatial smoothed MUSIC [9] Root-MUSIC [10], Minimum Norm (MN) [11] and the Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) [12] also result in the high resolution of multiple targets, but they are largely applicable to narrow-band sources. CONTACT Sreekumar G. ggsreekumar80@gmail.com However in the case of broad-band signals, the theories are not that established as compared to narrow-band. Traditional methods of DoA estimation cannot be used as such because the steering vector of the broad-band signal is changing with time. The usual approach for incoherent broad-band signal is to separate the data into several narrow-band frequency bins, perform DoA estimation for each frequency bin and then construct the final estimate by statistical average [13,14]. For coherent source, each frequency is focused to a reference focusing frequency and then a single correlation matrix is obtained by adding these focused correlation matrices, to which the narrow-band DoA estimation methods could be applied [15][16][17]. The performance of incoherent method suffers severely at low SNR, while the performance of coherent method is likely to depend on the selection of the focusing frequency [18]. The DS beam-forming implemented in time domain fails to correct for all the frequencies in a broad-band. On the other hand, the Fourier based beam-forming necessitates the steering corresponding to each frequency in the band, where the frequencies are assumed to be stable during the time window selected. However, in the case of chirping signals, the instantaneous nature of the frequency does not permit the steering of the beam for each frequency in the chirping band, for a given observation window. In this context, for the chirping signals, the beam-forming using fractional Fourier transform (FrFT), where the FrFT parameters are chosen based on the chirp signal, can improve the performance of DoA estimation in terms of computational efficiency, accuracy and resolution. Also by using FrFT, the detection at low SNR gets improved, since the FrFT with the matching parameters, converts the broad-band signal into its narrow-band/sparse representation in the fractional domain [19,20]. Even though FrFT beam-forming method has been studied in the past [21,22], its performance in respect to spatial resolution at low SNRs has not been analyzed with clear cut conviction.
Chirps are signals that exhibit a change in instantaneous frequency with time and are of particular interest in active and passive systems of radars and sonars [23]. Active sonar system using chirp signals are designed for detecting and real-time imaging of objects buried under seabed [24]. The applications also span other areas that include acoustic communications, seismic surveying and ultrasonic imaging. Chirp signals occur in nature in the form of birds and bat sounds, bottleneck dolphin and minke whale whistles. For moving source problems in passive sensor arrays, the chirp signals are produced due to Doppler effect when the sinusoidal source is accelerating [21]. Some active systems of radars and sonars use chirp signal for transmission, as they provide good detection sensitivity of low Doppler targets with high range resolution [23,25].
Active sensor array uses the chirp signals for transmission as different waveforms like linear frequency modulated (LFM), quadratic frequency modulated (QFM) and hyperbolic frequency modulated (EFM) signals [26], which when reflected from a target provide the receiver with very good immunity to noise and reverberation, thereby improving the basis for detection and estimation of range, direction and velocity [25,27,28]. Since the receiver system has the advantage of the knowledge of the transmitted waveform which is the chirp, an FrFT processor becomes the apt choice as the receiver. This is because chirps represented at its optimum FrFT order would produce impulses. However the use of FrFT is restricted to linear chirps as it only rotates the axes of the time-frequency plane. Sahay et al. [29] introduced the idea of modified FrFT for the analysis of higher order chirps by the proper selection of the kernel. The resulting narrow-banded form of the chirp signal then produces detection output with low side lobes, with accurate DoA measurement. In forming the beam, the challenge lies in the formulation of the array steering vector for the fractional Fourier domain and also the transformation equation from natural frequency to fractional frequency.
Passive systems intercept active chirp transmissions from enemy targets and estimate the parameters of the signal and hence are also called intercept systems. Unlike the active case where the received signal is known, a one-dimensional search is necessary to find the optimum order with which the peak position in the corresponding FrFT output is obtained. The peak position in the FrFT output depends on three parameters which are the chirp start position, start frequency and the bandwidth [30]. Traditional uniform arrays such as uniform linear array (ULA), uniform planar array (UPA) and uniform circular array (UCA) can only detect N-1 number of sources using N number of physical sensors. Extension of FrFT and modified FrFT to non-uniform arrays such as nested and coprime arrays are also recently studied which have the attractive capability of providing enhanced degrees of freedom [26].
In this paper the problem of estimating the DOAs of multiple broad-band far-field moving chirp targets using both active and passive sensor arrays is addressed. In active systems the echoes received from various directions of a transmitted chirp pulse is collected at the array and FrFT beam-forming is done to estimate the directions. In passive scenario distinct chirps from different directions are collected to estimate the multiple DoAs. The traditional and subspace beam-forming methods are used to study the effect of chirp on active and passive sensor arrays. The novel part of the work is the derivation of the array steering vector for FrFT domain and the transformation equation from natural frequency to fractional frequency. Also an improvement in performance of the FrFT method is observed over the fast Fourier transform (FFT) method in terms of computational efficiency, accuracy and resolution with low signal-to-noise ratio (SNR), limited snapshots and sensors, and spatially coherent/multi-path signals. The performance matrices used for the study are root mean square error (RMSE), 3-dB beam-width and CPU time for accuracy, resolution and computation time, respectively.
The rest of the paper is organized as follows. In Section 2, the chirp signal model for ULA is presented, in terms of the frequency components of the chirp. Then the relationship between linear chirp and FrFT is derived in Section 3. DoA estimation methods, forming beams based on FFT and FrFT are discussed in Section 4. The theory of multi-path effect and the modifications to the DoA estimation methods are proposed in Section 5. In Section 6, the numerical results of the detailed evaluation of the performance of a variety of DoA estimation methods based on both FFT and FrFT are presented and to show the relative variations in performance, for active and passive operations. The conclusions and directions for further work are summarized in Section 7 and Section 8, respectively.

Notations
Matrices are denoted by uppercase letters in boldface (e.g. A). Vectors are denoted by lowercase letters in boldface (e.g. a). Superscript H denotes Hermitian or transpose conjugate. xðtÞ denotes a continuous time signal and xðf Þ its FT. x p ðuÞ represents the FrFT of the signal xðtÞ with p denotes the fractional order of the transform. The symbol j ¼ ffiffiffiffiffiffi ffi À1 p denotes the imaginary part of a complex number.

Chirp Signal and Array Modeling
Consider a general uniform linear antenna array of N isotropic elements for a narrow-band model. Let a (θ; f ) be the N × 1 array manifold (steering) vector corresponding to the direction θ represented as follows: aðθ; f Þ ¼ ð1 exp½Àj2πf τðθÞ ::::: exp½Àj2πf ðN À 1ÞτðθÞÞ T (1) where f is the center frequency of the narrow-band signal and the inter element delay, τðθÞ ¼ dsinðθÞ=v. The parameters d and v represent, respectively, the distance between two sensors and the signal speed. d is taken as λ=2 to prevent aliasing in the spatial domain where λ is the wavelength.
Assume D narrow-band signals centered on frequency f impinging on the array from directions θ 1 ; θ 2 ; θ 3 ; :::θ D located at far-field regimes. The time series of data at the q th sensor is given by where τ q ðθÞ = (q − 1)τðθÞ is the time delay of the narrow-band signal at the q th sensor. The noise n q ðtÞ is assumed to be additive white Gaussian noise and uncorrelated from the sources. However the signal of interest in the present discussion is not narrowbanded; but the broad-banded linear chirp which is represented by the following equation: where b is the chirp rate, f 0 is the start frequency and c is the initial phase of the chirp. The term bt 2 þ f 0 t þ c defines the phase and its first derivative 2bt þ f 0 represents the instantaneous frequency. Therefore three parameters define a chirp completely viz. start frequency, chirp rate and the duration over which t is defined. In frequency domain representation the received data vector in (2) with the broad-banded signal as input has the form where Aðθ; f k Þ ¼ ½aðθ 1 ; f k Þ aðθ 2 ; f k Þ::::aðθ D ; f k Þ represents the array steering matrix for the k th bin and sðf k Þ ¼ ½s 1 ðf k Þ s 2 ðf k Þ:::s D ðf k Þ T represents the k th FFT coefficient of the signal vectors. nðf k Þ ¼ ½n 1 ðf k Þ n 2 ðf k Þ:::n N ðf k Þ T is the N × 1 noise Fourier coefficient vector. The term Aðθ; f k Þsðf k Þ represents inverse beam-forming vector which delays the plane wave to reach the spatially separated sensors.

Linear Chirp and FrFT
FrFT is a linear transformation mostly suited for chirps as it decomposes the signal into an orthonormal basis of chirps. It is a one parameter generalization of classical FT. Namias introduced the idea of FrFT in the area of quantum mechanics for the solution of differential equation problems [31]. The discrete implementation of FrFT was put forward by Ozaktas [32] and the code used in the simulations was derived from that. Since then, a number of applications have been developed, mostly in the field of optics [33]. Recently it is being applied to the area of sensor arrays [34,35]. The FrFT of a signal s(t) is defined as follows: where p is the order of FrFT and can be any real number between 0 and 1. When p = 1, FrFT reduces to the classical FT and the inverse FrFT is obtained by substituting p = −1. K p ðt; uÞ is the kernel of the FrFT defined as follows: The relationship between the chirp rate b and the optimum order p opt depends on the digital sampling scheme [36,37] used and is given by where f s is the sampling frequency, ϕ is the anticlockwise rotation angle of the transform and M is the processed data length. Conversely this can be used to compute the chirp rate given the true order. A onedimensional search is required to find the optimum order for passive systems whereas for active systems the complete information is available with the transmitter.
Using (5), the FrFT of the chirp signal represented in (3) with zero initial slope can be obtained as follows: where T is the pulse width of the chirp. As the chirp signal represented at its optimum order produces an impulse which is given by the peak point of (8).
Corollary: The theoretical peak point of (8) which transforms the natural frequency to fractional frequency is given by u ¼ ðf 0 þ bTÞsinðϕÞ ¼ f mid sinðϕÞ where f mid is the center frequency of the chirp.
Proof: Equation (8) can also be expressed as follows: The peak value of the expression is obtained when the term 0:5T cot ϕ þ bT þ f 0 À u cscϕ goes to zero which implies Neglecting 0:5T cot ϕ in comparison with f 0 þ bT for practical cases implies Equation (8) also represents the FrFT of the incoming chirp signal at the first sensor which is the reference sensor. For the q th sensor, the FrFT of the signal which is delayed by τ q ðθÞ can be obtained as per the property of FrFT as follows: F p ½sðt À τ q ðθÞÞ ¼ S p ðu À τðθÞ cos ϕÞÂ exp j2πðq À 1Þ τ 2 ðθÞ 2 sin ϕ cos ϕ À uτðθÞ sin ϕ h i As the time delay between two sensors τðθÞ is a small value and hence so its square, the above equation can be approximated as Therefore the array manifold vector in the FrFT domain can be written in terms of ψ ¼ 2πuτðθÞ sin ϕ as follows: a p ðθ; uÞ ¼ ð1 exp½Àjψ exp½Àj2ψ ::::: exp½ÀjðN À 1ÞψÞ T Observation 1: a p ðθ; uÞ ¼ aðθ; f Þ for u ¼ f mid sinϕ and p ¼ 1: Proof: The last term of a p ðθ; uÞ is exp½Àj ðN À 1Þψ = exp½Àj2πðN À 1ÞuτðθÞ sin ϕ. Substituting u ¼ f mid sinϕ and ϕ ¼ pπ=2, implies exp½Àj2πðN À 1Þ f mid sinðpπ=2Þ τðθÞ sinðpπ=2Þ. When p = 1, the expression reduces to a p ðθ; uÞ ¼ ½1 exp½Àj2πf mid τðθÞ ::::: exp½Àj2πðN À 1Þf mid τðθÞ T . This resembles aðθ; f Þ as in (1) where f mid ¼ f . Proof: As the chirp rate b is negative for a down-chirp, the rotation angle ϕ ¼ tan À1 ðf 2 s =2bMÞ becomes negative which in turn makes u negative. Therefore the steering vector a p ðθ; uÞ ¼ ð1 exp½jψ exp½j2ψ ::::: exp ½jðN À 1ÞψÞ T is the complex conjugate of the same band up-chirp.

Methods of Direction of Arrival Estimation
The DoA estimation methods are classified into conventional or non-subspace methods which includes CBF [6], MVDR [7] and the subspace approaches are the eigen value decomposition methods which includes MUSIC [8,9] and MN [10] methods.
The spatial energy spectrum for CBF, MVDR, MUSIC and MN, respectively, can be re-casted in the following way for chirp signals represented in FrFT domain as follows: where e is a vector of all zeros except a 1 at the 1 st position and the MVDR weight mðθÞ is given by U n is the noise eigenvector obtained from the eigen decomposition of R x ðuÞ, the spatial covariance matrix. R x ðuÞ is obtained from the array snapshots as follows: where L is the number of data blocks taken having T duration each which is the chirp pulse width and x l ðuÞ is the FrFT of the l th chirp block. It is well known that in the conventional frequency domain, the energy given by (15)(16)(17)(18) is calculated for all frequencies in the transform domain and summing up the total energy to mark the DoA. However as the FrFT of the chirp signal at its optimum order produces an impulse at u opt irrespective of the bandwidth, it is worth enough to take the peak point xðu opt Þ alone to calculate R x ðuÞ. This makes the FrFT beam-forming a single sample processing method gaining better time complexity for chirp signals whereas the FFT may not yield any meaningful result due to the non-stationary nature of the chirp signal. The peak is clearly detected even at −15 dB for FrFT case which leads to accurate bearing estimation whereas in the case of FFT even at 0 dB the spectrum is not visible [22] which reveals the inefficiency of FFT method to use for the beam-forming of chirp signals at low SNRs. This is because at low SNR even a single outlier from one narrow-band component can potentially lead to inaccurate estimates through the averaging process thus affecting the performance severely. Also as the energy of FFT beamformer is obtained by repeating over each frequency bin, the computational time required is O(K) whereas it becomes O(K = 1) for FrFT where K is the total frequency bins.
As beam-forming is a block processing method, the number of blocks taken should be much more than the number of sensors to require full rank for the spatial covariance matrix. This enables to perform matrix inversion which is needed for MVDR.
Result: Minimum number of data blocks in the beam-forming process required to get a full rank spatial covariance matrix is given by the formula Proof: The sensor spacing d is chosen to be λ min =2 where λ min = v=f max is the wavelength corresponding to the highest frequency of the incoming broad-band signals. This spacing ensures that there is no aliasing in the spatial domain.
Substituting d in (21) gives As per Nyquist-Shannon sampling theorem, f s ! 2f max . Therefore, the minimum number of data blocks L min ¼ N þ 1 which implies the number of blocks should be greater than the number of sensors. This makes the column full rank which is D, provided the sources are uncorrelated yielding a positive definite matrix R x ðuÞ which can be inverted.

Multi-Path Signal and DoA Estimation
Multi-path signal is often considered as an interference to be removed. It is generated due to reflection and diffraction of signal by various obstacles. The modified array steering vector with R multi-path components is expressed as follows: where aðθ r ; f Þ and c r are the steering vector and attenuation coefficient for the r th multi-path component, respectively. The DoA estimation methods cannot be used as such due to the coherence between the direct path and the multi-path signal. This is because the multipath components will get merged into the signal source resulting in a single peak value at the FrFT output which may not map to the direction of the true source.
The above problem can be solved by first decoherencing the coherent signals before giving to the DoA estimation algorithms. Spatial smoothing (SS) is one such technique which is extensively used which includes forward spatial smoothing or forward-backward spatial smoothing. In the present work forwardbackward spatial smoothing [38] is used. After decoherencing the coherent signals, the true source direction corresponds to the one with the maximum peak in the spatial spectrum as the multi-path components will get more attenuated. As SS is done on FrFT domain, the number of calculations is proportional to the number of sub-arrays selected whereas for FFT method it is proportional to the number of sub-arrays multiplied by the number of frequency bins which once again highlight the computational power of FrFT.
The sensor array works equally well for radar and sonar applications, with only difference in the wave properties leading to the rearrangement of sensor positions as shown in Table 1. In this case the sonar array is considered with the signal propagation speed of 1500 m/s.

Single Source
To demonstrate the power of FrFT method over FFT a single chirp source having f start and f end as 3000 Hz and 5000 Hz, respectively, arriving at 10 o with −15dB SNR is used for the comparison. Figures 1 and Figure 2 show the DoA estimation using CBF, MVDR, MUSIC and MN beam-formers for FFT and FrFT methods, respectively. The sampling frequency taken is 12,800 Hz with a total of 25,610 time samples and the pulse width of chirp as 200 ms. The samples of the impinging signal are divided into L = 10 blocks with each block having M = 2561 snapshots. In each block, the 2561 samples are converted into frequency domain or fractional frequency domain which are then processed using four different beam-forming algorithms. The array consists of N = 8 elements and the signal is corrupted with an AWGN noise of SNR = −15 dB.  FrFT beam-forming method is able to clearly locate the source with high accuracy and resolution whereas the FFT method fails for the broad-band chirp source. The performance analysis of the FrFT beam-former using CBF, MVDR, MUSIC and MN DoA estimation methods is evaluated using three parameters that are the resolution, accuracy and computational power. The resolution comparison using 3 dB beam-width is shown in Table 2. It is seen from the table that MUSIC algorithm performs the best followed by MVDR, MN and CBF has poor resolution. This is due to the strong orthogonality principle employed in the formulation of MUSIC algorithm.
As the resolution is clearly visible from the energy spectrum plot, accuracy is measured using RMSE plot. The variation of RMSE as a function of SNR for the various FrFT beam-formers with single source angle at 10 o is shown in Figure 3, averaged over 100 Monte Carlo simulations. The error is seen less for MUSIC and CBF method and more for MVDR and MN. The third and most relevant metric is the computational power which is represented in Table 4. The average CPU time for all FrFT beam-formers is found to be less than one second whereas FFT beam-formers producing error estimate need more than a minute to finish the same task. The peak side lobe level is also shown in Table 3 to evaluate the performance of various FrFT beam-formers. Lower the side-lobe level, better the beam-former performance that is observed best for MUSIC and worse for CBF and MN. Taking into account the various performance metrics used, MUSIC outperforms other FrFT beamforming techniques.

Multi-Path Sources
A single chirp source having f start and f end as 3000 Hz and 4000 Hz, respectively, arriving from À 20 o with −15dB SNR is reflected from two paths with incident angles À 40 o and 10 o . The attenuation coefficient for the multipath components are 0.4+ j0.8 and −0.3-j0.7. Figure 4 shows the DoA estimation using various SS-FrFT beamformers in the multi-path environment. The SS applied to subspace methods (SS-MUSIC and SS-MN) increases the rank of the signal subspace and is able to clearly resolve both the direct and multi-path components. Subspace methods detects all impinging directions corresponding to the signal eigenvalues. As the multi-path components are more attenuated compared to the direct, the one with maximum energy corresponds to the true direction. But in the case of conventional methods such as SS-CBF and SS-MVDR, the signal with maximum energy content will emerge as the source whereas the multi-path components with feeble energy is treated as equivalent as noise which gets attenuated in the overall process.

Passive Systems
In the case of passive systems, as no information about the received chirp signal is available, a one-dimensional binary search is done to estimate the value of optimum order. On the given block of data, FrFT is done for different values of p from 0 to 1 and the maximum peak value is selected. The algorithm implemented below can generate up to an accuracy of four decimals [29].    Algorithm: Iterative computation of optimum order for passive arrays using FrFT.
Input: Chirp signal with f start ¼ f 1 and f end ¼ f 2 .
Initializations: Select two orders p 1 = 0.25 and p 2 = 0.75 1: Compute FrFT for the selected orders. 2: Out of p 1 and p 2 select the optimum FrFT order as p opt corresponding to the highest peak.
For the simulation of multiple passive sources, four chirp signals (Chirp 1-Chirp 4) of different bandwidths/ chirp rates and center frequencies are used as given in Table 5. A down chirp (Chirp 5) is also included in the study with the same band as that of Chirp 1. The sources are located at far field impinging on an 8 element ULA from DoAs ( À 30 o ; À20 o ; À10 o ; 0 o ; 10 o ) with all equal power. Due to the change in the bandwidth, p opt value is calculated from (7) for each chirp. Then parallel filtering for each optimum order is performed and summed to obtain the effective FrFT output. The FrFT scheme involving MVDR, MUSIC and MN beam-formers is able to clearly resolve five sources with different chirp rates which are closely located at very low SNR as illustrated in Figure 5 with poor resolution only for CBF. The high spatial energy is seen for MUSIC and MVDR and is reduced when it comes to MN.

Active Systems
In the case of active sensor arrays, a single transmitted chirp with f start = 3000 Hz and f end = 4000 Hz is reflected from multiple moving targets located at different positions. The chirp signal undergoes a shift in the frequency band due to Doppler effect [21] and the reflected chirps are shown in Table 6. The frequency shifts are found different by using the Doppler equation for the three moving targets with velocities 5, 10 and 15 m/s. Same FrFT order was retained for the reflected chirps as the shift in the bandwidth was observed only within 6 to 20 Hz. But due to the shift in the frequency band, the FrFT bin value gets shifted for the received chirps as shown in Figure 6. Three distinct peaks shown in Figure 6.b corresponds to the three active moving chirp targets with overlapping frequency bins. This results in the resolution of the closely spaced sources for the FrFT beam-forming methods as illustrated in Figure 7, whereas the conventional FFT method fails.

Conclusion
In this paper, the efficacy and performance of applying FrFT beam-former for broad-banded chirp signals in active and passive mode using linear array is compared for a variety of processing techniques, restricting the   beam-former to the point where the FrFT peaks. It is seen that accuracy, resolution and computation power improves remarkably over the conventional method in the light of FrFT. The performance of algorithm has been tested over a large range of SNR.

Future Work and Application
The proposed method using FrFT scheme can overcome the conventional FFT method of target detection and localization that arises in various military and biomedical   applications. Future research in this area will be to apply the concept of FrFT to two dimensional arrays for the joint estimation of azimuth and elevation angles [39]. Also it will be of considerable interest to reduce the effect of mutual coupling as the spacing between the array elements becomes too small. More research in FrFT imaging can lead to better underwater object classification and imaging application as compared to conventional FFT based matched filtering scheme. Additionally, beam-forming with hyperbolic frequency modulated (HFM) chirps having Doppler invariant property, which goes beyond FrFT, is another open research area.