Mathematical Analysis and Improvement of the Maximum Spatial Eigenfilter for Direction of Arrival Estimation

Abstract

Maximum spatial eigenfiltering improves the accuracy of maximum likelihood direction-of-arrival estimators for closely-spaced signal sources but may interchangeably attenuate widely-spaced signal sources, producing a severe performance degradation. Although this behavior has been observed experimentally, it still lacks a mathematical explanation. In our previous work, we overcame these limitations using a differential spectrum-based spatial filter but this still caused a small degradation in the DOA estimate. In this paper, we develop a mathematical analysis of how the signal source separation and the Karhunen-Loève expansion affect the passbands of the maximum spatial eigenfilter. The farther the sources, the less significant is the maximum eigenvalue of the spatial correlation matrix and its corresponding eigenvector. Then, the magnitude response of the maximum spatial eigenfilter no longer approximates the spatial power spectrum and is not guaranteed to place multiple passbands around the signal sources. Consequently, we propose a spatial filter built from the eigenvectors of the entire signal subspace. This filter showed an overall runtime smaller than that of our previous work. It also provides a significant reduction in the threshold signal-to-noise ratio for closely-spaced signal sources and does not hamper the estimation for widely-spaced signal sources.

Index Terms
Direction of arrival; Maximum likelihood estimation; Noise reduction; Spatial filtering

I. INTRODUCTION

Localization techniques enable several types of state-of-the-art commercial and government applications and services [¹]. These techniques are regarded as a core technology in fourth and fifth generation (4G/5G) wireless telecommunications systems as they allow forming beams to individual mobile users in a multiple-input multiple-output (MIMO) urban environment [²], increasing spectral efficiency and the capacity of those systems [³]. They are also useful to estimate the position of indoor wireless users [¹], for target localization and imaging using MIMO radars [⁴], and even for localization of long-distance underwater acoustic sources [⁵].

In this context, direction-of-arrival (DOA) estimation plays a key role, since all these applications need to know a priori the signal source position. Currently, there are several DOA estimators but the most representative ones are: Multiple Signal Classification (MUSIC) [⁶], Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) [⁷], and Method of Direction Estimation (MODE) [⁸].

MODE has advantages to MUSIC and ESPRIT, making it a competitive DOA estimator, such as maximum likelihood (ML) estimation performance for a large number of snapshot data, but without requiring a computationally-intensive search for the global minimum. It produces DOA estimates through polynomial root finding and does not show convergence problems as in iterative methods. Also, differently from MUSIC and ESPRIT, MODE can handle both uncorrelated and correlated signals [⁹].

Additionally, a class of MODE-based estimators called Principal-eigenvector Utilization for Modal Analysis (PUMA) [¹⁰] and Enhanced PUMA (EPUMA) [¹¹] has been proposed in the literature. However, according to [¹²], the PUMA criterion is exactly equivalent to the MODE criterion and the algorithm proposed for EPUMA corresponds to MODE with Extra Roots (MODEX) [¹³].

MODEX was proposed to overcome the threshold breakdown effect present in MODE, which is a serious performance degradation due to low signal-to-noise ratio (SNR) and/or insufficient snapshot data. This is done by increasing the number of generated estimates and employing an additional ML procedure to select the best set among them, what improves the threshold performance. Although, some of the extra roots produced by MODEX are associated with the noise subspace and can lead to wrong results. So, in order to produce better estimates, Modified MODEX [¹⁴] generates extra roots by computing MODE three times, each one with a distinct non-triviality constraint.

Another way to improve even further the threshold performance is by filtering the data in the ML procedure, such as Krummenauer et al. proposed in [¹⁵]. They achieved a significant reduction on the threshold SNR of the conventional Modified MODEX for closely-spaced signal sources by employing the maximum spatial eigenfilter [¹⁶], [¹⁷] during the ML procedure.

On the other hand, we observed that that filter may significantly attenuate any widely-spaced signal source [¹⁸], leading to a severe performance degradation as the signal source separation increases. This separation relates with the eigendecomposition of the spatial correlation matrix so that the closer the signal sources, the larger the difference between its two largest eigenvalues [¹⁹]–[²¹].

Then, we proposed in [¹⁸] a moving average multiband filter based on the differential spectrum [¹⁹]–[²¹] to overcome the maximum spatial eigenfilter limitations. This filter improved the DOA estimation performance for closely-spaced signal sources but it was capable of preserving signal passbands also for widely-spaced signal sources, performing significantly better than the maximum spatial eigenfilter. Even so, it still performed slightly worse than the conventional Modified MODEX in the latter case.

Thus, in the present paper, we provide a mathematical demonstration of how signal source separation affects the maximum spatial eigenfilter passbands in light of the Karhunen-Loève expansion [¹⁷], [²²] to solve definitely this problem. As a consequence, we derive a new, less-time-consuming, multiband spatial filter that preserves the passbands around the DOA of every signal source, regardless of whether they are closely or widely-spaced from each other.

This work is organized as follows: Section II sets up the signal model for DOA estimation, provides an overview on the conventional Modified MODEX, and depicts the maximum spatial eigenfilter. Section III analyses the filter behavior, and Section IV presents our new proposition. Section V brings the experimental results, and Section VI draws conclusions and suggestions for future works.

II. Theoretical Framework

Let us consider M narrowband far-field signal sources impinging at DOAs θ_m, m = 1,…,M, on a uniform linear array (ULA) formed by K half-wavelength-spaced sensors, with K > M N. The set of N snapshots of the array output is modeled as [¹³]:

(1) $$\mathbf{y}\left(n\right)=\mathbf{A}\mathbf{s}\left(n\right)+\mathbf{n}\left(n\right),$$

for n = 1,…,N; y(n) ∈ ℂ^K×1 is the noisy output vector; s(n) ∈ ℂ^K×1 is the vector containing all M signals and n(n) ∈ ℂ^K×1 is the additive white gaussian noise vector. Both s(n) and n(n) consider the unconditional model [²³]; A = [a(ω₁),…, a(ω_M)] is the K × M matrix of steering vectors $\mathbf{a}\left({\omega }_m\right)={\left[1,\hspace{0.17em}{e}^{-j{\omega }_m},\dots ,\hspace{0.17em}{e}^{-j\left(K\hspace{0.17em}-\hspace{0.17em}1\right){\omega }_m}\right]}^{\mathrm{T}}$ for each one of the M frequencies ω_m = π sin θ_m.

The spatial correlation matrix R from the array output is given by R = ACA^H + σ² I, where I is an identity matrix; (·)^H is the conjugate-transpose operator; σ² is the unknown noise power. C ∈ ℝ^M×M is the signal correlation matrix. Since the signal sources are considered to be equally-powered and uncorrelated to each other, C = I. Actually, R can be estimated and eigendecomposed as [¹⁴]:

(2) $$\widehat{R}=\frac{1}{N}{\displaystyle \sum }_{n=1}^N\mathbf{y}\left(n\right){\mathbf{y}}^{\mathrm{H}}\left(n\right)={\displaystyle \sum }_{k=1}^K{\widehat{\lambda }}_k{\widehat{q}}_k{\widehat{q}}_k^{\mathrm{H}},$$

where ${\widehat{\lambda }}_k$ are its eigenvalues arranged in descending order of magnitude and ${\widehat{q}}_k$ are their corresponding eigenvectors. The first $\overline{M}=\min \left[M,\hspace{0.17em}\text{rank}\left(\mathbf{C}\right)\right]$ eigenvalues along with ${\widehat{q}}_1,\hspace{0.17em}\dots ,\hspace{0.17em}{\widehat{q}}_{\overline{M}}$ span the signal subspace of R. The remaining $\left(K-\overline{M}\right)$ eigenvalues and eigenvectors span its orthogonal or noise subspace [²⁴].

The conventional Modified MODEX uses MODE to compute three solutions to the DOA problem in order to keep good asymptotic estimation performance [¹⁴]. The first one is calculated from the complete MODE algorithm, estimating the roots $\left\{e^{j{\omega }_1},\dots ,e^{j{\omega }_{\overline{M}}}\right\}$ of the polynomial $b\left(z\right)=b_0{z}^{\overline{M}}+b_1{z}^{\left(\overline{M}-1\right)}+\cdots +b_{\overline{M}}$ by solving the following non-linear optimization problem:

(3) $$\mathbf{b}=\text{arg\hspace{0.17em}min\hspace{0.17em}tr}\left[\widehat{B}{\left({\widehat{B}}^{\mathrm{H}}\widehat{B}\right)}^{-1}{\widehat{B}}^{\mathrm{H}}{\widehat{Q}}_S\mathbf{W}{\widehat{Q}}_S^{\mathrm{H}}\right],$$

where $\mathbf{b}=\left[b_0,\hspace{0.17em}\dots ,b_{\overline{M}}\right]$ contains the polynomial coefficients and is subject to the unit-norm constraint b ≠ 0 to avoid a trivial solution, and to the conjugate-symmetry constraint $b_{\overline{M}}=b_{\overline{M}-m}^{*},m=1,\dots ,\overline{M}$, where (·)* is the complex conjugate operator [¹⁴]; ${\widehat{Q}}_S=\left[{\widehat{q}}_{\overline{1}},\dots ,\hspace{0.17em}{\widehat{q}}_{\overline{M}}\right]$; $\mathbf{W}=\left({\widehat{\Lambda }}_S-{\widehat{\sigma }}^2\mathbf{I}\right){\widehat{\Lambda }}_S^{-1}$, ${\widehat{\sigma }}^2=\text{tr}\left({\widehat{\Lambda }}_N\right)/(K-\overline{M})$; ${\widehat{\Lambda }}_S=\text{diag}\left({\widehat{\lambda }}_1,\dots ,{\widehat{\lambda }}_{\overline{M}}\right)$; ${\widehat{\Lambda }}_N=\text{diag}\left({\widehat{\lambda }}_{\left(\overline{M}+1\right)},\dots ,{\widehat{\lambda }}_K\right)$ and:

(4) $$\mathbf{B}={\left[\begin{array}{cccccc}{b}_{\overline{M}}& \dots & b_1& b_{\mathbf{0}}& \hspace{0.17em}& 0\\ \hspace{0.17em}& \ddots & \hspace{0.17em}& \ddots & \ddots & \hspace{0.17em}\\ \mathbf{0}& \hspace{0.17em}& b_{\overline{M}}& \dots & b_1& b_0\end{array}\right]}^{\mathrm{H}}\in {ℂ}^{K\times \left(K-\overline{M}\right)},$$

such that B^HA = 0. The other two solutions come respectively from the imposition of the real linear constraint ℜ(b₀) = 1 and imaginary linear constraint $\Im \left(b_{\overline{M}}\right)=1$, both computed using the first solution as a parameter.

However, this procedure generates $3\overline{M}$ estimates. Thus, an ML procedure must take place to combine those estimates into M-tuples that are included in the set Ω to select the best one among them:

(5) $$\widehat{\omega }=\mathrm{arg\hspace{0.17em}}\underset{\mathbf{\Omega }}{\min }\text{\hspace{0.17em}tr}\left\{\left[\mathbf{I}-\widehat{A}{\left({\widehat{A}}^{\mathrm{H}}\widehat{A}\right)}^{-1}{\widehat{A}}^{\mathrm{H}}\right]\widehat{R}\right\},$$

where $\widehat{A}$ is a candidate matrix of steering vectors for each M-tuple in Ω.

Even though the noise will not affect distinct MODE solutions simultaneously and/or the same estimate in each solution [¹⁴], the threshold breakdown effect may occur as the SNR gets below a critical value [¹⁵]. Then, the noise present in $\widehat{R}$ worsens the ML procedure performance in Eq. (5) such that the conventional Modified MODEX may not choose the best M-tuple [¹⁵].

Since the maximum eigenfilter, i.e. the eigenvector associated with the largest eigenvalue of $\widehat{R}$ [¹⁷], maximizes the output SNR, Krummenauer et al. [¹⁵] proposed using it as a spatial filter [¹⁶] to reduce the noise influence in Eq. (5):

(6) $$\tilde{\omega }=\mathrm{arg\hspace{0.17em}}\underset{\mathbf{\Omega }}{\min }\text{\hspace{0.17em}tr}\left\{\left[\mathbf{I}-\tilde{A}{\left({\tilde{A}}^{\mathrm{H}}\tilde{A}\right)}^{-1}{\tilde{A}}^{\mathrm{H}}\right]\tilde{R}\right\},$$

where $\tilde{A}=H\widehat{A}$; $\tilde{R}=H\widehat{R}{\mathbf{H}}^{\mathrm{H}}$ and H is the following (L + K) × K convolution matrix:

(7) $$\mathbf{H}={\left[\begin{array}{cccccc}{h}_L& \dots & h_1& h_0& \hspace{0.17em}& \mathbf{0}\\ \hspace{0.17em}& \ddots & \hspace{0.17em}& \ddots & \ddots & \hspace{0.17em}\\ \mathbf{0}& \hspace{0.17em}& h_L& \dots & h_1& h_0\end{array}\right]}^{\mathrm{H}}.$$

The coefficient vector h = [h₀,…,h_L], 0 < L < K, corresponds to the maximum eigenvector of the matrix $\overline{\mathbf{R}}$ given by the average of the (K – L) submatrices ${\widehat{R}}_l\in {ℂ}^{\left(L+1\right)\times \left(L+1\right)}$ along the main diagonal of $\widehat{R}$ [¹⁵]:

(8) $$\overline{\mathbf{R}}=\frac{1}{K-L}{\displaystyle \sum }_{l=\left(L+1\right)}^K{\widehat{R}}_l.$$

Then, we can find h from the eigendecomposition of $\overline{\mathbf{R}}$:

(9) $$\overline{\mathbf{R}}={\displaystyle \sum }_{i=1}^{L+1}{\overline{\lambda }}_i{\overline{\mathbf{q}}}_i{\overline{\mathbf{q}}}_i^{\mathrm{H}},$$

where the eigenvalues ${\overline{\lambda }}_i$ are arranged in descending order of magnitude and ${\overline{\mathbf{q}}}_i$ are their corresponding eigenvectors. In other words, the maximum spatial eigenfilter is simply given by $\mathbf{h}={\overline{\mathbf{q}}}_1$ [¹⁵].

III. Problem Statement and Diagnosis

We had already stated in [¹⁸] that, despite the optimality of this spatial eigenfilter in maximizing the SNR at its output, its magnitude response may interchangeably attenuate signal sources widely-spaced from the others even in absence of noise. In such cases, the filter may cause Modified MODEX to take a wrong decision in the ML selection procedure, yielding a severe degradation in estimation performance. Although this unexpected behavior has been observed experimentally, it still lacks an analytical explanation.

To derive a mathematical reasoning about the selectivity of that filter, we must bring our discussion to the frequency domain. Then, the magnitude response H(e^jω) of the maximum spatial eigenfilter is given by:

(10) $$H\left(e^{j\omega }\right)\triangleq {\displaystyle \sum }_{k=0}^L{h}_k{e}^{-jk\omega }={\displaystyle \sum }_{k=1}^{L+1}{\overline{q}}_{k,1}{e}^{-j\left(k-1\right)\omega },$$

where ${\overline{q}}_{k,1}$ is the k-th element of ${\overline{\mathbf{q}}}_1$ and ω = π sin θ. To assess the effectiveness of H(e^jω), we first estimate the spatial power spectrum of the signals impinging on the array by taking the Fourier transform of the first column of $\widehat{R}$ as:

(11) $$\widehat{S}\left(\omega \right)\triangleq {\displaystyle \sum }_{m=0}^{K\hspace{0.17em}-\hspace{0.17em}1}\widehat{r}\left(m\right)e^{-jm\omega }={\displaystyle \sum }_{k=1}^K{\widehat{r}}_{k,\hspace{0.17em}1}{e}^{-j\left(k\hspace{0.17em}-\hspace{0.17em}1\right)\omega },$$

where ${\widehat{r}}_{k,\hspace{0.17em}1}$ is its k-th element and $\widehat{r}\left(m\right)$ is the estimated spatial correlation function at lag m.

According to [¹⁹]–[²¹], the difference between ${\overline{\lambda }}_1$ and ${\overline{\lambda }}_2$ becomes larger whenever signal sources are more closely spaced. In this case, ${\overline{\lambda }}_1\gg {\overline{\lambda }}_2$ so that ${\overline{\lambda }}_1$ holds almost all signal energy, as shown by the upper dashed line in Fig. 1 for {θ₁, θ₂} = {10°, 15°}. Thus, we can write Eq. (9) as:

(12) $$\begin{array}{ll}\overline{\mathbf{R}}\hfill & \approx {\overline{\lambda }}_1{\overline{\mathbf{q}}}_1{\overline{\mathbf{q}}}_1^{\mathrm{H}}\hfill \\ \hspace{0.17em}\hfill & \approx \left[{\overline{\lambda }}_1{\overline{\text{q}}}_{1,1}^{*}{\overline{\mathbf{q}}}_1,\dots ,{\overline{\lambda }}_1{\overline{\text{q}}}_{\left(L+1\right),1}^{*}{\overline{\mathbf{q}}}_1\right]\hfill \\ \hspace{0.17em}\hfill & \approx \left[{\overline{\mathbf{r}}}_1,\dots ,{\overline{\mathbf{r}}}_{\left(L+1\right)}\right],\hfill \end{array}$$

where ${\overline{\mathbf{r}}}_i$, 1 ≤ i ≤ (L + 1), corresponds to the i-th column of $\overline{\mathbf{R}}$.

Fig. 1
Percentage of signal energy in ${\overline{\lambda }}_1$ as a function of the separation between two signal sources with DOAs θ₁ = 10° and −40° ≤ θ₂ ≤ 60° for SNR ranging from −15 to 15 dB. The dashed lines regard {θ₁, θ₂} = {10°, 15°} and {θ₁, θ₂} = {10°, 45°}.

Therefore, $\mathbf{h}=\hspace{0.17em}{\overline{\mathbf{q}}}_1\approx {\overline{\mathbf{r}}}_1/{\overline{\lambda }}_1{\overline{q}}_{1,1}^{*}$ such that H(e^jω) corresponds to the spatial power spectrum calculated using only (L + 1) samples:

(13) $$H\left(e^{j\omega }\right)\approx \frac{1}{{\overline{\lambda }}_1{\overline{q}}_{1,1}^{*}}{\displaystyle \sum }_{k=1}^{L+1}{\overline{r}}_{k,\hspace{0.17em}1}{e}^{-j\left(k\hspace{0.17em}-\hspace{0.17em}1\right)\omega }.$$

As a consequence, the filter passband encompasses both DOAs of the signal sources, what explains why the maximum spatial eigenfilter performed so well for {θ₁, θ₂} = {10°, 15°} in [¹⁵]. Fig. 2a and 2b respectively depict the average spatial power spectrum of the impinging signals on 100 independent experiments and 100 magnitude responses of the maximum spatial eigenfilter. In both figures, vertical dashed lines indicate the actual DOAs. One can observe that all the magnitude responses of that filter resemble the average spatial power spectrum of the signal.

Fig. 2
a) Average of 100 signal spatial power spectra and b) superimposition of 100 normalized magnitude responses of the maximum spatial eigenfilter for {θ₁, θ₂} = {10°, 15°}, both in absence of noise

On the other hand, this is not true for widely-spaced signal sources. In this case, ${\overline{\lambda }}_1$ and ${\overline{\lambda }}_2$ get close to each other such that ${\overline{\lambda }}_1$ holds only half of signal energy, as indicated by the lower dashed line in Fig. 1 for {θ₁, θ₂} = {10°, 45°}. So, Eq. (12) no longer holds, since the first column of $\overline{\mathbf{R}}$ depends more heavily on ${\overline{\mathbf{q}}}_2,\dots ,{\overline{\mathbf{q}}}_{L+1}$:

(14) $${\overline{\mathbf{r}}}_1={\overline{\lambda }}_1{\overline{q}}_{1,1}^{*}{\overline{\mathbf{q}}}_1+{\displaystyle \sum }_{i=2}^{L+1}{\overline{\lambda }}_i{\overline{q}}_{i,1}^{*}{\overline{\mathbf{q}}}_i.$$

Then, the filter coefficient vector h becomes rather different from the spatial correlation vector and the magnitude response |H(e^jω)| of the maximum spatial eigenfilter may not sufficiently resemble S(ω). Consequently, |H(e^jω)| is not guaranteed to place multiple passbands around widely-spaced signal sources. Instead, the filter may interchangeably attenuate them even in absence of noise, as shown in Fig. 3a and 3b for {θ₁, θ₂} = {10°, 45°}.

Fig. 3
a) Average of 100 spatial power spectra and b) superimposition of 100 normalized magnitude responses of the maximum spatial eigenfilter for {θ₁, θ₂} = {10°, 45°}, both in absence of noise.

In both cases, Fig. 1 shows that the percentage of signal energy in ${\overline{\lambda }}_1$ gets smaller as the SNR decreases. So, the maximum eigenvector ${\overline{\mathbf{q}}}_1$ becomes even less suitable to approximate ${\overline{\mathbf{r}}}_1$ such that the maximum spatial eigenfilter may not be able to preserve all signal sources at low SNRs.

IV. Maximizing the output signal power with a signal subspace filter

The observations presented in the last section lead us to conclude that the maximum spatial eigenfilter may, in principle, sacrifice some of the signal sources to maximize the SNR at its output, causing Modified MODEX to make an incorrect decision in the ML procedure.

In [¹⁸], we introduced a Moving Average eigenvalue-based multiband spatial filter computed from the differential spectrum. That filter is able to place passbands on the DOAs of the signal sources at the expense of performing multiple eigendecompositions of $\widehat{R}$. However, to avoid a prohibitive growth in runtime, we had to reduce both the spectral resolution and filter order, what weakened its ability to attenuate noise outside the passbands. Then, for widely-spaced signal sources, even though that filter performed significantly better than the maximum spatial eigenfilter, it still performed slightly worse than the conventional Modified MODEX.

Following from Eq. (14), if a spatial filter whose coefficient vector $\tilde{h}$ was made dependent on ${\overline{\mathbf{q}}}_1,\dots ,{\overline{\mathbf{q}}}_{\left(L+1\right)}$ in the same way as ${\overline{\mathbf{r}}}_1$:

(15) $$\tilde{h}={\displaystyle \sum }_{i=1}^{L+1}{c}_i{\overline{\mathbf{q}}}_i,$$

then its magnitude response would be supposed to place multiple passbands around the signal sources regardless of whether they are closely or widely spaced. Additionally, this eigenvector-based multiband filter would be faster to calculate than that differential spectrum-based filter. So, in order to preserve all signal sources at the filter output and reduce the noise influence on the filter calculation, we propose to derive a finite impulse response (FIR) filter that maximizes the output SNR subject to its coefficient vector $\tilde{h}$ being dependent on the eigenvectors of only the signal subspace of $\overline{\mathbf{R}}$. According to [²⁶], the output SNR of a FIR filter is given by:

(16) $$\text{SNR}=\frac{{\tilde{h}}^{\mathrm{H}}\overline{R}\tilde{h}}{{\sigma }^2{\tilde{h}}^{\mathrm{H}}\tilde{h}}.$$

So, the optimization of the filter coefficient vector can be written as:

(17) $$\tilde{h}=\mathrm{arg\hspace{0.17em}}\underset{\tilde{h}}{\max }\left\{\frac{{\tilde{h}}^{\mathrm{H}}\overline{R}\tilde{h}}{{\sigma }^2{\tilde{h}}^{\mathrm{H}}\tilde{h}}\right\},\text{subject\hspace{0.17em}to\hspace{0.17em}}\tilde{h}={\displaystyle \sum }_{k=1}^{\overline{M}}{c}_k{\overline{\mathbf{q}}}_k.$$

Since σ² does not influence the maximization, the problem simply becomes:

(18) $$\tilde{h}=\mathrm{arg\hspace{0.17em}}\underset{\tilde{h}}{\max }\left\{\frac{{\tilde{h}}^{\mathrm{H}}\overline{R}\tilde{h}}{{\tilde{h}}^{\mathrm{H}}\tilde{h}}\right\},\text{subject\hspace{0.17em}to\hspace{0.17em}\hspace{0.17em}}\tilde{h}={\displaystyle \sum }_{k=1}^{\overline{M}}{c}_k{\overline{\mathbf{q}}}_k.$$

Then, if we replace $\overline{\mathbf{R}}$ with Eq. (9), the objective function can be written as:

(19) $$\begin{array}{ll}\tilde{h}\hfill & \begin{array}{c}=\mathrm{arg\hspace{0.17em}}\underset{\tilde{h}}{\max }\left\{\hspace{0.17em}\frac{{\tilde{h}}^{\mathrm{H}}\hspace{0.17em}\left({{\displaystyle \sum }}_{i=1}^{L+1}{\overline{\lambda }}_i{\overline{\mathbf{q}}}_i{\overline{\mathbf{q}}}_i^{\mathrm{H}}\right)\tilde{h}}{{\tilde{h}}^{\mathrm{H}}\tilde{h}}\right\}\end{array}\hfill \\ \hfill & =\mathrm{arg\hspace{0.17em}}\underset{\tilde{h}}{\max }\left\{\hspace{0.17em}\hspace{0.17em}\frac{{{\displaystyle \sum }}_{i=1}^{L+1}\left[{\tilde{h}}^{\mathrm{H}}\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i\right]{\left[{\tilde{h}}^{\mathrm{H}}\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i\right]}^{\mathrm{H}}}{{\tilde{h}}^{\mathrm{H}}\tilde{h}}\right\}\hfill \\ \hfill & =\mathrm{arg\hspace{0.17em}}\underset{\tilde{h}}{\max }\left\{\hspace{0.17em}\frac{{{\displaystyle \sum }}_{i=1}^{L+1}{\left|{\tilde{h}}^{\mathrm{H}}\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i\right|}^2}{{||\tilde{h}||}^2}\right\}\hfill \\ \hfill & =\mathrm{arg\hspace{0.17em}}\underset{\tilde{h}}{\max }\left\{{\displaystyle \sum }_{i=1}^{L+1}{\left|\frac{<\tilde{h},\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i>}{||\tilde{h}||}\right|}^2\right\}.\hfill \end{array}$$

where 〈·,·〉 stands for the inner product. But the coefficient vector $\tilde{h}$ that maximizes Eq. (19) is the same that maximizes each of its addends, so:

(20) $$\begin{array}{cc}\tilde{h}& \begin{array}{c}=\mathrm{arg\hspace{0.17em}}\underset{\tilde{h}}{\max }\left\{\frac{<\tilde{h},\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i>}{||\tilde{h}||}\right\}\end{array}\end{array}$$

Since ${\overline{\lambda }}_i$ is not a function of $\tilde{h}$, we can do:

(21) $$\begin{array}{cc}\tilde{h}& \begin{array}{c}=\mathrm{arg\hspace{0.17em}}\underset{\tilde{h}}{\max }\left\{\frac{<\tilde{h},\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i>}{||\tilde{h}||\left|\sqrt{{\overline{\lambda }}_i}\right|}\right\}\end{array}\end{array}$$

such that the expression inside curly brackets is the cosine similarity between $\tilde{h}$ and $\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i$. Now, by imposing the constraint, Eq. (21) becomes:

(22) $${\text{c}}_k=\arg \underset{c_k}{\max }\left\{\frac{\langle {\displaystyle {\sum }_{k=1}^{\overline{M}}{c}_k{\overline{\text{q}}}_k},\sqrt{{\overline{\lambda }}_i}{\overline{\text{q}}}_i\rangle }{\Vert {\displaystyle {\sum }_{k=1}^{\overline{M}}{c}_k{\overline{\text{q}}}_k}\Vert \left|\sqrt{{\overline{\lambda }}_i}\right|}\right\}$$

Considering the cosine similarity is maximum for identical vectors, we conclude that:

(23) $${\displaystyle \sum }_{k=1}^{\overline{M}}{c}_k{\overline{\mathbf{q}}}_k=\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i$$

Pre-multiplying Eq. (23) with ${\overline{\mathbf{q}}}_i^{\mathrm{H}}$, we get:

(24) $${\displaystyle \sum }_{k=1}^{\overline{M}}{c}_k{\overline{\mathbf{q}}}_i^{\mathrm{H}}{\overline{\mathbf{q}}}_k=\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i^{\mathrm{H}}{\overline{\mathbf{q}}}_i$$

since ${\overline{\mathbf{q}}}_i^{\mathrm{H}}{\overline{\mathbf{q}}}_k=\delta \left(i-k\right)$, we find $c_i=\sqrt{{\overline{\lambda }}_i},\hspace{0.17em}i=1,\dots ,\overline{M}$, such that the solution to Eq. (17) is given by:

(25) $$\tilde{h}={\displaystyle \sum }_{i=1}^{\overline{M}}\sqrt{{\overline{\lambda }}_i}{\overline{\mathbf{q}}}_i.$$

Actually, this solution corresponds to a parallel combination of multiple filters ${\tilde{h}}_i=\hspace{0.17em}{\overline{\mathbf{q}}}_i,\hspace{0.17em}i=1,\hspace{0.17em}\dots ,\hspace{0.17em}\overline{M}$, each one with magnitude response:

(26) $${\tilde{H}}_i\left(e^{j\omega }\right)={\displaystyle \sum }_{k=1}^{L+1}{\overline{q}}_{k,i}{e}^{-j\left(k\hspace{0.17em}-\hspace{0.17em}1\right)\omega },i=1,\dots ,\overline{M}.$$

As shown in Fig. 4a, the magnitude responses $\left|{\tilde{H}}_i\left(e^{j\omega }\right)\right|$ are complementary and their combination place passbands around the actual DOAs. However, the corresponding phase responses $\angle {\tilde{H}}_i\left(e^{j\omega }\right)$ do not match, as can be seen in Fig. 4b. So, before combining the filters, we can equalize the phase response of each one of them by simply convolving its impulse response with its corresponding time-reversed, time-shifted, conjugated version, in order to avoid magnitude and phase distortion. Then, this allows us to derive a filter of order 2L whose coefficient vector $\widehat{h}$ is given by the following weighted sum:

(27) $$\begin{array}{c}\widehat{h}={\displaystyle \sum }_{i=1}^{\overline{M}}{\overline{\text{λ}}}_i\text{conv}\left({\tilde{h}}_i,\mathbf{P}{\tilde{h}}_i^{\text{*}}\right)={\displaystyle \sum }_{i=1}^{\overline{M}}{\overline{\text{λ}}}_i\text{conv}\left({\overline{\mathbf{q}}}_i,\mathbf{P}{\overline{\mathbf{q}}}_i^{\text{*}}\right),\end{array}$$

Fig. 4
a) Magnitude responses and b) phase responses for ${\tilde{H}}_1\left(e^{j\omega }\right)$ and ${\tilde{H}}_2\left(e^{j\omega }\right)$ in a single experiment in absence of noise considering {θ₁, θ₂} = {10°, 45°}.

where conv(·) stands for the convolution operator, and P is the (L + 1) × (L + 1) co-identity matrix [²⁵].

In frequency domain, the frequency response of this linear phase filter corresponds to:

(28) $$\widehat{H}\left(e^{j\omega }\right)={\displaystyle \sum }_{i=1}^{\overline{M}}{\overline{\text{λ}}}_i{\left|{\tilde{H}}_i\left(e^{j\omega }\right)\right|}^2{e}^{-j\omega L},$$

In order to evaluate the effectiveness of this filter, we applied it to the cases {θ₁, θ₂} = {10°, 15°} and {θ₁, θ₂} = {10°, 45°}, as shown respectively in Fig. 5a and 5b. The magnitude response of the proposed filter was able to place multiple passbands around the DOAs of the signal sources, preserving the power of both signals regardless of their spacing.

Fig. 5
Superimposition of 100 normalized magnitude responses of the proposed filter in absence of noise for a) {θ₁, θ₂} = {10°, 15°} and b) {θ₁, θ₂} = {10°, 45°}.

However, using the entire signal subspace of $\overline{\mathbf{R}}$ increases the probability of subspace swap [²⁶], [²⁷]. So, spurious peaks are more likely to arise in the filter magnitude response for closely-spaced signal sources at low SNR values. This may compromise the ML procedure and cause performance degradation under severe noise conditions [²⁸].

V. Results

In this section, we compare the estimation performance of Modified MODEX using our present proposition during the ML procedure to those using the maximum spatial eigenfilter [¹⁵], and the differential spectrum-based filter [¹⁸].

To do so, we performed T = 1000 Monte-Carlo experiments for each SNR ranging from −15 to 15 dB in steps of 1.25 dB, considering a ULA composed of K = 10 sensors that take N = 100 snapshots of M = 2 narrowband signal sources. We adopted order L = 7 as in [¹⁵] for all filters evaluated in this work. Additionally, the differential spectrum was calculated using 32 samples. Since the filters were evaluated according to the signal source separation, we considered one signal source fixed at θ₁ = 10°, while the other varied in the interval −40° ≤ θ₂ ≤ 60° in steps of 2.5°. Additionally, the case where both signal sources have θ₁ = θ₂ = 10° was unconsidered.

The estimation performance was evaluated using Root Mean Square Error (RMSE) for each SNR value. RMSE is compared to Cramér-Rao Lower Bound (CRLB) [²⁹] and the threshold SNR is discussed.

Figs. 6 to 8 respectively show gray-shaded RMSE surfaces for Modified MODEX using the maximum spatial eigenfilter, the differential spectrum-based filter, and the filter proposed in this work. In those figures, the dashed lines represent the reference scenarios {θ₁, θ₂} = {10°, 15°} and {θ₁, θ₂} = {10°, 45°}. The dashed surface corresponds to the CRLB. As a reference, the threshold SNR of the conventional Modified MODEX is −6.25 dB for {θ₁, θ₂ dB for {θ₁, θ₂} = {10°, 15°} and −10 dB for {θ₁, θ₂} = {10°, 45°}.

Fig. 6
RMSE surface of Modified MODEX using the maximum spatial eigenfilter, and the corresponding CRLB. The dashed lines regard {θ₁,θ₂} = {10°, 15°} and {θ₁, θ₂} = {10°, 45°}.

Fig. 7
RMSE surface of Modified MODEX using the filter calculated using the differential spectrum-based filter, and the corresponding CRLB. The dashed lines regard {θ₁, θ₂} = {10°, 15°} and {θ₁, θ₂} = {10°, 45°}.

Fig. 8
RMSE surface of Modified MODEX using the proposed filter, and the corresponding CRLB. The dashed lines regard {θ₁, θ₂} = {10°, 15°} and {θ₁, θ₂} = {10°, 45°}.

The maximum spatial eigenfilter significantly improved the estimation performance for closely-spaced signal sources such that its RMSE surface was below the CRLB, as presented in Fig. 6. The threshold SNR for {θ₁, θ₂} = {10°, 15°} was −12.50 dB. However, the maximum spatial eigenfilter caused a severe degradation on estimation performance for widely-spaced signal sources. For {θ₁, θ₂} = {10°, 45°}, the threshold SNR was 3.75 dB, far above that of the conventional Modified MODEX.

The differential spectrum-based filter performed a little worse than the maximum spatial eigenfilter for closely-spaced signal sources, as shown in Fig. 7. For instance, its threshold SNR was −11.25 dB for {θ₁, θ₂} = {10°, 15°}. On the other hand, it significantly improved the estimation performance for widely-spaced signal sources compared to the maximum spatial eigenfilter. For {θ₁, θ₂} = {10°, 45°}, its threshold SNR was −8.75 dB which is much below that of the maximum spatial eigenfilter but still slightly above that of the conventional Modified MODEX.

Finally, Fig. 8 shows that our proposition performed similarly to the other two filters for closely-spaced signal sources. However, it is worth noting that the RMSE surface attained the CRLB and the threshold SNR was mostly kept at −10 dB regardless of the signal source separation. For instance, it achieved −10 dB for both {θ₁, θ₂} = {10°, 15°} and {θ₁, θ₂} = {10°, 45°}. Therefore, the proposed spatial filter in no way hampered the estimation performance of the conventional Modified MODEX and significantly improved it for closely-spaced signal sources, keeping the threshold SNR at −10 dB for almost any signal source separation.

As a total of one million experiments of Modified MODEX were carried out to calculate the RMSE surfaces for each filter, we took this opportunity to record their respective computation times. Regarding the computational effort, filtering techniques clearly increase the total computation time, since additional calculations are needed. Then, we measured the runtime of each experiment for Modified MODEX taking into account both the DOA estimation time and filter calculation time when applicable. As a reference, the mean runtime for the conventional Modified MODEX was 3.19 milliseconds (ms) per experiment. Modified MODEX using maximum spatial eigenfiltering took 3.49 ms on average. With the differential spectrum-based filter, this time was 4.11 ms, and using our present proposition, the mean runtime was 3.62 ms. In relative terms, Modified MODEX using our present proposition was 13% slower than the conventional version, 4% slower than using the maximum spatial eigenfilter but 12% faster than using the differential spectrum-based filter.

We also must note that our proposition makes use of the convolution operation – Eq. (27) – and, with L = 7, its impulse response vector has 15 coefficients instead of only 8 for the other filters. This causes $\tilde{R}$ to assume larger dimensions and, consequently, the amount of time spent to calculate Eq. (6) in the ML procedure is larger. Even so, our proposition performed faster than the differential spectrum-based filter, while delivering better results.

VI. Conclusions and remarks

Although the maximum spatial eigenfilter allowed a significant improvement in the DOA estimation of closely-spaced signal sources, we observed that its magnitude response may interchangeably attenuate widely-spaced signal sources, leading to a severe performance degradation and increasing the threshold SNR. We pointed that out firstly in [¹⁸], but without an accompanying analytical explanation. At that time, we proposed the differential spectrum-based filter that, despite largely overcoming the limitations of the maximum spatial eigenfilter, it still caused a slightly increase in the threshold SNR for widely-spaced signal sources.

From Karhunen-Loève expansion, we mathematically explained this unexpected behavior in terms of the relation between the energy of the principal eigenvalues of the spatial correlation matrix and the separation between the signal sources, as pointed out in [²¹]. We proved that, for closely-spaced signal sources, the maximum eigenvalue concentrates almost all signal energy and the coefficients of the maximum spatial eigenfilter approximate the signal spatial correlation vector. In this case, the filter magnitude response places passbands around the DOAs of the signal sources. On the other hand, for widely-spaced signal sources, half of signal energy spreads through the signal subspace such that the maximum spatial eigenfilter no longer approximates the spatial power spectrum and can even suppress any of the signal sources.

As a corollary, to circumvent that problem, we derived a new multiband filter that maximizes signal power at its output by combining the eigenvectors of the entire signal subspace. To avoid phase distortion, we performed a phase equalization procedure and the resulting filter was able to place multiple passbands around the DOAs of all signal sources, regardless of their separation.

Experimental results showed that using our proposition during the ML procedure in no way hampered the estimation performance of the conventional Modified MODEX and significantly improved it for closely-spaced signal sources. In addition, it kept the threshold SNR almost independent of signal source separation at the expense of a small increase in runtime.

As future work, we plan extending the scope of spatial filtering to other important DOA estimation methods beyond MODEX-based estimators. In addition, we intend to evaluate the distortion effects caused by the non-flat gain of the filter passbands on the DOA estimation in two or three dimensions for wideband systems.

ACKNOWLEDGEMENT

This work was supported by the Instituto Federal de Goiás under grant number 23380.000425/2014-09.

REFERENCES

• ^[1]
  R. Zhang, W. Xia, F. Yan, L. Shen, “A Single-Site Positioning Method Based on TOA and DOA Estimation Using Virtual Stations in NLOS Environment”, China Communications, vol. 16, no. 2, pp. 146–159, Feb. 2019.
• ^[2]
  Charles Schroeder, “The Road to Realizing 5G Technologies”, in Telecom Leadership Forum 2018 (TLF 2018), 2018.
• ^[3]
  S. Chen, S. Sun, G. Xu, X. Su, Y. Cai, “Beam-space Multiplexing: Practice, Theory, and Trends-From 4G TD-LTE, 5G, to 6G and Beyond”, IEEE Wireless Communications, vol. 27, no. 2, pp. 162–172, Apr. 2020.
• ^[4]
  W. G. Tang, H. Jiang, Q. Zhang, “ADMM For Gridless DOD and DOA Estimation In Bistatic MIMO Radar Based On Decoupled Atomic Norm Minimization With One Snapshot”, in 2019 IEEE Global Conference on Signal and Information Processing (GlobalSIP), Ottawa, Canada, 2019.
• ^[5]
  Y. Chen, F. Wang, J. Wan, K. Xu, “DOA Estimation for Long-distance Underwater Acoustic Sources based on Signal Self-nulling”, in 2017 IEEE 2nd Advanced Information Technology, Electronic and Automation Control Conference (IAEAC), Chongquing, China, Mar. 2017.
• ^[6]
  R. O. Schmidt, “Multiple emitter location and signal parameter estimation”, IEEE Transactions on Antennas and Propagation, vol. 34, no. 3, pp. 276–280, 1986.
• ^[7]
  R. H. Roy and T. Kailath, “ESPRIT-Estimation of signal parameters via rotational invariance techniques”, IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 37, no. 7, pp. 984–995, 1989.
• ^[8]
  P. Stoica, K. C. Sharman, “Novel eigenanalysis method for direction estimation”, IEE Proceedings F (Radar and Signal Processing), vol. 137, no. 1, pp. 19–26, 1990.
• ^[9]
  C. Qian, L. Huang, M. Cao, J. Xie and H. C. So, “PUMA: An Improved Realization of MODE for DOA Estimation”, IEEE Transactions on Aerospace and Electronic Systems, vol. 53, no. 5, pp. 2128-2139, Oct. 2017.
• ^[10]
  C. Qian, L. Huang, H. C. So, N. D. Sidiropoulos, and J. Xie, “Unitary PUMA algorithm for estimating the frequency of a complex sinusoid”, IEEE Transactions on Signal Processing, vol. 63, no. 20, pp. 5358–5368, 2015.
• ^[11]
  C. Qian, L. Huang, N. Sidiropoulos, and H. C. So, “Enhanced PUMA for direction-of-arrival estimation and its performance analysis”, IEEE Transactions on Signal Processing, vol. 64, no. 16, pp. 4127–4137, Aug. 2016.
• ^[12]
  D. Zachariah, P. Stoica and M. Jansson, “Comments on ‘Enhanced PUMA for Direction-of-Arrival Estimation and Its Performance Analysis’”, in IEEE Transactions on Signal Processing, vol. 65, no. 22, pp. 6113-6114, Nov.15, 2017.
• ^[13]
  A. B. Gershman, P. Stoica, “New MODE-based techniques for direction finding with an improved threshold performance”, Signal Processing, vol. 76, pp. 221–235, 1999.
• ^[14]
  A. Lopes, I. S. Bonatti, P. L. D. Peres, C. A. Alves, “Improving the MODEX algorithm for direction estimation”, Signal Processing, vol. 83, no. 9, pp. 2047–2051, Sep. 2003.
• ^[15]
  R. Krummenauer, M. Cazarotto, A. Lopes, P. Larzabal, P. Forster, “Improving the threshold performance of maximum likelihood estimation of direction of arrival”, Signal Processing, vol. 90, no. 11, pp. 1582–1590, Nov. 2010.
• ^[16]
  P. Forster, G. Vezzosi, “Application of spheroidal sequences to array processing”, in IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP ‘87), Dallas, USA, 1987, pp. 2268–2271.
• ^[17]
  S. Haykin, Adaptive Filter Theory, 4th ed., Englewood Cliffs, NJ, USA: Prentice-Hall, 2001.
• ^[18]
  R. P. Lemos, H. V. L. Silva, E. L. Flôres, J. A. Kunzler, D. F. B. Beltrán, “Spatial filtering based on differential spectrum for improving ML DOA estimation performance”, IEEE Signal Processing Letters, vol. 23, no. 12, pp. 1811–1815, Aug. 2016.
• ^[19]
  Y. R. Ferreira, R. P. Lemos, “A new DOA estimation algorithm based on angle search through the difference between the principal singular values”, in Proceedings of the International Microwave and Optoelectronics Conference, Brasília, Brazil, 2005, pp. 283–286.
• ^[20]
  Y. R. Ferreira, R. P. Lemos, “A new DOA estimation algorithm based on differential spectrum”, in Proceedings of the Eighth International Symposium on Signal Processing and Its Applications, Sydney, Australia, 2005, pp. 303–307.
• ^[21]
  H. V. L. Silva, R. P. Lemos, Y. R. Ferreira, L. G. R. Guedes, “A branch-and-bound inspired technique to improve the computational efficiency of DOA estimation”, Signal Processing, vol. 93, no. 4, pp. 947–956, Apr. 2013.
• ^[22]
  S. M. Kay, Modern Spectral Estimation - Theory and Application, Englewood Cliffs, NJ, USA: Prentice Hall, 1988.
• ^[23]
  P. Stoica, A. Nehorai, “Performance study of conditional and unconditional direction-of-arrival estimation”, IEEE Transactions on Acoustic, Speech and Signal Processing, vol. 38, no. 10, pp. 1783–1795, Oct. 1990.
• ^[24]
  F. Li, R. J. Vaccaro, “Unified analysis for DOA estimation algorithms in array signal processing”, Signal Processing, vol. 25, pp. 147–169, 1991.
• ^[25]
  M. G. Bellanger, Adaptive Digital Filters, 2nd ed., New York, NY, USA: Marcel Dekker, 2001.
• ^[26]
  J. K. Thomas, L. L. Scharf, D. W. Tufts, “The probability of a subspace swap in the SVD”, IEEE Transactions on Signal Processing, vol. 43, no. 3, pp. 730–736, Mar. 1995.
• ^[27]
  M. Hawkes, A. Nehorai, P. Stoica, “Performance breakdown of subspace-based methods: Prediction and cure”, in Proceedings of 2001 IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP ‘01), Salt Lake City, USA, 2001, pp. 4005–4008.
• ^[28]
  M. Shaghaghi, S. A. Vorobyov, “Subspace leakage analysis and improved DOA estimation with small sample size”, IEEE Transactions on Signal Processing, vol. 63, no. 12, pp. 3251–3265, Jun. 2015.
• ^[29]
  H. L. V. Trees, Optimum Array Processing. Part IV of Detection, Estimation and Modulation Theory, New York, NY, USA: John Wiley and Sons, 2002.

Publication Dates

• Publication in this collection
  03 Mar 2021
• Date of issue
  Mar 2021

History

• Received
  12 May 2020
• Reviewed
  14 May 2020
• Accepted
  03 Dec 2020