Wideband Subspace Estimation Through Projection Matrix Approximation

TL;DR

This paper introduces a polynomial approximation approach for wideband subspace estimation, reducing parameters and improving DOA estimation accuracy by modeling the signal subspace's projection matrix across frequencies.

Contribution

It proposes a novel polynomial fitting method for the projection matrix to enhance wideband subspace estimation and DOA accuracy, with theoretical bounds and numerical validation.

Findings

Reduces the number of parameters needed for wideband subspace representation.

Improves the accuracy of wideband DOA estimators like IC-MUSIC and MTOPS.

Provides asymptotic bounds for estimation bias and RMS error.

Abstract

In this paper, we present a wideband subspace estimation method that characterizes the signal subspace through its orthogonal projection matrix at each frequency. Fundamentally, the method models this projection matrix as a function of frequency that can be approximated by a polynomial. It provides two improvements: a reduction in the number of parameters required to represent the signal subspace along a given frequency band and a quality improvement in wideband direction-of-arrival (DOA) estimators such as Incoherent Multiple Signal Classification (IC-MUSIC) and Modified Test of Orthogonality of Projected Subspaces (MTOPS). In rough terms, the method fits a polynomial to a set of projection matrix estimates, obtained at a set of frequencies, and then uses the polynomial as a representation of the signal subspace. The paper includes the derivation of asymptotic bounds for the bias and root-mean-square (RMS) error of the projection matrix estimate and a numerical assessment of the method and its combination with the previous two DOA estimators.

Tables

Table 1. TABLE I: Chebyshev interpolation error using measure in ( 63 ).

	1	2	3	4
0	0.1486	0.2624	0.3313	0.3766
1	0.06228	0.09993	0.1106	0.114
2	0.01739	0.02558	0.03154	0.03407
3	0.003732	0.007473	0.01038	0.01325
4	0.0006299	0.00211	0.003046	0.00411
5	0.00008891	0.0007461	0.001155	0.00163

Table 2. TABLE II: Chebyshev interpolation error using measure in ( 64 ).

	1	2	3	4
0	0.5051	0.6319	0.6578	0.6794
1	0.1552	0.1986	0.2247	0.2454
2	0.03203	0.0652	0.07654	0.08669
3	0.005908	0.02142	0.02551	0.03059
4	0.0008885	0.006993	0.00854	0.01067
5	0.0001158	0.00227	0.002881	0.003813

Equations

$$\mbox{\boldmath$P$\unboldmath}(f)=\mbox{\boldmath$A$\unboldmath}(f)(\mbox{\boldmath$A$\unboldmath}(f)^{H}\mbox{\boldmath$A$\unboldmath}(f))^{-1}\mbox{\boldmath$A$\unboldmath}(f)^{H}.$$

$$\mbox{\boldmath$P$\unboldmath}(f)\approx\sum_{q=0}^{Q}\mbox{\boldmath$G$\unboldmath}_{q}f^{q}$$

$$\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})\equiv\widehat{}\mbox{\boldmath$Q$\unboldmath}_{K}(f_{r})\widehat{}\mbox{\boldmath$Q$\unboldmath}_{K}(f_{r})^{H}$$

$$s_k(t)\equiv e^{-j2\pi f_{o}t}{s}_{bp,k}(t).$$

$$(\stackrel{~}{}{a}_{m,k}\ast s_k)(t).$$

$$\stackrel{~}{}{a}_{m,k}(t)=e^{-j2\pi f_o{\tau }_{m,k}}\delta (t-{\tau }_{m,k}).$$

$$x_m(t)=\sum _{k=1}^K(\tilde{}{a}_{m,k}\ast s_k)(t)+{\eta }_m(t),$$

$$|a_{m,k}^{\prime }(f)|\le \pi {\tau }_{max}.$$

$$\begin{gathered} (h(\cdot ;f_r)\ast x_m)(t) \\ =\sum _{k=1}^K(h(\cdot ;f_r)\ast \tilde{}{a}_{m,k}\ast s_k)(t)+(h(\cdot ;f_r)\ast {\eta }_m)(t) \end{gathered}$$

$$(h(\cdot ;f_r)\ast \tilde{}{a}_{m,k})(t)\approx a_{m,k}(f_r)h(t;f_r),$$

$$\begin{gathered} (h(\cdot ;f_r)\ast x_m)(t)= \\ \sum _{k=1}^K{a}_{m,k}(f_r)(h(\cdot ;f_r)\ast s_k)(t)+(h(\cdot ;f_r)\ast {\eta }_m)(t). \end{gathered}$$

$$\displaystyle[\mbox{\boldmath$A$\unboldmath}(f)]_{m,k}$$

$$\displaystyle[\mbox{\boldmath$x$\unboldmath}(t;f_{r})]_{m}$$

$$\displaystyle[\mbox{\boldmath$s$\unboldmath}(t;f_{r})]_{k}$$

$$\displaystyle[\mbox{\boldmath$\epsilon$\unboldmath}(t;f_{r})]_{m}$$

$$\mbox{\boldmath$x$\unboldmath}(t;f_{r})=\mbox{\boldmath$A$\unboldmath}(f_{r})\mbox{\boldmath$s$\unboldmath}(t;f_{r})+\mbox{\boldmath$\epsilon$\unboldmath}(t;f_{r}),r=1,\,2,\ldots,\,R.$$

$$\mbox{\boldmath$R$\unboldmath}(f_{r})\equiv\mbox{\boldmath$A$\unboldmath}(f_{r})\mbox{\boldmath$R$\unboldmath}_{s}(f_{r})\mbox{\boldmath$A$\unboldmath}(f_{r})^{H}+\sigma^{2}\mbox{\boldmath$I$\unboldmath}_{M},$$

$$\widehat{}\mbox{\boldmath$R$\unboldmath}(f_{r})\equiv\frac{1}{N}\sum_{n=0}^{N-1}\mbox{\boldmath$x$\unboldmath}(nT;f_{r})\mbox{\boldmath$x$\unboldmath}(nT;f_{r})^{H}.$$

$$\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})\equiv\widehat{}\mbox{\boldmath$Q$\unboldmath}_{K}(f_{r})\widehat{}\mbox{\boldmath$Q$\unboldmath}_{K}(f_{r})^{H}.$$

$$[\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})]_{m,m^{\prime}}\approx\sum_{q=0}^{Q}[\mbox{\boldmath$G$\unboldmath}_{q}]_{m,m^{\prime}}f_{r}^{q},\;\;r=1,\,2\ldots,\,R,$$

$$\displaystyle\mathcal{E}\{\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})\}$$

$$\displaystyle\text{Var}\{\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})\}$$

$$\sum_{r=1}^{R}w_{r}^{2}\Big{|}[\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})]_{m,m^{\prime}}-\sum_{q=0}^{Q}[\widehat{}\mbox{\boldmath$G$\unboldmath}_{q}]_{m,m^{\prime}}f_{r}^{q}\Big{|}^{2},$$

$$[\widehat{}\mbox{\boldmath$G$\unboldmath}_{q}]_{m,m^{\prime}}\equiv\sum_{r=1}^{R}[\mbox{\boldmath$B$\unboldmath}]_{q,r}[\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})]_{m,m^{\prime}},$$

$$B$$

$$\displaystyle[\mbox{\boldmath$\Lambda$\unboldmath}]_{r,r^{\prime}}$$

$$\widehat{}\mbox{\boldmath$G$\unboldmath}_{q}\equiv\sum_{r=1}^{R}[\mbox{\boldmath$B$\unboldmath}]_{q,r}\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})$$

$$\widehat{}\mbox{\boldmath$P$\unboldmath}_{1}(f)\equiv\sum_{q=0}^{Q}\widehat{}\mbox{\boldmath$G$\unboldmath}_{q}f^{q}.$$

$$\widehat{}\mbox{\boldmath$P$\unboldmath}_{1}(f)\equiv\sum_{r=1}^{R}{d}_{r}(f)\widehat{}\mbox{\boldmath$P$\unboldmath}_{0}(f_{r})$$

$${d}_{r}(f)=\sum_{q=0}^{Q}[\mbox{\boldmath$B$\unboldmath}]_{q,r}f^{q}.$$

$$\mbox{\boldmath$P$\unboldmath}(f)=\sum_{r=1}^{R}{d}_{r}(f)\mbox{\boldmath$P$\unboldmath}(f_{r}).$$

$$\displaystyle\mathcal{E}\{\widehat{}\mbox{\boldmath$P$\unboldmath}_{1}(f)\}$$

$$\displaystyle\text{Var}\{\widehat{}\mbox{\boldmath$P$\unboldmath}_{1}(f)\}$$

Full text

I Introduction

In array processing, the estimation of the subspace spanned by several sources is a fundamental step in DOA estimation [1, Ch. 9]. This estimation relies on the so-called narrowband condition, i.e, the array response must be constant in the spectral band covered by the sources. In practice, however, this condition is often unrealistic in applications involving high data rates or acoustic or seismic signals, in which the array response varies with frequency significantly. In the literature on DOA estimation, these cases are classified as “wideband” and addressed by dividing the signals’ band into bins in which the array response is approximately constant. The problem is then the way the data from all bins should be combined in order to produce a single set of DOA estimates. In the literature, there are fundamentally two approaches for this combination. The first is the coherent approach in which the data from all bins are linearly combined [2, 3, 4, 5], and the second is the incoherent approach in which the combination is performed by other means, such as averaging the DOA narrowband estimates from each bin or adding up the bin pseudo-spectrum functions (as is done in IC-MUSIC), [6]. Additionally, there exist other ways to process the bin data derived from general principles such as Maximum Likelihood (ML) [7, 8, 9, 10, 11], polynomial matrix decompositions [12, 13, 14, 15], or group sparsity [16, 17, 18].

There is a relevant feature in this estimation problem that seems to be overlooked in the literature, which is the fact that the array response varies smoothly along the band covered by the impinging signals, and the same happens with the corresponding signal subspace. This smooth variation is a consequence of two basic facts. First, for any sensor array, there is an upper bound for the distance among sensors that can be measured as a delay $\tau_{max}$ at the propagation speed. And second, if the array is viewed as a vector Linear Time-Invariant system, then its spectrum is a band-limited function of bandwidth at most $\tau_{max}$. Thus, the array response variation with frequency is smooth (a band-limited function) and we may even bound any of its derivatives using standard results; (see Bernstein’s inequality [19, Th. 6.7]). Additionally, if we consider a number of impinging waves from different DOAs, then the subspace spanned by the array response spectrum varies smoothly with frequency. To be more precise, the orthogonal projection matrix of the subspace varies smoothly with frequency and can be approximated in any given frequency band using a polynomial. This is a consequence of the fact that this projection matrix is analytic [20, Ch. 2, Th. 1.10] and Weiertrass theorem [21, Th. 2.4.1].

The purpose of this paper is to present a method for exploiting this last subspace smoothness in order to improve wideband subspace estimation. The method is based on a polynomial approximation and its input is the set of projection matrices that is usually computed in the binning approach. Fundamentally, the method consists of fitting a polynomial to these last matrices using weighted least squares, and then evaluating the polynomial at any frequency for obtaining an approximate projection matrix.

The next section is a extended introduction in which we justify the smooth variation of the signal subspace, outline the proposed method, and describe the organization of the paper.

I-A Notation and main symbols

We use the following notation and basic concepts:

• •

We write vectors in lower case ($a$, $x$) and matrices in upper case, ($A$, $X$).

• •

\mbox{\boldmathI\unboldmath}_{M} is an identity matrix of size ${M\times M}$ and \mbox{\boldmath0\unboldmath}_{M} the ${M\times M}$ zero matrix.

• •

[\mbox{\boldmatha\unboldmath}]_{m} and [\mbox{\boldmathA\unboldmath}]_{m,k} respectively denote the $m$th and $(m,k)$ components of $a$ and $A$. Also, [\mbox{\boldmathA\unboldmath}]_{\cdot,k} denotes the $k$th column of $A$ and [\mbox{\boldmathA\unboldmath}]_{\cdot,m:k} denotes the matrix form by its columns $m$ to $k$.

• •

\mbox{\boldmathA\unboldmath}^{H} is the Hermitian of $A$ and \mbox{\boldmathA\unboldmath}^{\dagger} its pseudo-inverse.

• •

The operator ’$\equiv$’ introduces new symbols.

• •

’$*$’ denotes convolution: $({a}*{b})(t)$ is the convolution of ${a}(t)$ and ${b}(t)$.

• •

’$\mathcal{E}\{\cdot\}$’ and ’$\text{Var}\{\cdot\}$’ denote the expectation and variance operators.

• •

${\delta}(t)$ and ${\delta}_{n}$ respectively denote Dirac and Kronecker delta.

• •

In the paper, a given ${M\times M}$ matrix $P$ is said to be a projection matrix if \mbox{\boldmathP\unboldmath}=\mbox{\boldmathP\unboldmath}^{H} and \mbox{\boldmathP\unboldmath}^{2}=\mbox{\boldmathP\unboldmath}.

• •

\|\mbox{\boldmathB\unboldmath}\|_{F} and \|\mbox{\boldmathB\unboldmath}\|_{2} denote the Frobenius and 2-norm of a given matrix $B$ respectively.

• •

${O}(\cdot)$ is the big-O notation. A function ${g}(N)$ is ${O}(1/N)$ if there are positive numbers $N_{o}$ and $A$ such that $|{g}(N)|\leq A/N$ if $N>N_{o}$.

• •

${o}(\cdot)$ is the little-o notation. A function ${g}(N)$ is ${o}(1/N)$ if for any $\epsilon>0$ there is positive $N_{o}$ such that $|{g}(N)|\leq\epsilon/N$ if $N>N_{o}$.

• •

In some contexts, ’$\cong$’ marks an equality that holds if an ${o}(1/N)$ term is added.

The main symbols in the paper are the following matrix functions:

• •

\mbox{\boldmathR\unboldmath}(f_{r}): expected covariance matrix at the $r$th bin, (14).

• •

\widehat{}\mbox{\boldmathR\unboldmath}(f_{r}): sample covariance matrix at $r$th bin, (15).

• •

\mbox{\boldmathP\unboldmath}(f): expected signal projection matrix function, (1).

• •

\widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}): initial signal projection matrix estimate at $r$th frequency bin, (16).

• •

\widehat{}\mbox{\boldmathP\unboldmath}_{1}(f): polynomial estimate of signal projection matrix, (25).

• •

\widehat{}\mbox{\boldmathP\unboldmath}_{2}(f): correction of \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) to the closest projection matrix, (30).

II Characterization of the signal subspace for wideband subspace estimation, method’s outline, and paper organization

Given a sensor array, the estimation of the subspace spanned by several impinging waves is, in principle, a well-posed problem only in narrow spectral bands, given that the array response varies with frequency. This implies that any extension to wide bands of subspace estimation must take into account the response variation in some way. A relevant feature of this variation is that it is smooth along the band covered by the incoming signals, and the same happens with the subspace spanned by several DOAs. Let us justify this assertion by performing the following analysis, valid for a generic sensor array.

Consider an array of $M$ sensors and $K$ impinging waves. Also, place the origin of coordinates at the center of the smallest circle or sphere containing all sensors, and let $\tau_{max}$ denote its diameter but measured as a delay at the propagation speed. The response of any of the sensors to any of the DOAs consists of the convolution with a response of the form $a{\delta}(t-\tau)$ for some factor $a$ and delay $\tau$ following ${|\tau|\leq\tau_{max}/2}$. Thus, the spectrum of this response is $a{e}^{-j2\pi\tau f}$ which is a band-limited signal in the $f$ variable with spectrum inside the interval $[-\tau_{max}/2,\tau_{max}/2]$. This is a smooth function of $f$ and there exist bounds on its derivatives of any order, (Bernstein’s inequality [19, Th. 6.7]). So, if \mbox{\boldmathA\unboldmath}(f) is the ${M\times K}$ matrix whose $(m,k)$ component is the response of the $m$th sensor to the $k$th DOA, the whole matrix \mbox{\boldmathA\unboldmath}(f) is formed by band-limited functions with spectrum inside $[-\tau_{max}/2,\tau_{max}/2]$, i.e, by smooth functions. In turn, this smoothness of \mbox{\boldmathA\unboldmath}(f) translates into a smooth variation of its span with $f$. To see this point, we must resort to a result in Perturbation Theory, [20, Ch.2, Th. 1.10]. The span of \mbox{\boldmathA\unboldmath}(f) is uniquely represented by an orthogonal projection matrix \mbox{\boldmathP\unboldmath}(f), whose components are analytic functions. 111\mbox{\boldmathP\unboldmath}(f) is the eigenprojector associated with the $K$ eigenvalues of \mbox{\boldmathA\unboldmath}(f)\mbox{\boldmathA\unboldmath}(f)^{H}. Since this last matrix is Hermitian and all its components are entire functions, we may apply Theorem 1.10 (Chapter 2) in [20] to conclude that \mbox{\boldmathP\unboldmath}(f) is analytic. Note than \mbox{\boldmathP\unboldmath}(f) is the projection matrix that is usually obtained through the normal equations if \mbox{\boldmathA\unboldmath}(f) has full column rank, i.e,

[TABLE]

However, since \mbox{\boldmathP\unboldmath}(f) is analytic, it is well defined by continuity even at the isolated frequencies at which the column rank of \mbox{\boldmathA\unboldmath}(f) is smaller than $K$ and, additionally, it is approximable in any band $[f_{A},f_{B}]$ by a polynomial of the form

[TABLE]

where $Q$ is the polynomial order and \mbox{\boldmathG\unboldmath}_{q} are ${M\times M}$ coefficient matrices (Weiertrass theorem, [21, Th. 2.4.1]). Besides, there is an order $Q$ for which the mismatch in (2) is negligible for any possible selection of the DOAs. This is so, because the DOAs can be parameterized using a finite set of variables, say angles of arrival $\theta_{1}$, $\theta_{2},\ldots$, $\theta_{K}$, and the domains of these variables are closed sets (usually finite intervals). The order $Q$ required depends on the specific array geometry and tolerable mismatch and, in principle, must be determined numerically. We may expect the value of $Q$ to increase with ${f_{B}-f_{A}}$, being equal to zero for short bands (narrowband case). [In Sec. VII-B, we assess the selection of $Q$ for a uniform linear array (ULA) formed by 10 sensors.]

In order to see the relevance of the smooth variation of \mbox{\boldmathP\unboldmath}(f), consider a value of $Q$ for which the mismatch of (2) is negligible, and recall the usual binning approach in wideband subspace estimation. In this approach, there is a set of sample covariance matrices \widehat{}\mbox{\boldmathR\unboldmath}(f_{r}), computed at a set of $R$ frequency bins with central frequencies $f_{r}$ lying inside $[f_{A},f_{B}]$, ${r=1,\,2,\ldots}$, $R$, and each of them provides a subspace estimate at its corresponding frequency. More precisely, the span of the $K$ eigenvectors associated with the $K$ largest eigenvalues of each \widehat{}\mbox{\boldmathR\unboldmath}(f_{r}) approximates the span of \mbox{\boldmathP\unboldmath}(f_{r}). But we can describe this approximation in terms of projection matrices. Specifically, if \widehat{}\mbox{\boldmathQ\unboldmath}_{K}(f_{r}) is the $M\times K$ matrix formed by these last eigenvectors and following \widehat{}\mbox{\boldmathQ\unboldmath}_{K}(f_{r})^{H}\widehat{}\mbox{\boldmathQ\unboldmath}_{K}(f_{r})=\mbox{\boldmathI\unboldmath}_{K}, then the matrix

[TABLE]

approximates \mbox{\boldmathP\unboldmath}(f_{r}) component-wise, i.e, [\widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r})]_{m,m^{\prime}} is an estimate of [\mbox{\boldmathP\unboldmath}(f_{r})]_{m,m^{\prime}} for all $m,\,m^{\prime}=1,\,2,\ldots,\,M$. Now if \mbox{\boldmathP\unboldmath}(f) is oversampled, i.e, if $R$ is larger than the number of coefficient matrices in (2), $Q+1$, then we may estimate \mbox{\boldmathG\unboldmath}_{q} in (2) using a simple procedure, such as least squares. This would give a set of estimates \widehat{}\mbox{\boldmathG\unboldmath}_{q} that could then be used to estimate the signal subspace at any frequency in $[f_{A},f_{B}]$, simply by evaluating (2) with \widehat{}\mbox{\boldmathG\unboldmath}_{q} in place of \mbox{\boldmathG\unboldmath}_{q}.

In the sections that follow, we develop this approach in detail. First, we present the signal model in the next section, including the usual binning approach whose output is a set of projection matrix estimates. Then, we present the method in Sec. IV, which is based on weighted least squares, and show that it is asymptotically consistent. We justify this last fact in Sec. V, where we present asymptotic expressions for the bias and for a bound on the RMS error of the projection matrix function estimate. Finally, we assess the method’s quality in Sec. VII numerically, both for the estimation of \mbox{\boldmathP\unboldmath}(f) and for DOA estimation.

III Signal model for wideband subspace estimation

Consider the generic scenario, presented in the previous section, consisting of $M$ sensors and $K$ waves impinging from different directions of arrival. Let ${s}_{bp,k}(t)$ denote the received band-pass signal from the $k$th direction and assume that the receiver’s demodulators operate at a frequency $f_{o}$, so that the lowpass equivalent of ${s}_{bp,k}(t)$ is

[TABLE]

For any array geometry, the effect of the $m$th sensor on ${s}_{k}(t)$ can be described through the convolution with an unknown impulse response $\tilde{}{a}_{m,k}(t)$; i.e, the signal received at the $m$th sensor is

[TABLE]

$\tilde{}{a}_{m,k}(t)$ models models the geometry of the $k$th impinging wave relative to the array and is a time-limited response with support lying inside $[-\tau_{max}/2,\tau_{max}/2]$ for the delay $\tau_{max}$ discussed in the previous section. Besides, since we may transfer any constant factor from $\tilde{}{a}_{m,k}(t)$ to ${s}_{k}(t)$ in the convolution in (5), we may assume that the Fourier transform of $\tilde{}{a}_{m,k}(t)$ is bounded by one; i.e, denoting the Fourier transform of $\tilde{}{a}_{m,k}(t)$ by ${a}_{m,k}(f)$, we have $|{a}_{m,k}(f)|\leq 1$. Note that this model of the sensor array includes, as a particular case, the usual scenario in which the effect of the $m$th sensor on the $k$th signal is just a delay $\tau_{m,k}\in[-\tau_{max}/2,\tau_{max}/2]$, just by taking

[TABLE]

Now, if ${x}_{m}(t)$, ${m=1,\,2,\ldots}$, $M$, denotes the lowpass signal from the $m$th sensor, we have by superposition

[TABLE]

where we have included a set of noise processes ${\eta}_{m}(t)$ that we model as circularly-symmetric, complex-white and of zero mean and variance $\sigma^{2}$.

Next, we proceed to take this last signal model to the frequency domain. For doing this, we assume the following:

• •

The signals ${s}_{k}(t)$ have significant spectral content inside a band $[f_{A},f_{B}]$.

• •

At each sensor, the receiver filters the signals ${x}_{m}(t)$ using $R$ pass-band filters with impulse response ${h}(t;f_{r})$, ${r=1,\,2,\ldots}$, $R$, where $f_{r}$ is the filter’s central frequency. The frequencies $f_{r}$ are distinct, lie in $[f_{A},f_{B}]$, and the filters’ pass bands are $[f_{r}-\Delta f/2,f_{r}+\Delta f/2]$ for a fixed bandwidth $\Delta f$. Besides, the frequencies form an increasing sequence with spacing at least $\Delta f$, ($f_{r}\leq f_{r+1}$, ${r=1,\,2,\ldots}$, $R-1$). For simplicity, we assume unit-energy responses ${h}(t,f_{r})$, so that the noise power at the output of the corresponding filters have the same power, denoted $\sigma^{2}$.

• •

$\Delta f$ is selected to ensure the narrowband condition relative to the possible spectra ${a}_{m,k}(f)$ in the following sense: for any possible DOA and any sensor, it must be ${{a}_{m,k}(f_{1})\approx{a}_{m,k}(f_{2})}$ for any two frequencies following ${|f_{2}-f_{1}|\leq\Delta f}$. Since the responses ${a}_{m,k}(f)$ are time-limited to $[-\tau_{max}/2,\tau_{max}/2]$, from Bernstein’s inequality [22, Th. 11.1.2] we have

[TABLE]

Besides this bound is attained by any sensor with delay $\tau_{max}/2$ as can be readily inferred from (6). Thus, $\Delta f$ must be much smaller than $1/(\pi\tau_{max})$ to ensure a negligible variation of ${a}_{m,k}(f)$.

• •

The scenario is wideband in the sense that $f_{B}-f_{A}>\Delta f$.

From (7), the output of the $r$th filter at the $m$th sensor follows the model

[TABLE]

and, since ${a}_{m,k}(f_{r})$ has a negligible variation in any band of width $\Delta f$, we may approximate the convolution in (9) in the following way

[TABLE]

and, in turn, write (9) as

[TABLE]

Here is where we can see the subspace structure of ${x}_{m}(t)$, given that the signal component on the right-hand side is a linear combination of array response signatures \mbox{\boldmatha\unboldmath}_{m,k}(f_{r}) for any $m$. We may write this model employing the usual notation in narrowband subspace estimation. For this define (${m=1,\,2,\ldots}$, $M$, ${k=1,\,2,\ldots}$, $K$),

[TABLE]

Then, (11) can be written as

[TABLE]

This is a set of $R$ independent narrowband models that allow for the computation of subspace estimates [1, Sec. 9.3]. Specifically, for each output \mbox{\boldmathx\unboldmath}(t;f_{r}), we have the following processing:

If we model the components of \mbox{\boldmaths\unboldmath}(t;f_{r}) as wide-sense stationary processes which are independent of \mbox{\boldmath\epsilon\unboldmath}(t;f_{r}), then the expected covariance of \mbox{\boldmathx\unboldmath}(t,f_{r}), \mathcal{E}\{\mbox{\boldmathx\unboldmath}(t;f_{r})\mbox{\boldmathx\unboldmath}(t;f_{r})^{H}\}, is given by

[TABLE]

where [\mbox{\boldmathR\unboldmath}_{s}(f_{r})]_{k,k^{\prime}} is the cross-power density spectrum of components $k$ and $k^{\prime}$ of \mbox{\boldmaths\unboldmath}(t;f_{r}), (${k,k^{\prime}=1,}\,2,\ldots\,K$). 2. 2.

The ML estimator of \mbox{\boldmathR\unboldmath}(f_{r}) can be readily computed from $N$ samples of \mbox{\boldmathx\unboldmath}(t;f_{r}) and is given by

[TABLE] 3. 3.

Let \widehat{}\mbox{\boldmathQ\unboldmath}_{K}(f_{r}) denote the ${M\times K}$ matrix formed by the $K$ eigenvectors associated with the $K$ largest eigenvalues of \widehat{}\mbox{\boldmathR\unboldmath}(f_{r}). The span of \widehat{}\mbox{\boldmathQ\unboldmath}_{K}(f_{r}) is the ML estimate of the corresponding span of \mbox{\boldmathR\unboldmath}(f_{r}) which, in turn, is the span of \mbox{\boldmathA\unboldmath}(f_{r}) if \mbox{\boldmathR\unboldmath}_{s}(f_{r}) is non-singular.

\widehat{}\mbox{\boldmathQ\unboldmath}_{K}(f_{r}) provides an estimate of \mbox{\boldmathP\unboldmath}(f_{r}) which is

[TABLE]

The set of matrices \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}), ${r=1,\,2,\ldots}$, $R$, is the data from which we proceed to estimate the whole function \mbox{\boldmathP\unboldmath}(f) in $[f_{A},f_{B}]$ in the next section.

IV Estimation of the projection matrix function

Let us assume that an order $Q$ has been selected such that there exist polynomials of the form in (2) with negligible mismatch and that ${Q<R}$. (2) and the narrowband estimates \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}) in (16) allow us to pose a linear regression model for an arbitrary component $(m,m^{\prime})$ of \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}),

[TABLE]

in which [\mbox{\boldmathG\unboldmath}_{q}]_{m,m^{\prime}}, ${q=0,}$ $1,\ldots,$ $Q$, is the set of unknown parameters. The estimates [\widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r})]_{m,m^{\prime}} are independent given that $f_{r+1}-f_{r}\geq\Delta f$, ${r=1,\,2,\ldots}$, $R-1$. Besides, we prove in the next section that \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}) is a consistent estimator of \mbox{\boldmathP\unboldmath}(f_{r}), given that

[TABLE]

where the variance is taken component-wise. All this suggests to pose a weighted least-squares cost function for the estimation of the coefficients [\mbox{\boldmathG\unboldmath}_{q}]_{m,m^{\prime}} with weights $w_{r}$, i.e, the function

[TABLE]

[For a possible set of coefficients $w_{r}$, see (45).] Its minimizer is the estimate

[TABLE]

where (${r,r^{\prime}=1,\,2,\ldots}$, $R$; $q=0,\,1,\ldots,\,Q$)

[TABLE]

Since $B$ is independent of $m$ and $m^{\prime}$, we may write (20) as

[TABLE]

and the estimate of \mbox{\boldmathP\unboldmath}(f) at any $f$ is

[TABLE]

\widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) can be written in terms of \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}) by substituting (23) into (24). We have

[TABLE]

where

[TABLE]

Let us now analyze the quality of the estimate \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) by considering its mean and variance. Since the mismatch of the polynomial approximation is negligible, ${R>Q}$ and the weights $w_{r}$ are positive, we have that (25) holds for \mbox{\boldmathP\unboldmath}(f),

[TABLE]

So, from (18), (25) and this last equation, we have

[TABLE]

Thus, \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) is a consistent estimate of \mbox{\boldmathP\unboldmath}(f).

For any $f$ in $[f_{A},f_{B}]$, \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) is a projection matrix only approximately and it may be necessary to provide an exact projection matrix in some applications. This drawback can be solved by taking the rank-$K$ projection matrix lying closest to \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) as the final signal projection matrix estimate at each frequency; i.e, the final signal projection matrix estimate is

[TABLE]

It can be easily checked that \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f) is just the signal projection matrix of \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f), [i.e, the projection matrix associated with the $K$ largest eigenvalues of \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f)]. Though the computation of \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f) from \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) requires an additional eigenvalue decomposition, we may expect that the number of such decompositions is small in practice, given that \widehat{}\mbox{\boldmathP\unboldmath}(f) can be well approximated in $[f_{A},f_{B}]$ by a $Q$-order polynomial. We assess this point in Secs. VII-E and VII-F numerically.

V Assessment of the first two moments of \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r})

In this section, we present asympotic expressions for the mean and for a bound on the variance of \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}). These expressions prove (18) as a corollary and provide a possible set of coefficients $w_{r}$ in (19). For simplicity, let us suppress the dependency on $f_{r}$ in writing in the rest of this section, i.e, let us write $R$ rather than \mbox{\boldmathR\unboldmath}(f_{r}) and \widehat{}\mbox{\boldmathR\unboldmath} rather than \widehat{}\mbox{\boldmathR\unboldmath}(f_{r}), etc. Let us start by recalling the existing results on the perturbation of the eigenvalues of a sample covariance matrix [23, 24]. For this, write the eigenvalue decomposition of $R$ as

[TABLE]

where $\lambda_{m}$ is the $m$th eigenvalue in decreasing order and \mbox{\boldmathq\unboldmath}_{m} is a corresponding eigenvector, \mbox{\boldmathq\unboldmath}_{m}^{H}\mbox{\boldmathq\unboldmath}_{m^{\prime}}={\delta}_{m-m^{\prime}}, ($m,\,m^{\prime}=1,\ldots,\,M$). Due to the model in (14), we have ${\lambda_{m}=\sigma^{2}}$, $m=K+1,,\ldots,\,M$. For the sample covariance matrix \widehat{}\mbox{\boldmathR\unboldmath} in (15), this same decomposition takes the form

[TABLE]

where “$\widehat{\;\;}$” marks the estimates of the corresponding parameters. Let \widehat{}\mbox{\boldmath\epsilon\unboldmath}_{m} denote the estimation error for \mbox{\boldmathq\unboldmath}_{m}, i.e,

[TABLE]

In order to recall the asymptotic expressions for the first two moments of \widehat{}\mbox{\boldmath\epsilon\unboldmath}_{m}, define first the coefficients

[TABLE]

From Eqs. (12) to (14) in [24], we have the following asymptotic expressions, where “$\cong$” denotes ${o}(1/N)$ equalities (${m,\,m^{\prime}=1,\,2,\ldots}$, $M$):

[TABLE]

[TABLE]

[TABLE]

In App. A , we use these expressions to prove the formula

[TABLE]

where {[\mbox{\boldmathQ\unboldmath}]_{\cdot,m}\equiv\mbox{\boldmathq\unboldmath}_{m}}, (${m=1,\,2,\ldots}$, $M$) and \mbox{\boldmath\Lambda\unboldmath}_{b} is the diagonal matrix

[TABLE]

Since $b_{k,\ell}$ is ${O}(1/N)$, (38) implies that \widehat{}\mbox{\boldmathP\unboldmath}_{0} is asymptotically unbiased provided there is some separation between the signal and noise subspaces, i.e, ${\lambda_{K}-\lambda_{K+1}>0}$. To see this last point, note that all coefficients $b_{k,\ell}$ in (39) depend on one signal eigenvalue and one noise eigenvalue and, therefore, we have

[TABLE]

The variance of the components of \widehat{}\mbox{\boldmathP\unboldmath}_{0} can be computed using the results in [23, 24] through a laborious derivation and the result seems to depend on $P$. However, we may proceed in an indirect way by bounding \|\widehat{}\mbox{\boldmathP\unboldmath}_{0}-\mbox{\boldmathP\unboldmath}\|_{F}^{2}, given that

[TABLE]

In App. B, we prove the following asymptotic bound

[TABLE]

Recalling that ${\lambda_{m}=\sigma^{2}}$ if ${K<m\leq M}$, we may combine the last two inequalities to obtain a bound on the quadratic error for any component $(m,m^{\prime})$,

[TABLE]

This inequality provides a possible choice of coefficients $w_{r}$ in (19), if we eliminate the factors independent of either the eigenvalues $\lambda_{k}$ or $\sigma^{2}$. Specifically, we have the possible coefficients

[TABLE]

where $\widehat{\sigma}$ is a deviation estimate obtained from the covariance matrix \widehat{}\mbox{\boldmathR\unboldmath} of one or more frequency bins. For a single frequency bin, this estimate would be

[TABLE]

Since, from (38), the mean of [\widehat{}\mbox{\boldmathP\unboldmath}_{0}-\mbox{\boldmathP\unboldmath}]_{m,m^{\prime}} is ${O}(1/N)$, (43) is also valid for the variance, i.e,

[TABLE]

The bound in (42) has an interesting interpretation. Since ${\lambda_{K}-\sigma^{2}}$ is the power gap between the signal and noise subspaces, the factor

[TABLE]

is the signal power measured in number of gaps and

[TABLE]

is the noise power again measured in number of gaps. Therefore, the bound in (46) is proportional to the product of these two relative powers and decreases as $1/N$.

VI Application of the proposed method to wideband DOA estimation

In this section, we apply the method in Sec. IV to wideband DOA estimation in a ULA. Let us first particularize the signal model in Sec. III to this type of array, then recall two wideband DOA estimators, IC-MUSIC and MTOPS, and finally adapt these estimators to the method proposed in this paper.

In a ULA, the time-domain response of the $m$th sensor to the $k$th DOA in (5) is

[TABLE]

where

• •

$f_{o}$ is the array’s central frequency, so that the sensor spacing is $c/(2f_{o})$ where $c$ is the propagation speed,

• •

$\gamma_{k}\equiv\sin(\theta_{k})$ and $\theta_{k}$ is the angle of arrival relative to the broadside,

• •

$\tau_{m}$ is the delay associated with the $m$th sensor along the array. If $\tau_{max}$ denotes the array diameter, measured as a delay, then

[TABLE]

The array response \mbox{\boldmathA\unboldmath}(f) in (12) is just the Fourier transform of (49), (${m=1,\,2,\ldots}$, $M$, ${k=1,\,2,\ldots}$, $K$),

[TABLE]

where

[TABLE]

and where we have written \mbox{\boldmathA\unboldmath}(f,\mbox{\boldmath\gamma\unboldmath}) rather than \mbox{\boldmathA\unboldmath}(f) to show the dependency on the parameters $\gamma_{k}$ explicitly. Finally, we may write the model in (14) for \mbox{\boldmathR\unboldmath}(f_{r}) as

[TABLE]

The problem of estimating the angles of arrival $\theta_{k}$ can now be cast as the problem of estimating the parameters $\gamma_{k}$, given that there is a one-to-one relationship between $\theta_{k}$ and $\gamma_{k}$, $\gamma_{k}=\sin(\theta_{k})$.

In this paper, we consider the estimation of the DOA parameters $\gamma_{k}$ by means of the following two methods:

• •

IC-MUSIC, [25, Sec. 4.4.3]. In this estimator, the $K$ DOA estimates are given by the abscissa of the main $K$ local maxima of the pseudo-spectrum

[TABLE]

where

[TABLE]

• •

MTOPS, [26]. This estimator starts by computing a set of $M\times K$ matrices \widehat{}\mbox{\boldmathU\unboldmath}(f_{r}) and another set of $M\times(M-K)$ matrices \widehat{}\mbox{\boldmathV\unboldmath}(f_{r}), ${r=1,\,2,\ldots}$, $R$, whose columns are ortho-normal bases of the signal and noise subspaces respectively; i.e, the columns of \widehat{}\mbox{\boldmathU\unboldmath}(f_{r}) span the subspace associated with \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}) and

[TABLE]

Note that we may view the matrix pair (\widehat{}\mbox{\boldmathU\unboldmath}(f_{r}),\widehat{}\mbox{\boldmathV\unboldmath}(f_{r})) as a function of \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}), given that one such pair can be easily computed from \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}). The MTOPS estimator uses the pseudo-spectrum

[TABLE]

where ($r=2,\,3,\ldots,R$)

[TABLE]

The MTOPS DOA estimates are the smallest $K$ local minima of $\mu(\gamma)$.

These estimators can be easily adapted to the method proposed in this paper, simply by applying them to a set of samples of the projection function estimate \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f); i.e, the matrices \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}), (${r=1,\,2,\ldots}$, $R$), in IC-MUSIC and MTOPS would be replaced by a set of matrices \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f^{\prime}_{r}), (${r=1,\,2,\ldots}$, $R^{\prime}$), where the frequencies $f^{\prime}_{r}$ are regularly spaced in $[f_{A},f_{B}]$ and $R^{\prime}$ is close to $Q+1$. Note that the computation of the matrices \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f^{\prime}_{r}) involves $R^{\prime}$ eigenvalue decompositions. However, this additional complexity is small if $R\gg Q$, given that we may expect $R^{\prime}$ to be close to $Q$.

VII Numerical examples

We have carried out several numerical examples for a 10-sensor ULA following the model in Sec. VI. The sub-sections that follow contain the main results:

• •

Sec. VII-A contains a list of the main parameters in the numerical examples.

• •

In Sec. VII-B, we assess the selection of $Q$ using Chebyshev interpolation in order to upper-bound the mismatch of the polynomial approximation in (2) for several values of $Q$.

• •

In Sec. VII-C, we evaluate the approximation error of one component of \mbox{\boldmathP\unboldmath}(f) using the corresponding component of \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f).

• •

In Sec. VII-D, we evaluate the same error but for the whole matrix \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) using the error measure in (65) and in the RMS sense.

• •

Finally, in Secs. VII-E and VII-F, we combine the proposed method with IC-MUSIC and MTOPS for DOA estimation.

VII-A Main parameters in the numerical examples

The parameters in the numerical examples were the following:

Sensor array. Linear array of ${M=10}$ sensors with half wavelength spacing.

Central frequency. $f_{o}=2.4$ GHz.

Maximum delay along the array.

[TABLE]

DOA parameters. The DOA was parameterized in terms of $\gamma$ rather than $\theta$, where $\gamma=\sin(\theta)$, following the approach in Sec. VI.

Received signals. Linearly-modulated signals of the form

[TABLE]

where

• •

$a_{n}$ are zero-mean, independent complex Gaussian noise samples of variance equal to 1.

• •

${g}(t)$ is a raised-cosine pulse with chip period $T_{ch}\equiv 2.6$ nsec and roll-off factor $\beta\equiv 0.25$.

Sampling period. $T\equiv T_{ch}/2$.

Signals’ bandwidth relative to $f_{o}$. The signals’ two-sided bandwidth $B$ followed ${B/f_{o}=0.2}$, i.e, ${B=0.48}$ GHz. However, in the numerical examples, only the band in which ${g}(t)$ has flat spectrum was used, i.e, the band $[-B_{1}/2,B_{1}/2]$, where ${B_{1}\equiv}{(1-\beta)/T_{ch}}{=0.288}$ GHz. So, the relative bandwidth employed was ${B_{1}/f_{o}}{=0.12}$.

Number of slots. $N_{sl}=50$.

Number of samples per slot. $N=1024$.

Number of frequency bins. The number of frequency bins was fixed to ${R=41}$ and they were equally spaced in $[-B_{1}/2,B_{1}/2]$.

Directions of arrival (DOAs). Two cases have been assessed:

• •

Three DOAs given by

[TABLE]

• •

Two DOAs of the form

[TABLE]

where the increment $\Delta\gamma$ is a simulation parameter.

Signal-to-noise ratio (SNR). The SNRs in the numerical examples are equal to the total signal energy at frequency $f_{o}$ divided by the corresponding noise energy.

Least squares weights $w_{r}$. Equal to 1.

Generation method for numerical trials. The numerical trials were generated using the frequency domain model in (13) for reducing the simulation time. However, the results were validated by generating these trials in the time domain and checking that the results are consistent if $N\rightarrow\infty$.

Number of Monte Carlo trials. 2200.

VII-B Polynomial order $Q$ versus mismatch for a ten-sensor uniform linear array

We have computed an upper bound for the approximation error in (2) by means of Chebyshev interpolation applied to \mbox{\boldmathP\unboldmath}(f) in range ${f\in[f_{A},f_{B}]}$ with ${Q+1}$ nodes, where the array response matrix is given by (51), [27, Ch. 6]. This interpolation scheme delivers a polynomial approximation in the $f$ variable of the form

[TABLE]

for $f$ in $[f_{A},f_{B}]$ and fixed $Q$, $\gamma_{1}$, $\ldots$, $\gamma_{K}$. Specifically, we have evaluated the error in (62) for any value of the variables involved. Table I shows the maximum error in (62) for several values of $Q$ and $K$, i.e the error measure

[TABLE]

Note that, since the components of \mbox{\boldmathP\unboldmath}(f) are bounded by one, this interpolation scheme provides significant accuracy even for small values of $Q$. Table II shows the error measure

[TABLE]

where, again, the supremum is taken over all variables involved. This is the maximum error that would be incurred if \mbox{\boldmathP\unboldmath}(f) were replaced by its polynomial approximation in the computation of \mbox{\boldmathP\unboldmath}(f)\mbox{\boldmathx\unboldmath} for unit-norm vectors $x$. The conclusion is the same.

VII-C Approximation of \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) and \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f) to \mbox{\boldmathP\unboldmath}(f)

In this section, we assess qualitatively the error in approximating the true projection matrix \mbox{\boldmathP\unboldmath}(f) using either \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f), \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) or \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f) for a single component of these matrices and the DOAs in (60). The objective of this assessment is to show the improvement that can be achieved graphically.

Fig. 1 shows the real and imaginary parts of component [\mbox{\boldmathP\unboldmath}(f)]_{1,10} and the state-of-the-art estimate [\widehat{}\mbox{\boldmathP\unboldmath}_{0}(f)]_{1,10} in one realization of the numerical example for ${\textrm{SNR}=8~{}\textrm{dB}}$. The smooth thick curves are the components of [\mbox{\boldmathP\unboldmath}(f)]_{1,10} and the noisy thin curves the corresponding components of [\widehat{}\mbox{\boldmathP\unboldmath}_{0}(f)]_{1,10}. Note that [\widehat{}\mbox{\boldmathP\unboldmath}_{0}(f)]_{1,10} approximates [\mbox{\boldmathP\unboldmath}(f)]_{1,10} but with some error, fundamentally due to the variation of the received signals’ sample spectra. Though this figure only shows one component of \mbox{\boldmathP\unboldmath}(f) and \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f), the trends also hold for the rest of components, i.e, the whole matrix \mbox{\boldmathP\unboldmath}(f) is a smooth function of $f$ and {\mbox{\boldmathP\unboldmath}(f)\approx\widehat{}\mbox{\boldmathP\unboldmath}_{0}(f)}.

Fig. 2 shows the result of fitting a polynomial of order ${Q=3}$ to the real part of ${R=41}$ equally-spaced samples of \mathrm{Re}\{[\widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r})]_{1,10}\}. The continuous curve is the true value [\mbox{\boldmathP\unboldmath}(f)]_{1,10} and the dashed curve the fitted value [\widehat{}\mbox{\boldmathP\unboldmath}_{1}(f)]_{1,10}. Note that the fitted value is a significantly better estimate of [\mbox{\boldmathP\unboldmath}(f)]_{1,10} along the frequency band than the initial estimates (dots).

The final correction in (30) for obtaining \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f) from \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) produces a slight variation, that can be readily seen in Fig. 3 for the real part of component $(1,10)$. This figure shows the error in approximating \mbox{\boldmathP\unboldmath}(f) using either \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) or \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f) for component $(1,10)$. Note that the curve is smooth for \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) and \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f) and that the approximation error is small in both cases.

VII-D RMS approximation error of \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) versus the number of sample covariance matrices $R$

Fig. 4 shows the approximation error for the whole projection matrix in the example of the previous sub-section, where the error norm is

[TABLE]

Note that, except for $Q=0$, \widehat{}\mbox{\boldmathP\unboldmath}_{1}(f) outperforms \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f) by a significant margin that can reach $17~{}\textrm{dB}$ for a high number of bins $R$.

Fig. 5 shows the curves in Fig. 4 for ${Q\neq 2}$ minus the curve for ${Q=2}$ in dBs. This figure allows us to see what value of $Q$ performs best versus the number of bins $R$. As can be readily seen, ${Q=1}$ is the best choice up to ${R=10}$ (value below zero in “${Q=1}$” curve), while ${Q=2}$ is the best choice between $R=11$ and ${R=250}$, and ${Q=3}$ is the best choice for ${R>250}$.

VII-E Improvement in DOA separation

In order to assess the effect of the proposed method on the separation of DOA estimates, we have evaluated the RMS error in the estimation of $\gamma_{1}$ for varying $\Delta\gamma$ in (61) using several variants of IC-MUSIC and MTOPS for ${\textrm{SNR}=-10~{}\textrm{dB}}$. The variants were the following:

• •

Standard IC-MUSIC. Standard IC-MUSIC estimator using the pseudo-spectrum in (54).

• •

41-bin IC-MUSIC. Standard IC-MUSIC estimator but with projection matrices \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}) replaced with their estimates \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f_{r}).

• •

5-bin IC-MUSIC. IC-MUSIC estimator applied to $R^{\prime}=5$ projection matrices \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f^{\prime}_{r}) as proposed at the end of Sec. VI. The frequencies $f^{\prime}_{r}$ formed a regular grid covering the band $[f_{o}-B_{1}/2,f_{o}+B_{1}/2]$ and the IC-MUSIC pseudo-spectrum was (54) but with frequencies $f^{\prime}_{r}$, ${r=1,\,2,\ldots}$, $R^{\prime}$.

• •

1-bin IC-MUSIC. The same estimator but computed from the single projection matrix \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f_{o}), i.e, with $f^{\prime}_{1}=f_{o}$ and $R^{\prime}=1$.

And the MTOPS estimators were:

• •

Standard MTOPS. MTOPS estimator computed from the pseudo-spectrum in (57). Since this pseudo-spectrum has numerous local peaks, the peak lying closest to 1-bin IC-MUSIC was selected as estimate. (See [26] for comments on this drawback of MTOPS.)

• •

41-bin MTOPS. Standard MTOPS estimator but with projection matrices \widehat{}\mbox{\boldmathP\unboldmath}_{0}(f_{r}) replaced with their estimates \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f_{r}).

• •

5-bin MTOPS. MTOPS estimator applied to ${R^{\prime}=5}$ projection matrices \widehat{}\mbox{\boldmathP\unboldmath}_{2}(f^{\prime}_{r}), (Sec. VI), where the frequencies $f^{\prime}_{r}$ formed a regular grid covering the band $[f_{o}-B_{1}/2,f_{o}+B_{1}/2]$.

Fig 6(a) shows the performance of 41-bin IC-MUSIC versus that of standard IC-MUSIC. Note that the proposed method improves the smallest $\Delta\gamma$ for which the two DOAs are separable. With the standard estimator, we have $\Delta\gamma=0.15$, while the adapted method with $Q=1$ provides $\Delta\gamma=0.11$ roughly. Larger values of $Q$ imply a larger threshold for $\Delta\gamma$. Fig. 6(b) shows the same comparison for 5-bin IC-MUSIC and the behavior is the same, but with a more noticeable degradation with increasing $Q$.

Figs. 6(c), 7(a) and 7(b) show the same comparison for 5-bin IC-MUSIC, 41-bin MTOPS and 5-bin MTOPS, respectively, and the conclusions that can be drawn are similar.

VII-F Performance improvement in SNR threshold

We have repeated the numerical example in the previous sub-section, but for the three DOAs in (60), and assessed the RMS error of $\gamma_{1}$ for varying SNR. Fig. 8(a) shows the performance of 41-bin IC-MUSIC. Note that the behavior is similar to the one in the previous sub-section, i.e, the method provides an improvement of roughly 7 dBs in the SNR threshold. Additionally, we can see that the mismatch for the polynomial approximation can be preceived in the high RMS error for high SNRs and $Q=0$. Figs 8(b), 8(c), 9(a), 9(b) show the result of the same assessment for the remaining estimators in the previous section, 5-bin and 1-bin IC-MUSIC, and 41-bin and 5-bin MTOPS, respectively, with similar conclusions.

VIII Conclusions

We have presented a method for characterizing the signal subspace through a polynomial approximation of the signal projection matrix, valid in a given frequency band. Fundamentally, the method consists of fitting a polynomial to a set of sample signal projection matrices, obtained through the usual binning approach in wideband subspace estimation. The resulting polynomial provides an approximate signal projection matrix at any frequency in the band, which can then be used to improve the quality of wideband DOA estimators such as IC-MUSIC and MTOPS. We have presented asymptotic bounds for the bias and RMS error of the polynomial estimate and we have assessed its performance in several numerical examples.

Appendix A First- and second-order moments of sample covariance matrix eigenvectors

Let us first compute the first- and second-order moments of \widehat{}\mbox{\boldmathq\unboldmath}_{m}. From (35), we have

[TABLE]

where we define the coefficients

[TABLE]

From (36) and (66), we have

[TABLE]

If $m=m^{\prime}$ this expression reduces to

[TABLE]

Adding up (69) for ${m=1,\,2,\ldots}$, $K$, we obtain a formula for \mathcal{E}\{\widehat{}\mbox{\boldmathP\unboldmath}_{0}\}, noting that \mbox{\boldmathP\unboldmath}=\sum_{k=1}^{K}\mbox{\boldmathq\unboldmath}_{k}\mbox{\boldmathq\unboldmath}_{k}^{H}:

[TABLE]

And doing the same for $m=K+1,\,\ldots,\,M$, we obtain

[TABLE]

The last two formulas involve coefficients $b_{k,\ell}$ in which both $k$ and $\ell$ are either smaller than ${K+1}$ or larger than $K$, i.e, coefficients computed from pairs of eigenvalues associated to either the signal or noise subspace. Such coefficients $b_{k,\ell}$ can be arbitrarily large, given that the only condition on the eigenvalues $\lambda_{m}$ is that there is a significant gap $\lambda_{K}-\lambda_{K+1}$ between the signal and noise eigenvalues. We solve this drawback by writing \mathcal{E}\{\widehat{}\mbox{\boldmathP\unboldmath}_{0}\} in terms of \mathcal{E}\{(\mbox{\boldmathI\unboldmath}_{M}-\mbox{\boldmathP\unboldmath})\widehat{}\mbox{\boldmathP\unboldmath}_{0}\} and \mathcal{E}\{\mbox{\boldmathP\unboldmath}(\mbox{\boldmathI\unboldmath}_{M}-\widehat{}\mbox{\boldmathP\unboldmath}_{0})\}. We have from (A),

[TABLE]

and from (71),

[TABLE]

Finally, combining the last two equations, we have

[TABLE]

The matrix form of this expression is (38).

Appendix B Derivation of bound on expected quadratical error

The eigenvalue decomposition of $R$ in (31) can be written as

[TABLE]

where $Q$ is unitary, [\mbox{\boldmathQ\unboldmath}]_{m}\equiv\mbox{\boldmathq\unboldmath}_{m}, and $\Lambda$ is a diagonal matrix with components [\mbox{\boldmath\Lambda\unboldmath}]_{m,m^{\prime}}={\delta}_{m-m^{\prime}}\lambda_{m}, ${m,m^{\prime}=1,\,2,\ldots}$, $M$. The sample covariance matrix \widehat{}\mbox{\boldmathR\unboldmath} can be described as the average of $N$ independent, complex normal ${M\times 1}$ vectors \mbox{\boldmathx\unboldmath}_{n}, of zero mean and covariance $R$,

[TABLE]

Next, consider the vectors \mbox{\boldmaths\unboldmath}_{n}\equiv\mbox{\boldmathQ\unboldmath}^{H}\mbox{\boldmathx\unboldmath}_{n} which are also independent, complex normal and of zero mean, but of covariance matrix $\Lambda$. The sample covariance matrix of these vectors is

[TABLE]

Given a realization \widehat{}\mbox{\boldmathR\unboldmath} and its corresponding \widehat{}\mbox{\boldmathC\unboldmath}, let us bound the error in approximating $P$ using \widehat{}\mbox{\boldmathP\unboldmath}_{0} by resorting to a perturbation theory result in [28]. First, define the following measure for the dissimilarity between the spans of $P$ and \widehat{}\mbox{\boldmathP\unboldmath}_{0}:

[TABLE]

where \widehat{}\mbox{\boldmathC\unboldmath}_{sn} is the block formed by the intersection of the first $K$ rows and last $M-K$ columns of \widehat{}\mbox{\boldmathC\unboldmath}, i.e, \widehat{}\mbox{\boldmathC\unboldmath}_{sn}\equiv[\widehat{}\mbox{\boldmathC\unboldmath}]_{1:K,K+1:M}. Combining theorems 2.1 and 3.1 of [28], we have that if $\widehat{\eta}<1/4$ then

[TABLE]

(See also comments on page 232 of [28].)

In order to turn (79) into an asymptotic inequality, let us first compute the first two moments of the components of \widehat{}\mbox{\boldmathC\unboldmath}_{sn} and then the mean and variance of $\widehat{\eta}$. For simplicity, let $c_{k,\ell}$ and $s_{n,k}$ denote [\widehat{}\mbox{\boldmathC\unboldmath}]_{k,\ell} and [\mbox{\boldmaths\unboldmath}_{n}]_{k} respectively. For any indices $k$, $\ell$, $p$ and $q$, lying between 1 and $M$ and following $k\neq\ell$ and $p\neq q$, we have:

• •

$\mathcal{E}\{c_{k,\ell}\}=0$ given that $k\neq\ell$ and

[TABLE]

• •

The formula

[TABLE]

To prove this result, recall the formula for the expectation of the product of four complex normal random variables $a_{k}$, [29]:

[TABLE]

The proof is the following. We have

[TABLE]

In this parenthesis, we have:

• –

The first term is zero because \mathcal{E}\{s_{n,k}s_{n,\ell}^{*}\}=[\mbox{\boldmath\Lambda\unboldmath}]_{k,\ell}=0.

• –

If $n\neq n^{\prime}$ the second term is zero because \mbox{\boldmaths\unboldmath}_{n} and \mbox{\boldmaths\unboldmath}_{n^{\prime}} are independent and \mathcal{E}\{\mbox{\boldmaths\unboldmath}_{n}\}=0. If $n=n^{\prime}$ this term is zero because \mathcal{E}\{\mbox{\boldmaths\unboldmath}_{n}\mbox{\boldmaths\unboldmath}_{n}^{T}\}=\mbox{\boldmath0\unboldmath}_{M}.

• –

The third term is zero if ${n\neq n^{\prime}}$ because \mbox{\boldmaths\unboldmath}_{n} is independent of \mbox{\boldmaths\unboldmath}_{n^{\prime}} and \mathcal{E}\{\mbox{\boldmaths\unboldmath}_{n}\}=0. And, if $n=n^{\prime}$, then it is also zero if $k\neq q$ or $\ell\neq p$, because \mathcal{E}\{\mbox{\boldmaths\unboldmath}_{n}\mbox{\boldmaths\unboldmath}_{n}^{H}\}=\mbox{\boldmath\Lambda\unboldmath} is a diagonal matrix. Thus, we have that the third term is equal to

[TABLE]

• –

The fourth term is zero because \mathcal{E}\{\mbox{\boldmaths\unboldmath}_{n}\}=0.

So, in summary, we have

[TABLE]

• •

We also have

[TABLE]

The proof is the following. We have

[TABLE]

Comparing this expression with the second line of (83), we can readily see that

[TABLE]

Therefore, from (85), we obtain (86).

Let us now compute the mean of $\widehat{\eta}$. From (78) and (86), we have

[TABLE]

Regarding the second-order moment, it follows the formula:

[TABLE]

Expanding the summand using (82), we have

[TABLE]

Note that this is a sum of products consisting of factors of the form in either (81) or (86). As a consequence, all these products are ${O}(1/N^{2})$ and we have

[TABLE]

Finally, let us derive the asymptotic inequality. Start by decomposing the expectation of \|\widehat{}\mbox{\boldmathP\unboldmath}_{0}-\mbox{\boldmathP\unboldmath}\|_{F}^{2} by conditioning on $\widehat{\eta}^{2}<1/4$:

[TABLE]

From (79), we have that the first term follows

[TABLE]

Regarding the second term, the expectations inside the parenthesis are bounded, because projection matrices have components bounded by one. Besides, we may apply Markov’s inequality and use (92) to obtain

[TABLE]

So, we have that the whole second term in (B) is ${O}(1/N^{2})$ and, therefore, is ${o}(1/N)$. In summary, recalling (B), we obtain the asymptotic inequality

[TABLE]