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.

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 = phased.ULA with properties: Element: [1x1 phased.IsotropicAntennaElement] NumElements: 10 ElementSpacing: 0.5000 ArrayAxis: 'y' 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

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 = phased.RootMUSICEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 NumSignalsSource: 'Property' NumSignals: 2 ForwardBackwardAveraging: 0 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 = phased.ESPRITEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 NumSignalsSource: 'Property' NumSignals: 2 SpatialSmoothing: 0 Method: 'TLS' ForwardBackwardAveraging: 1 RowWeighting: 1

ang = step(hdoaEsprit,x)

ang = 21.9988 16.9748

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

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

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 = phased.ESPRITEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 NumSignalsSource: 'Auto' NumSignalsMethod: 'AIC' SpatialSmoothing: 0 Method: 'TLS' ForwardBackwardAveraging: 1 RowWeighting: 4 ang = 21.9884 17.0003

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 = phased.RootMUSICEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 NumSignalsSource: 'Property' NumSignals: 4 ForwardBackwardAveraging: 1 SpatialSmoothing: 2 ang = 40.0010 -22.9959 12.1376 0.1843

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.

Was this topic helpful?