Phased Array System Toolbox

High Resolution Direction of Arrival Estimation

This example illustrates several high resolution direction of arrival (DOA) estimation techniques. It introduces variants of the root MUSIC, ESPRIT and root WSF algorithms and discusses their respective merits in the context of far-field, narrowband signal sources received by a uniform linear array (ULA) antenna.

Modeling the Received Array Signals

Define a uniform linear array (ULA) composed of 10 isotropic antennas. The array element spacing is 0.5 meters.

N = 10;
ha = phased.ULA('NumElements',N,'ElementSpacing',0.5)
ha = 

  System: phased.ULA 

  Properties:
           Element: [1x1 phased.IsotropicAntennaElement]
       NumElements: 10                                  
    ElementSpacing: 0.5                                 
             Taper: 1                                   
                                                        

We simulate the array output for two incident signals. Both signals are incident from 90 degrees in azimuth. Their elevation angles are 73 and 68 degrees respectively. In this example, we assume that these two directions are unknown and need to be estimated. We now simulate the baseband received signal demodulated at operating frequency 300 MHz.

% Model the multichannel received signals at the array
fc = 300e6;                               % Operating frequency
fs = 8192;                                % Sampling frequency
lambda = physconst('LightSpeed')/fc;      % Wavelength
pos = getElementPosition(ha)/lambda;      % Element position in wavelengths
ang1 = [90;73]; ang2 = [90;68];           % Direction of the signals
angs = [ang1 ang2];
Nsamp = 1024;                             % Number of snapshots
noisePwr = 0.01;                          % Noise power

rs = rng(2012);                           % Set random number generator
x = sensorsig(pos,Nsamp,angs,noisePwr);

Because a ULA is symmetric around its axis, a DOA algorithm cannot uniquely determine azimuth and elevation. Therefore, the results returned by these high resolution DOA estimators are in the form of broadside angles. An illustration of broadside angle can be found in the following figure.

Here we calculate the broadside angles corresponding to the two incident angles, which are 17 and 22 degrees:

ang_true = az2broadside(angs(1,:),angs(2,:))
ang_true =

   17.0000   22.0000

Estimating the Direction of Arrival (DOA)

We first assume that we know a priori that there are two sources. To estimate the DOA, we'll use the root MUSIC technique, so we construct a DOA estimator using the root MUSIC algorithm.

hdoaMusic = phased.RootMUSICEstimator('SensorArray',ha,...
            'OperatingFrequency',fc,...
            'NumSignalsSource','Property','NumSignals',2)
hdoaMusic = 

  System: phased.RootMUSICEstimator 

  Properties:
                 SensorArray: [1x1 phased.ULA]
            PropagationSpeed: 299792458       
          OperatingFrequency: 300000000       
            NumSignalsSource: 'Property'      
                  NumSignals: 2               
    ForwardBackwardAveraging: false           
            SpatialSmoothing: 0               
                                              

Since the array response vector of a ULA is conjugate symmetric, we can use forward-backward (FB) averaging to reduce the computational complexity because the computations are all real. These FB based estimators also have a lower variance, and reduce the correlation between signals.

To apply forward-backward averaging, we set the ForwardBackwardAveraging property of the root MUSIC DOA estimator to true. In this case, the root MUSIC algorithm is also called the unitary root MUSIC algorithm.

hdoaMusic.ForwardBackwardAveraging = true;

We can perform the DOA estimation by calling the step method as follows:

ang = step(hdoaMusic,x)
ang =

   16.9960   21.9964

We can also use an ESPRIT DOA estimator. As in the case of root MUSIC, we set the ForwardBackwardAveraging property to true. This algorithm is also called unitary ESPRIT.

hdoaEsprit = phased.ESPRITEstimator('SensorArray',ha,...
             'OperatingFrequency',fc,'ForwardBackwardAveraging',true,...
             'NumSignalsSource','Property','NumSignals',2)
hdoaEsprit = 

  System: phased.ESPRITEstimator 

  Properties:
                 SensorArray: [1x1 phased.ULA]
            PropagationSpeed: 299792458       
          OperatingFrequency: 300000000       
            NumSignalsSource: 'Property'      
                  NumSignals: 2               
            SpatialSmoothing: 0               
                      Method: 'TLS'           
    ForwardBackwardAveraging: true            
                RowWeighting: 1               
                                              
ang = step(hdoaEsprit,x)
ang =

   21.9988   16.9748

Estimating the Number of Signal Sources

In practice, we generally do not know the number of signal sources, and need to estimate the number of sources from the received signal. The number of signal sources can be estimated by the DOA estimator by specifying 'Auto' for the NumSignalsSource property and choosing either 'AIC' or 'MDL' for the NumSignalsMethod property. For AIC, the Akaike information criterion (AIC) is used, and for MDL, the minimum description length (MDL) criterion is used.

Before we set NumSignalsSource, we need to release the DOA object because it is locked to improve efficiency during processing.

release(hdoaEsprit);
hdoaEsprit.NumSignalsSource = 'Auto';
hdoaEsprit.NumSignalsMethod = 'AIC';
ang = step(hdoaEsprit,x)
ang =

   21.9988   16.9748

Reducing Computational Complexity

In addition to the forward-backward averaging mentioned before, there are other methods to reduce computational complexity. One of these approaches is to solve an equivalent problem with reduced dimensions in beamspace. While the ESPRIT algorithm performs the eigenvalue decomposition (EVD) of a 10x10 real matrix in our example, the beamspace version can reduce the problem to the EVD of a 3x3 real matrix. This technique uses a priori knowledge of the sector where the signals are located to position the center of the beam fan. For this case, the beam fan is pointed along 20 degrees in azimuth.

hdoaBSEsprit = phased.BeamspaceESPRITEstimator('SensorArray',ha,...
               'OperatingFrequency',fc,...
               'NumBeamsSource','Property','NumBeams',3,...
               'BeamFanCenter',20);
ang = step(hdoaBSEsprit,x)
ang =

   21.9875   16.9943

Another technique is the root weighted subspace fitting (WSF) algorithm. This algorithm is iterative and the most demanding in terms of processing power. However, you can set the maximum number of iterations by specifying the MaximumIterationCount property to maintain the cost below a specific limit.

hdoaWSF = phased.RootWSFEstimator('SensorArray',ha,...
          'OperatingFrequency',fc,'MaximumIterationCount',2);
ang = step(hdoaWSF,x)
ang =

   21.9962   16.9961

Optimizing Performance

In addition to FB averaging, row weighting can be used to improve the statistical performance of the element-space ESPRIT estimator. Row weighting is a technique that applies different weights on the rows of the signal subspace matrix. The row weighting parameter determines the maximum weight. In most cases, it is chosen as large as possible. However, its value can never be greater than (N-1)/2 where N is the number of elements of the array.

release(hdoaEsprit);
hdoaEsprit.RowWeighting = 4
ang = step(hdoaEsprit,x)
hdoaEsprit = 

  System: phased.ESPRITEstimator 

  Properties:
                 SensorArray: [1x1 phased.ULA]
            PropagationSpeed: 299792458       
          OperatingFrequency: 300000000       
            NumSignalsSource: 'Auto'          
            NumSignalsMethod: 'AIC'           
            SpatialSmoothing: 0               
                      Method: 'TLS'           
    ForwardBackwardAveraging: true            
                RowWeighting: 4               
                                              

ang =

   21.9884   17.0003

Estimating Coherent Sources in Multipath Environments

If several sources are correlated or coherent (as in multipath environments, for example), the spatial covariance matrix becomes rank deficient and certain subspace-based DOA estimation methods fail. We now model a received signal composed of 4 narrowband components. We assume that 2 of the first 3 signals are multipath reflections of the first source, whose magnitudes are a quarter and half of the first source, respectively.

scov = eye(4);
magratio = [1;0.25;0.5];
scov(1:3,1:3) = magratio*magratio';

All signals are in the zero elevation domain, with azimuth incident angles of -23, 0, 12, and 40 degrees.

% Incident azimuth
az_ang = [-23 0 12 40];
% When the elevation is zero, the azimuth within [-90 90] is the same as
% the broadside angle.
el_ang = zeros(1,4);

% The received signals
x = sensorsig(pos,Nsamp,[az_ang; el_ang],noisePwr,scov);
rng(rs);                                 % Restore random number generator

We study the performance of aforementioned DOA algorithms under this coherent sources scenario. To simplify the example, we only give one trial to each algorithm instead of running Monte Carlo simulations. Given the high SNR, the results can be a good indicator of estimation accuracy.

First, we verify that the AIC criterion underestimates the number of sources which causes the unitary ESPRIT algorithm to give incorrect estimates. The AIC estimates the number of sources to be only two because three of them are correlated.

release(hdoaEsprit);
hdoaEsprit.NumSignalsSource = 'Auto';
hdoaEsprit.NumSignalsMethod = 'AIC';
ang = step(hdoaEsprit,x)
ang =

  -15.3535   40.0024

The root WSF algorithm, which is very robust, is particularly interesting in the context of correlated signals. With the correct number of sources as an input, the algorithm does not require any modification.

release(hdoaWSF);
hdoaWSF.NumSignalsSource = 'Property';
hdoaWSF.NumSignals = 4;
ang = step(hdoaWSF,x)
ang =

   40.0016  -22.9919    0.0737   12.0693

ESPRIT and root MUSIC, however, fail to estimate the correct DOA of the different sources, even if we specify the number of sources and use the unitary implementations.

release(hdoaMusic);
hdoaMusic.NumSignalsSource = 'Property';
hdoaMusic.NumSignals = 4;
hdoaMusic.ForwardBackwardAveraging = true;
ang = step(hdoaMusic,x)
ang =

   40.0077  -22.8313    4.4976  -11.9038

We need to apply a spatial smoothing technique to estimate the DOA of the correlated signals. Using spatial smoothing, however, decreases the effective aperture of the array. Therefore, the variance of the estimators increases, since the subarrays are smaller than the original array.

release(hdoaMusic);
Nr = 2;    % Number of multipath reflections
hdoaMusic.SpatialSmoothing = Nr
ang = step(hdoaMusic,x)
hdoaMusic = 

  System: phased.RootMUSICEstimator 

  Properties:
                 SensorArray: [1x1 phased.ULA]
            PropagationSpeed: 299792458       
          OperatingFrequency: 300000000       
            NumSignalsSource: 'Property'      
                  NumSignals: 4               
    ForwardBackwardAveraging: true            
            SpatialSmoothing: 2               
                                              

ang =

   40.0010  -22.9959   12.1376    0.1843

Comparison of DOA Algorithms

In summary, ESPRIT, root MUSIC and root WSF are three important DOA algorithms that provide good performance and reasonable computational complexity for ULAs. Unitary ESPRIT, unitary root MUSIC, and beamspace unitary ESPRIT provide ways to significantly reduce the computational cost of the estimators, while also improving their performance. Root WSF is particularly attractive in the context of correlated sources because, contrary to the other methods, it does not require spatial smoothing to properly estimate the DOAs when the number of sources is known.