An Estimation of Angular Coordinates and Angular Distance of Two Moving Objects by Using the Monopulse, MUSIC, and Root-MUSIC Methods

The paper presents results of extensive simulations carried out in order to assess the precision and angular resolution of subspace methods in real radar system. It has been assumed that such a system uses the 32-element uniform linear array (ULA) and radiates only 3 bursts consisting of 8 pulses in given direction. In order to avoid blind Doppler frequencies, pulse repetition interval (PRI) is different in each burst. It has been shown that change of PRI is not only necessary to avoid blind Doppler frequencies but also allows to avoid false values of angular coordinates when two objects are visible in the same beam, in the same range gate, and their echoes attain maximal values in the same Doppler filter. It has also been shown that precision and angular resolution of both MUSIC and root-MUSIC method can be improved by appropriate preprocessing of signal samples used by these methods.


Introduction
The paper studies the problem of detection and estimation of angular coordinates of moving objects by means of MUSIC and root-MUSIC methods. MUSIC and root-MUSIC had been invented almost 30 years ago [1,2]; however, application of these methods in real radar system has been possible only recently thanks to the progress in the field of active electronically scanned arrays, digital beamforming, and multiprocessor systems containing clusters of high-speed general purpose PowerPC processors and FPGA devices connected by means of high-speed serial bus such as RapidIO. An important reason stimulating the studies on these methods were difficulties concerning effective estimation of angular coordinates of closely spaced moving objects by using the monopulse methods, [3][4][5]. For example, it is not possible to separate individual objects even with different radial speeds if they are illuminated by the same antenna beam and are in the same range gate. More appropriate to solve problems of this kind are superresolution methods, represented here by MUSIC and root-MUSIC, [6,7]. In application to estimation of angular coordinates of moving objects, these methods use the spatial correlation (covariance) matrix R formulated on the basis of complex samples of received signals. As a rule, these samples are obtained on the outputs of individual matched filters (MF) of the array antenna unit, which have functional scheme similar to that shown in Figure 1.
All computer simulations presented in this paper have been performed under the assumption that this linear equidistant array antenna contains M = 32 identical receiving elements spaced by the distance d = 0.7λ 0 , where λ 0 is the length of received wave.
This 32-element array antenna has been used to create new 24-element array. Samples x i (n), where 1 ≤ i ≤ 32, of signals received by the real, 32-element array antenna are grouped in the following way: where w l = exp[ jl(2πd/λ) sin(θ max )] is the lth coefficient of vector w shaping the directive gain characteristics of antenna elements, belonging to the 24-element array, and θ max determines the direction of maximum directivity of each antenna element. This kind of grouping is often called initial preprocessing. Samples x (e) k (n), where 1 ≤ k ≤ 24, corresponding to individual elements of the equivalent, 24element array are subsequently processed using the MUSIC and root-MUSIC methods. It has been assumed that the antenna array described above radiates 24 pulses, grouped into three 8-pulse bursts differing by the repetition time (frequency), as shown in Figure 2.
The signals (radar echoes) received at the outputs of elements of the array are used for digital beamforming, s(n) = w H sum · x(n) and are processed later according to the MTD (moving target detection) technique, [8]. Echo signals s(n) received after each pulse are subjected to coherent integration using predefined sets of weight coefficients w k where b is the number of range gate determining the distance between the object and the radar, B is the number of range gates in each scan, and k is the number of pulse belonging to a given pulse sequence called the burst. This integration is called Doppler filtration, because it permits to distinguish between the echo signals from the objects moving with different radial speeds. The signal samples from individual array elements, weighted according to (1), form the vector x (e) that is the basis for calculation of the following estimate: of the correlation matrix R, where x (e)H is the Hermitian conjugate with respect to the vector x (e) . In practical considerations, the factor 1/N appearing in (3) is usually neglected, because it has no influence on the values of the used eigenvectors, but only decides about the scale of the related eigenvalues of matrix R [2,9].

The Multiple Signal Classification (Music) Method
As it was mentioned in the introduction, the MUSIC method belongs to the group of subspace methods in which the ISRN Signal Processing spectral functions are determined on the base of eigenvectors of the space correlation matrix R. Similarly, the matrix R is formulated on the base of complex samples of the received signals, [1,6,10]. In order to explain the essence of this method in the radar application, assume that the Nelement antenna array receives P echo signals, reflected from P objects, where P < N. In a given moment, to each echo of the radiated pulse one can assign the vector of complex samples of the signals received in this moment by individual elements of the array antenna; see Figure 1. Therefore, we assume that to the first element of this array, from the direction θ k , where 1 ≤ k ≤ P, arrives the signal s 1 (t) = A cos(ωt + φ k ) represented in further considerations by its complex amplitude s 1 = A exp( jφ k ). It is easy to prove on the base of Figure 1 that complex amplitude of the signal coming from the same direction θ k to the i-element of the array can be written as follows: c ≈ 3 · 10 8 m/s is the speed of light, and f denotes frequency of the narrow-band received signal. The complex samples of this signal, appearing on the outputs of individual matched filters MF, will be shifted in phase in a similar way. The set of these samples constitutes a N -dimensional column vector . . .
where v(θ k ) = 1, e − jβd sin(θk) , e − jβ2d sin(θk) , . . . , e − j(N−1)βd sin(θk) T , is the so-called array steering vector, [10]. Each element of the array shown in Figure 1 receives simultaneously the signals from all directions θ k , where 1 ≤ k ≤ P, and noise signal with variance σ 2 n . It means that the resulting column vector of complex samples of all P signals and noise, before further digital processing, can be written as where is the N × P matrix, the columns of which are the elements of direction vectors (5), s = [s 1 (θ 1 ), s 1 (θ 2 ), s 1 (θ 3 ), . . . , s 1 (θ P )] T is the P-element, column vector representing the signals received by the first element of the array, and n is the Nelement column vector representing received noise signals. The correlation matrix R of samples of signals obtained on the outputs of antenna system, see Figure 1, can be written in the form of the sum, using the correlation matrix of desirable signals R s and the noise correlation matrix R n , [1,10], namely, According to [1], the correlation matrix of desirable signals is the diagonal matrix P × P, if the received desirable signals are not correlated. Under the assumption P < N, matrix VR s V H is the singular matrix, which means that It follows from (10) that σ 2 n is the eigenvalue of the matrix R, [9]. Space, in which the desirable signals are not defined has the dimension N − P. Hence, σ 2 n is (N − P)-order eigenvalue of the matrix R. The matrices R and VR s V H are not negative defined, and therefore, the matrix R has also P other eigenvalues λ k , satisfying condition λ k > σ 2 n > 0, where 1 ≤ k ≤ P. The eigenvectors q k , assigned to these eigenvalues are mutually orthogonal [9,10]. According to the general definition of matrix eigenvalues problem, (R − λ k )q k = 0, we can write Relation (12) shows that N -dimensional space of signals and noise can be divided into two mutually orthogonal subspaces, that is, subspace of signals Q s ≡ [q 1 , q 2 , q 3 , . . . , q P ] and subspace of noise According to this partition, the correlation matrix R can be written as the following sum: where λ k denotes the eigenvalue corresponding to the vector q k and σ 2 n is the variance of the white noise received from the individual antenna elements. After appropriate grouping of the factors of sum (13), we obtain the relation, in which one can distinguish P eigenvectors representing desirable signals and N − P eigenvectors belonging to the noise subspace Each of vectors (4) containing complex samples of desirable signal belongs to the signal subspace Q s , and therefore, it can be written in the form of sum of eigenvectors, defined in this subspace, namely, where b k is the k coefficient of a suitable value. It should be pointed out that each component q k of vector (15) is orthogonal with respect to an arbitrary eigenvector q m from the noise subspace Q n . Consequently, the whole vector (15) is orthogonal to q m .This unique property can be expressed as follows: Applying equation (16) to all eigenvectors of the noise subspace Q n , we find that the dot product of the vector v(θ k ), representing the signal received from direction θ k , and the sum of eigenvectors from the noise subspace Q n , will also take the value close to zero. In the ideal case (exactly defined correlation matrix R, exactly evaluated eigenvalues and their corresponding eigenvectors, and precise partition onto the vectors from the subspaces of signal and noise), we have Using the property described by (17), the following estimate of the spectral power density of the signal can be formulated: This estimate is usually called the spectrum of the MUSIC method, [1,6]. Placing all eigenvectors q k from the noise subspace into the columns of matrix Q n , spectrum (18) can be written in the equivalent, simpler form Function P(θ) attains local maximum values for angles θ k determining directions of arrival of signals being received.

The Root-Music Method
Determination of angular positions θ k on the base of spectrum (18) or (19), requires performing calculations for great number of discrete values of angle θ and next determination of all its maximum values in the given, relatively large scanning range θ min ≤ θ ≤ θ max . This task is especially laborious and time consuming when the angular resolution of the order of one tenth or one hundreds of degree is required. Therefore, in order to reduce the amount of calculations, the modified version of the MUSIC method, called root-MUSIC, has been elaborated. In this improved version, the problem of evaluation of the local maximum values of function (19) is replaced by the problem of finding the roots . Estimated values of angular coordinates of objects can be evaluated for the assumed number of roots that should be equal to the number of received desirable signals multiplied by 2. This number is usually determined by means of special criteria, among which the most known are AIC (akaike information criterion) and MDL (minimum description length), [11].
Denominator of function (19) is in general a polynomial, which can be written as is the Hermitian matrix of degree N . Coefficients c n of the polynomial (20) can be determined by summing elements p kl of matrix P placed on its nth diagonals, namely, According to this formula, As was mentioned earlier, P is the Hermitian matrix, and therefore, the coefficient of polynomial (20) with indexes −n and n are mutually conjugate, so c −n = c * n for 1 ≤ n ≤ N − 1. Equation C(z) = 0 is of degree 2N − 1 and has 2N − 1 roots, and to each root z n , there corresponds another root 1/z * n . These roots are most frequently evaluated by means of the companion matrix method, [12,13]. It follows from the literature that this method ensures sufficient accuracy of calculations for all roots, even when polynomial C(z) = 0 is of relatively high degree; for instance, 2N − 1 = 63. Due to this valuable property, the companion matrix method has been implemented in the computational environment MATLAB.
Estimates of P angular positions of objects being detected are evaluated on a basis of 2P roots situated nearest to the unitary circle determined on the complex plane z = Re(z) + j Im(z), namely, on the basis of P pairs (z n , 1/z * n ). Of course, each pair chosen in this manner determines only one location. With negligibly small power of the noise, σ 2 n ≈ 0, the roots lay exactly on the unitary circle mentioned above. From the substitution z = exp[− jβd sin(θ)] introduced above, it follows that estimates of angular positions θ n are where 1 ≤ n ≤ P, β = 2π/λ and d is the antenna element spacing.

Application of Music and Root-Music Methods to the Estimation of Angular Coordinates of Moving Objects
As was already mentioned in the introduction, exact estimation of the angular coordinates of moving objects by means of monopulse methods is in many cases impossible because of their limited angular resolution. For this reason, they cannot distinguish the objects illuminated by the same antenna beam and situated at the same slant distance. In order to compare angular resolution of the monopulse, MUSIC and root-MUSIC methods, some computer simulations have been carried out. In these simulations, angular coordinates of two objects (planes) moving with the speed of v = 100 m/s have been evaluated. It has been assumed also that both planes are moving at 45 • with respect to the north direction. The position of the first plane is determined by constant angle θ 1 = −14.1 • , while the second plane changes its position θ 2 gradually from −14.1 • to −10.1 • preserving the course and speed. In other words, angular separation Δθ = |θ 1 − θ 2 | between these planes (objects) changes in the range of Δθ = 0 • ÷ 4 • . A limited number of pulses radiated by the radar in given direction and consequently limited number of signals vectors x (e) cause that estimate of correlation matrix R is inaccurate and have to be modified using diagonal loading technique before eigendecomposition and estimation of angular coordinates by means of subspace methods where δ = 4σ 2 n − loading factor. The value of mean pulse repetition interval (mean PRI) has been set to T p = T p−mean = 2 ms. The real PRI in ith burst is equal to 0.9 · T p−mean ≤ T pi ≤ T p−mean and is modified in each burst in order to avoid blind Doppler frequencies (T p1 = 1.0 · T p−mean , T p2 = 0.95 · T p−mean , T p2 = 0.9 · T p−mean ).
The mean square errors of estimates of angular coordinates of both objects obtained by means of the monopulse method are illustrated in Figure 3.
As it is seen, the second plane changing its angular position causes estimation errors of angular coordinates of the first of them. This effect can be eliminated in some extent by means of the additional Doppler filtration. This filtration allows to attenuate echo signals coming from the second plane when both signals have substantially different Doppler frequency shifts normalized with respect to pulse repetition frequency (PRF) level located at −18 dB. Therefore, if echo of the first object attain its maximal value in nth Doppler filter, Doppler frequency of echo of the second object is not located in the stopband of this filter and cannot be sufficiently attenuated. The difference of Doppler frequencies Δ f d is proportional to the speed of both objects, and at some point, echoes of these objects can be separated using Doppler filtration. Unfortunately, at higher speeds (v = 1000 m/s) these echos may have similar normalized Doppler frequencies f d−norm = f d1 %F p = f d2 %F p despite the fact that their Doppler frequencies are different f d1 = f d−norm + n · F p , f d2 = f d−norm + m · F p . This effect causes errors of estimation of angular coordinates of the objects being detected.
The disadvantageous effect under discussion is well illustrated by the simulation results shown in Figure 4. These simulations have been performed for the planes, the speed of which has been increased to v = 1000 m/s. The pulse repetition interval (PRI) has been also increased to T p = 4 ms in order to lower the PRF and increase number of similar normalized Doppler frequencies.
From the results presented above, it follows that the additional Doppler filtration is not an universal solution ensuring proper distinction of the moving objects at all times.
The signal to noise ratio S/N noise , given in Figures 3 and  4, determines the value of this parameter on the outputs of matched filters MF, see Figure 1. The power ratio of desired signal and jamming signal S/N jamm is calculated for the inputs of pulse compression blocks. This ratio has been defined before digital compression, because the noise signal can be differently correlated to the radiated LFM pulse. The second, very important advantage of MUSIC and root-MUSIC methods is their sufficient immunity to relatively high power jamming signals. In order to confirm this conclusion, the additional narrowband strong signal, situated in the direction defined by θ jamm = 1.4 • has been introduced. The surface power density of this signal is 60 dB greater than the corresponding surface power density of both  desirable radar signals. In this situation, the radar systems with receiving antennas nonadapted to the interferences, will be jammed and eventual detections may have improper estimates of angular coordinates. Mean square errors of such improper estimates of angular coordinates are illustrated in Figure 5. Of course, a negative influence of the single strong jamming signal can be decreased by using the receiving antenna in a form of multielement array, which can be adapted to this undesirable signal, [14,15].
In other words, the jamming signal should be attenuated by the array antenna in the highest degree. The simulation results presented in Figure 6 confirm an effectiveness of application of the above approach to solve the radiolocation problem under consideration.  Next, the same radiolocation problem has been solved by using the MUSIC and root-MUSIC methods. Figures 7 and 9 show mean square errors of estimates obtained by means of MUSIC method for S/N noise = 0 dB and S/N noise = 30 dB in presence of jamming signal determined above.
Similarly, the corresponding mean square errors of estimates evaluated using the root-MUSIC method for S/N noise = 0 dB, S/N noise = 30 dB and the same jamming signal are shown in Figures 8 and 10.
The values of angular resolution evaluated for these methods are given correspondingly in Tables 1 and 2. For comparison, Tables 3 and 4 contain values of angular resolution obtained, when MUSIC and root-MUSIC methods determine angular coordinates on the basis of  In other words, in this process, the stage of initial preprocessing, see introduction, has been omitted. Comparing results given in Tables 1 and 2, as well in Tables 3 and  4, we can see that using the proposed initial preprocessing, see relation (1), gives significant amelioration of precision and angular resolution. According to these results, using of initial preprocessing permits to reduce the S/N noise ratio for about 6 dB conserving necessary precision and resolution. Thus, the error magnitudes given in Tables 1 and 2 seem to be acceptable for most similar radiolocation problems encountered in practice. As it has been mentioned in the beginnig of the paper, all simulations have been carried out 8 ISRN Signal Processing assuming that PRI is changing in each burst. The following values of PRI have been used: Estimation of correlation matrix on basis of signals vectors received for 3 different intervals allows MUSIC and root-MUSIC methods to estimate the correct values of angular coordinates despite the fact that signals could have the same or very close normalized Doppler frequency f d−norm = f d %F p and be correlated for given PRI. This technique is well known in radar literature, but it has been mainly used to avoid blind speeds. Presented simulation show, that it also helps to mitigate the problem of estimation of angular coordinates of closely spaced highly correlated signals. When objects are moving with very high velocities (v = 1000 m/s) their angular coordinates could be false when emitted signals have the same PRI in each burst as shown in Figures 11 and 12.
This effect is especially apparent when PRF is relatively small for instance when F p = 250 Hz (T p = 4 ms). The values of normalized Doppler frequencies and mean square errors of angular coordinates obtained for this case are illustrated in Figures 13 and 14.
The further study of this problem and comparison of these approach with well-known techniques of decoration of signals such as spatial smoothing or redundancy averaging is behind the scope of this paper.

Conclusions
The radiolocation problem defined at the beginning of the introduction and in Section 4 has been subsequently solved by different methods, that is, amplitude monopulse  method, amplitude monopulse method aided by the coherent Doppler filtration, MUSIC, and root-MUSIC. Diagrams shown in Figures 3 and 10, respectively, illustrate the results of computer simulations obtained by means of these methods. Thus, on the basis of diagrams shown in Figure 3, we can deduce that traditional amplitude monopulse method is inadequate for this purpose. Some amelioration can be obtained using the additional coherent Doppler filtration. Unfortunately, this solution is not always effective, because false estimates can appear in cases when echo signals after the Doppler filtration [8] attain similar values in the same frequency channel or have similar normalized Doppler frequencies. This effect is illustrated on diagram shown in Figure 4. Partial elimination of this undesirable effect is possible using multiburst radar signals with variable pulse repetition time (frequency), that is, similar to that shown in Figure 2. The  most effective and sufficiently precise for these applications, see Tables 1 and 2, proved to be MUSIC and root-MUSIC methods. This conclusion is well justified by the corresponding simulation results illustrated in Figures 7 and 10.
The radiolocation problem under consideration becomes especially difficult to solve after introducing the narrowband strong jamming interference. In the paper, it has been assumed that this jamming signal is incoming from the direction θ jamm = 1.4 • , and its power density is 60 dB greater than the corresponding surface power density of both desirable radar signals. The influence of this jamming signal on results of the similar simulations performed by using the MUSIC and root-MUSIC methods is illustrated by diagrams shown in Figures 7 and 10.
The presented results confirm that MUSIC and root-MUSIC methods are suitable for effective solution of radiolocation problems similar to that discussed here. Moreover,  they show that these methods can be treated as reliable for the radar signals, for which the power signal-noise ratio is not less than the signal-noise ratio, at which the signal surpasses a detection threshold. It has been also confirmed that the proposed initial preprocessing, see (1), makes possible significant amelioration of precision and angular resolution of MUSIC and root-MUSIC methods. Consequently, they  can be applied for relatively weak radar signals, for which S/N noise ≈ 0 dB. All simulation results presented in this paper have been obtained under the assumption that the number P of objects being detected is known exactly. This number is required to the appropriate partition of space, spanned on the correlation matrix R, into the useful signal subspace Q s and the noise subspace Q n . Of course, is not easy to determine this number for the majority of quasireal radiolocation scenarios, and for this reason, it can be a source of potential errors. Thus, the additional algorithm determining this number with highest possible precission, according to AIC or MDL criterion, is required.

Summary
The main subject of considerations are MUSIC and root-MUSIC methods used to estimation of the angular coordinates (directions of arrival) and angular distance of two moving objects in presence of uncorrelated noise signal and an external, relatively strong narrowband jamming interference. At the receiving antenna, the 32-element uniform linear array (ULA) is used. Extensive computer simulations have been carried out to demonstrate the sufficient accuracy and good spatial resolution properties of these methods in the scenario defined above. It is also shown that using the proposed initial preprocessing, we can increase the accuracy and angular resolution of the methods under discussion. Most of simulation results, presented mainly in a graphical form, have been compared with the corresponding simulation results obtained by using the traditional amplitude monopulse method and the amplitude monopulse method aided by the coherent Doppler filtration.