## High Resolution Direction of Arrival Estimation

This example illustrates several high-resolution direction of arrival (DOA) estimation techniques. It introduces variants of the MUSIC, 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; ula = phased.ULA('NumElements',N,'ElementSpacing',0.5)```
```ula = phased.ULA with properties: Element: [1x1 phased.IsotropicAntennaElement] NumElements: 10 ElementSpacing: 0.5000 ArrayAxis: 'y' Taper: 1 ```

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

```fc = 300e6; % Operating frequency fs = 8192; % Sampling frequency lambda = physconst('LightSpeed')/fc; % Wavelength pos = getElementPosition(ula)/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 signal = 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 angles can be found in the following figure. Calculate the broadside angles corresponding to the two incident angles.

`ang_true = az2broadside(angs(1,:),angs(2,:))`
```ang_true = 1×2 17.0000 22.0000 ```

The broadside angles are 17° and 22°.

### Estimating the Direction of Arrival (DOA)

Assume that we know a priori that there are two sources. To estimate the DOA, use the root-MUSIC technique. Construct a DOA estimator using the root-MUSIC algorithm.

```rootmusicangle = phased.RootMUSICEstimator('SensorArray',ula,... 'OperatingFrequency',fc,... 'NumSignalsSource','Property','NumSignals',2)```
```rootmusicangle = phased.RootMUSICEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 NumSignalsSource: 'Property' NumSignals: 2 ForwardBackwardAveraging: false SpatialSmoothing: 0 ```

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

To apply forward-backward averaging, 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.

`rootmusicangle.ForwardBackwardAveraging = true;`

Perform the DOA estimation:

`ang = rootmusicangle(signal)`
```ang = 1×2 16.9960 21.9964 ```

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

```espritangle = phased.ESPRITEstimator('SensorArray',ula,... 'OperatingFrequency',fc,'ForwardBackwardAveraging',true,... 'NumSignalsSource','Property','NumSignals',2)```
```espritangle = phased.ESPRITEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 NumSignalsSource: 'Property' NumSignals: 2 SpatialSmoothing: 0 Method: 'TLS' ForwardBackwardAveraging: true RowWeighting: 1 ```
`ang = espritangle(signal)`
```ang = 1×2 21.9988 16.9748 ```

Finally, use the MUSIC DOA estimator. MUSIC also supports forward-backward averaging. Unlike ESPRIT and root-MUSIC, MUSIC computes a spatial spectrum for specified broadside scan angles. Directions of arrival correspond to peaks in the MUSIC spatial spectrum.

```musicangle = phased.MUSICEstimator('SensorArray',ula,... 'OperatingFrequency',fc,'ForwardBackwardAveraging',true,... 'NumSignalsSource','Property','NumSignals',2,... 'DOAOutputPort',true)```
```musicangle = phased.MUSICEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 ForwardBackwardAveraging: true SpatialSmoothing: 0 ScanAngles: [1x181 double] DOAOutputPort: true NumSignals: 2 NumSignalsSource: 'Property' ```
`[~,ang] = musicangle(signal) `
```ang = 1×2 17 22 ```
`plotSpectrum(musicangle)` Directions of arrival for MUSIC are limited to the scan angles in the `ScanAngles` property. Because the true directions of arrival in this example coincide with the search angles in `ScanAngles`, MUSIC provides precise DOA angle estimates. In practice, root-MUSIC provides superior resolution to MUSIC. However, the MUSIC algorithm can also be used for DOA estimation in both azimuth and elevation using a 2-D array. See Direction of Arrival Estimation with Beamscan, MVDR, and MUSIC.

### Estimating the Number of Signal Sources

In practice, you generally do not know the number of signal sources, and need to estimate the number of sources from the received signal. You can estimate the number of signal sources 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 you can set NumSignalsSource, you must release the DOA object because it is locked to improve efficiency during processing.

```release(espritangle); espritangle.NumSignalsSource = 'Auto'; espritangle.NumSignalsMethod = 'AIC'; ang = espritangle(signal)```
```ang = 1×2 21.9988 16.9748 ```

### Reducing Computational Complexity

In addition to forward-backward averaging, other methods can 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. In this example, point the beam fan to 20° in azimuth.

```bsespritangle = phased.BeamspaceESPRITEstimator('SensorArray',ula,... 'OperatingFrequency',fc,... 'NumBeamsSource','Property','NumBeams',3,... 'BeamFanCenter',20); ang = bsespritangle(signal)```
```ang = 1×2 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 computational complexity. You can set the maximum number of iterations by specifying the MaximumIterationCount property to maintain the cost below a specific limit.

```rootwsfangle = phased.RootWSFEstimator('SensorArray',ula,... 'OperatingFrequency',fc,'MaximumIterationCount',2); ang = rootwsfangle(signal)```
```ang = 1×2 16.9961 21.9962 ```

### Optimizing Performance

In addition to FB averaging, you can use row weighting to improve the statistical performance of the element-space ESPRIT estimator. Row weighting is a technique that applies different weights to the rows of the signal subspace matrix. The row weighting parameter determines the maximum weight. In most cases, it is chosen to be 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(espritangle); espritangle.RowWeighting = 4```
```espritangle = phased.ESPRITEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 NumSignalsSource: 'Auto' NumSignalsMethod: 'AIC' SpatialSmoothing: 0 Method: 'TLS' ForwardBackwardAveraging: true RowWeighting: 4 ```
`ang = espritangle(signal)`
```ang = 1×2 21.9884 17.0003 ```

### Estimating Coherent Sources in Multipath Environments

If several sources are correlated or coherent (as in multipath environments), the spatial covariance matrix becomes rank deficient and subspace-based DOA estimation methods may fail. To show this, model a received signal composed of 4 narrowband components. Assume that 2 of the first 3 signals are multipath reflections of the first source, having magnitudes equal to 1/4 and 1/2 that of the first source, respectively.

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

All signals are incident at 0° elevation, with azimuth incident angles of -23°, 0°, 12°, and 40°.

```% 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 signal = sensorsig(pos,Nsamp,[az_ang; el_ang],noisePwr,scov); rng(rs); % Restore random number generator```

Compare the performance of the DOA algorithm when sources are coherent. To simplify the example, run only one trial per algorithm. Given the high SNR, the results will be a good indicator of estimation accuracy.

First, verify that the AIC criterion underestimates the number of sources, causing the unitary ESPRIT algorithm to give incorrect estimates. The AIC estimates the number of sources as two because three sources are correlated.

```release(espritangle); espritangle.NumSignalsSource = 'Auto'; espritangle.NumSignalsMethod = 'AIC'; ang = espritangle(signal)```
```ang = 1×2 -15.3535 40.0024 ```

The root-WSF algorithm is robust in the context of correlated signals. With the correct number of sources as an input, the algorithm correctly estimates the directions of arrival.

```release(rootwsfangle); rootwsfangle.NumSignalsSource = 'Property'; rootwsfangle.NumSignals = 4; ang = rootwsfangle(signal)```
```ang = 1×4 40.0016 -22.9919 12.0693 0.0737 ```

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

```release(rootmusicangle); rootmusicangle.NumSignalsSource = 'Property'; rootmusicangle.NumSignals = 4; rootmusicangle.ForwardBackwardAveraging = true; ang = rootmusicangle(signal)```
```ang = 1×4 40.0077 -22.8313 4.4976 -11.9038 ```

You can apply spatial smoothing to estimate the DOAs of correlated signals. Using spatial smoothing, however, decreases the effective aperture of the array. Therefore, the variance of the estimators increases because the subarrays are smaller than the original array.

```release(rootmusicangle); Nr = 2; % Number of multipath reflections rootmusicangle.SpatialSmoothing = Nr```
```rootmusicangle = phased.RootMUSICEstimator with properties: SensorArray: [1x1 phased.ULA] PropagationSpeed: 299792458 OperatingFrequency: 300000000 NumSignalsSource: 'Property' NumSignals: 4 ForwardBackwardAveraging: true SpatialSmoothing: 2 ```
`ang = rootmusicangle(signal)`
```ang = 1×4 40.0010 -22.9959 12.1376 0.1843 ```

### Comparison of DOA Algorithms

In summary, ESPRIT, MUSIC, root-MUSIC, and root-WSF are 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.